      SUBROUTINE XBTbiaslijing(depth,temp,maxlevel,maxcalc,
     * itsprobe,ixbttype,iyear0,numsec,numlevels,bmiss,
     * xswitch,iplat)

C     XBT BIAS CORRECTION FROM LIJING ET AL. 2011
C     SO FAR THE ONLY PROBE TYPES COVERED ARE
C      2 (xbt4 - 208) [Also 202, 203, 205, 206, 208
C                      209, 211, 212, 215, 216 ]
C     42 (xbt7 - 207) [Also 201, 213, 214, 235, 236,
C                      217, 232]
C    212 (xbt6t - 216)
C    222 (xbt7t - 217) [Also 232 TSK DB]
C    ASSUMPTIONS: 
C     UNKNOWN PROBE TYPE SHALLOW: SAME AS xbt4
C     UNKNOWN PROBE TYPE DEEP: SAME AS xbt7
C     T5 NOT CORRECTED
C     KNOWN PROBES NOT COVERED ABOVE USE
C     SHALLOW/DEEP SAME AS UNKNOWN
C     SHALLOW/DEEP CRITERA SWITCH AT 550 M.

      parameter (ishallowprobe=208,ideepprobe=207)
      parameter (ideepprobetsk=217)
      parameter (icmax=10,maxyear=3000)
      parameter (maxlevelloc=6000)

      character*13 firsttime
      character*30 ctext,cprobe
      character*80 text,filename

      dimension icomma(icmax,2)
      dimension iyearcount(maxyear,0:4)
      dimension xoutput(maxyear,5,4)
      dimension xoutsd(maxyear,5,3)
      dimension iprobetype(maxyear,5)
      dimension depth(maxlevel),temp(maxlevel,maxcalc)
      dimension t0(maxlevelloc)

      common /eachparm/ numeach(6),numtot2,numtot3,ntax,numsecset

      save firsttime
      save xoutput,xoutsd,iprobetype,iyearcount,
     * iyearmin,iyearmax

      if ( firsttime .ne. 'XBTbiaslijing' ) then
       firsttime='XBTbiaslijing'

c     read in bias corrections

       call extraname('sys.inf'//CHAR(0),
     *  'lijing_XBT_corrs.csv'//CHAR(0),
     *  filename)

       open(9,file=filename,status='old')

       read(9,'(a)') text
       iyearmin=maxyear
       iyearmax=0
       itype0=0

       do 50 n=1,maxyear

        call clearstring(text,80)
        read(9,'(a)',end=4) text
        call findlastchar(nsp,text,80)
c       write(6,*) text(1:nsp)

        call commafindx(nsp,icmax,ncomma,icomma,text)

c     value
        ivaluefound=0
        if ( text(icomma(4,1):icomma(4,1)+icomma(4,2)-1)
     *   .ne. 'NaN' ) then 
         read(text(icomma(4,1):icomma(4,1)+icomma(4,2)-1),*)
     *   xvalue
         ivaluefound=1
        endif
        if ( ivaluefound .eq. 0 ) goto 50

c     standard deviation
        if ( text(icomma(5,1):icomma(5,1)+icomma(5,2)-1)
     *   .ne. 'NaN' ) then 
         read(text(icomma(5,1):icomma(5,1)+icomma(5,2)-1),*)
     *    xsd
        endif

        call clearstring(ctext,30)
        ctext(1:icomma(1,2))=
     *   text(icomma(1,1):icomma(1,1)+icomma(1,2)-1)
        nspx=icomma(1,2)
        itype=0
        if ( nspx .eq. 1 .and. ctext(1:1) .eq.
     *   'a' ) then
         itype=1
        elseif ( nspx .eq. 1 .and. ctext(1:1) .eq.
     *   'b' ) then
         itype=2
        elseif ( (nspx .eq. 6 .and. ctext(1:6) .eq.
     *   'offset') .or. (nspx .eq. 1 .and.
     *   ctext(1:1) .eq. 'c' ) ) then
         itype=3
        elseif ( (nspx .eq. 12 .and. ctext(1:12) .eq.
     *   'thermal bias') .or. (nspx .eq. 7 .and.
     *    ctext(1:7) .eq. 'delta T')  ) then
         itype=4
        endif

c       if ( itype .ne. itype0 )
c    *   call arrayinit(iyearcount,1,maxyear)
        itype0=itype


c     Year
        read(text(icomma(2,1):icomma(2,1)+icomma(2,2)-1),*
     *   ) xyear
        nyear=xyear
        if ( iyearmax .lt. nyear) iyearmax=nyear
        if ( iyearmin .gt. nyear) iyearmin=nyear
c       if ( nyear .eq. 1976 ) write(6,*) 'n',itype,iyearcount(nyear)

        iyearcount(nyear,itype)=iyearcount(nyear,itype)+1
c       iyearcount(nyear+1)=iyearcount(nyear+1)+1
c       iyearcount(nyear+2)=iyearcount(nyear+2)+1
        iyc=iyearcount(nyear,itype)

c     Probe type
        cprobe=text(icomma(3,1):icomma(3,1)+icomma(3,2)-1)
c       getxcowl is in subroutine XBTbiascowlX.f
        call getxcowl(icomma(3,2),cprobe,iwodprobe)

c       write(6,*) nyear,xvalue,xsd,iwodprobe
        if ( itype .gt. 0 ) then
c        do 35 ny=nyear,nyear+2
         do 35 ny=nyear,nyear
          xoutput(ny,iyc,itype)=xvalue 
          xoutsd(ny,iyc,itype)=xsd
          iprobetype(ny,iyc)=iwodprobe
35       continue
        endif

50     continue 
4      continue

       close(9)
       do 555 n2=iyearmin,iyearmax
        iyearcount(n2,0)=iyearcount(n2,1)
        do 556 n3=2,3
556      if ( iyearcount(n2,n3) .gt. iyearcount(n2,0) )
     *        iyearcount(n2,0)=iyearcount(n2,n3)
555    continue

      endif

c     Determine year and probe type

      if ( itsprobe .ne. 2 ) return

c      Determine shallow or deep

      if ( numlevels .le. maxlevel ) then
       dpt=depth(numlevels)
      else
       call overmax(maxlevel,numsec,numlevels,0,ic,np)
       dpt=temp(ic,np)
      endif

      call xbtshallowdeep(ishallow,ixbttype,dpt,550.,iplat)
c     write(6,*) 'shallow',ishallow
c     call xbtshallowdeep(ishallow,ixbttype,xdepth)

c     Do not adjust if ishallow=-1 (different probe types)

      if ( ishallow .eq. -1 ) return

c      Check if this is a TSK

      itsk=istsk(ixbttype)

      ifound=0
      xa=0.
      xb=0.
      xc=0.
      tb=0.

      iyear=iyear0
      if ( iyear .lt. iyearmin ) iyear=iyearmin
      if ( iyear .gt. iyearmax ) iyear=iyearmax

      do 55 n=1,iyearcount(iyear,0)

       iput=0
       if ( ixbttype .eq. iprobetype(iyear,n) ) then
        ifound=2
        iput=1
c    Special case: TSK DB same eq. for TSK T7
       elseif ( iprobetype(iyear,n) .eq. ideepprobetsk
     * .and. itsk .gt. 0 .and. ixbttype .eq. 232 ) then
        ifound=2
        iput=1
       elseif ( ifound .lt. 2 ) then 

        iput=1
        if ( ishallow .eq. 1 .and. iprobetype(iyear,n)
     *       .eq. ishallowprobe ) then
         ifound=2
        elseif ( ishallow .eq. 0 .and. iprobetype(iyear,n)
     *   .eq. ideepprobe ) then
         ifound=2
        elseif ( ifound .gt. 0 ) then
         iput=0
        endif

        if ( ifound .eq. 0 ) ifound=1
         
       endif
      
c      write(6,*) 'n',n,iput,xoutput(iyear,n,1),
c    *  xoutput(iyear,n,2),xoutput(iyear,n,3),
c    * iyear,ishallow,
c    * ixbttype,iprobetype(iyear,n)
       if ( iput .eq. 1 ) then
        xa=xoutput(iyear,n,1)
        xb=xoutput(iyear,n,2)
        xc=xoutput(iyear,n,3)
        tb=xoutput(iyear,n,4)
       endif

55    continue

c     make depth/temperature corrections

      if ( xswitch .eq. 1. ) then
       ihanawa=1
      else
       ihanawa=0
      endif

c     write(6,*) 't0',ihanawa,xa,xb,xc,tb
      call timefromdepth(numlevels,t0,maxlevel,ihanawa,depth)
c     write(6,*) 'tfromd',numlevels,t0(numlevels)
      call droprate_apply(xa,xb,numlevels,depth,maxlevel,t0)

      do 80 n=1,numlevels

       if ( n .le. maxlevel ) then
        depth(n)=depth(n)-xc
        temp(n,1)=temp(n,1)-tb
       else
        call overmax(maxlevel,numsec,n,0,ic,np)
        temp(ic,np)=temp(ic,np)-xc
        temp(ic,np+1)=temp(ic,np+1)-tb
       endif

80    continue

      return
      end
