subroutine iztrsb(lamda,xlat,iprofl,pteran,pcloud,
     1 iofset,lprint,ezgr,ssgr,tgr,sbgr,ezcl,sscl,tcl,sbcl)
c
c***********************************************************************
c
c     iztrsb
c
c     september, 2001 by charlie wellemeyer of ssai
c
c     purpose
c       performs pressure interpolations of the table parameters i0
c       tr and sb, which are used by the calling routines to calculate
c       table radiances or n-values. uses lagrangian interpolation
c       except for high clouds for which a linear extrapolation is used. 
c
c     modification
c       Zack Zhang, SSAI, Mar. 1, 2002
c       Add checkup statement in function pintrp to avoid negative 
c       value from Lagrange interpolation
c     variables
c       name      type  i/o   description
c       ----      ----  ---   -----------
c      arguments
c       lamda     i*4    i    wavelength index
c       xlat      r*4    i    latitude
c       iprofl    i*4    i    profile index
c       pteran    r*4    i    terrain pressure
c       pcloud    r*4    i    cloud pressure
c       iofset    i*4    i    offset pointer.
c       lprint    l*4    i    if = .true., print
c       ezgr      r*4    o    total atmospheric term for ground
c       ssgr      r*4    o    single scattering atmospheric term for ground
c       tgr       r*4    o    radiance transmitted from ground
c       sbgr      r*4    o    atmospheric reflectance at ground
c       ezcl      r*4    o    total atmospheric term for cloud
c       sscl      r*4    o    single scattering total atmos term for cloud
c       tcl       r*4    o    radiance transmitted from cloud
c       sbcl      r*4    o    atmospheric reflectance at cloud
c
c      internal
c       indx1     i*4         table index
c       indx2     i*4         coefficient index
c       iptab1(4) i*4         pressure offset for table index
c       iptab2(4) i*4         pressure offset for coefficient index
c       lamtb1(9) i*4         wavelingth offset for table index
c       lamtb2(9) i*4         wavelength offset for coefficient index
c
c***********************************************************************
c
      implicit none
c
c     -- input parameters
      integer lamda,iprofl,iofset
      integer lamda1
      logical lprint(30)
      real xlat,pteran,pcloud
c     -- internal parameters
      integer iptab1(4),lamtb1(9),iptab2(4),lamtb2(9)
      integer layer1(12), layer2(12), ip, indx1, indx2
      real i0ofp(4), trofp(4), sbofp(4), i0, tr, sb, lgpcld,
     1    iss,issofp(4)
c     -- output parameters
      real ezgr,ssgr,tgr,sbgr,ezcl,sscl,tcl,sbcl
      real pintrp
c
      include 'pterp.com'
c
c     -- data statements of nodal indecese for pressure and
c     -- wavelength in 1) angle dependent and 2) angle
c     -- independent portions of the tables
c
      data iptab1/1,1891,3781,5671/
      data iptab2/1,190,379,568/
      data lamtb1/0,210,420,630,840,1050,1260,1470,1680/
      data lamtb2/0,21,42,63,84,105,126,147,168/
c
      if (lprint(18)) write (6,1000) lamda,
     1 iprofl+1,pteran,pcloud,iofset
c
      do ip = 1,4
c        -- compute offset addresses for each pressure node
         indx1=iptab1(ip)+lamtb1(lamda-4)+iprofl*10+iofset
         indx2=iptab2(ip)+lamtb2(lamda-4)+iprofl
c
c        -- perform table interpolation
c
         lamda1=lamda
         if (lamda .eq. 13) lamda1=lamda-2
         call interp(lamda1,indx1,indx2,lprint,i0,iss,tr,sb)
c     if(lamda.eq.10) then
c       print*,'iztrsb',ip,lamda,iprofl,iofset
c       print*,i0,iss,tr,sb
c     endif
c
c        -- store parameters for pressure interpolation
c
         i0ofp(ip) = i0
         issofp(ip) = iss
         trofp(ip) = tr
         sbofp(ip) = sb
         if (lprint(18)) write (6,1100) ip, indx1, indx2, sb
         if (lprint(18)) write (6,1200) ip, i0, iss, tr
      enddo
c
c     -- interpolate log of table parameters in log pressure
c
         ezgr = pintrp(cofpgr,i0ofp)
         ssgr = pintrp(cofpgr,issofp)
         tgr = pintrp(cofpgr,trofp)
         sbgr = pintrp(cofpgr,sbofp)
c     if(lamda.eq.10) then
c       print*,'ezgr,tgr,sbgr,cofpgr'
c       print*,ezgr,tgr,sbgr,cofpgr
c       print*,i0ofp
c     endif
         if(pcloud.ge.0.4) then
           ezcl = pintrp(cofpcl,i0ofp)
           sscl = pintrp(cofpcl,issofp)
           tcl = pintrp(cofpcl,trofp)
           sbcl = pintrp(cofpcl,sbofp)
c
c  use linear extrapolation for clouds at pressures less than 0.4
c
         else
           lgpcld=alog10(pcloud)
           ezcl=(lgpcld-logpr(2))/(logpr(3)-logpr(2))*
     1                 (i0ofp(3)-i0ofp(2))+i0ofp(2)
           sscl=(lgpcld-logpr(2))/(logpr(3)-logpr(2))*
     1                 (issofp(3)-issofp(2))+issofp(2)
           tcl=(lgpcld-logpr(2))/(logpr(3)-logpr(2))*
     1                 (trofp(3)-trofp(2))+trofp(2)
           if(tcl .lt. 0) tcl=(trofp(3)+trofp(2))*0.5
           sbcl=(lgpcld-logpr(2))/(logpr(3)-logpr(2))*
     1                 (sbofp(3)-sbofp(2))+sbofp(2)
         endif
c
      if (lprint(18)) write (6,1300) ezgr,ssgr,tgr,sbgr,ezcl,sscl,
     1     tcl,sbcl
c
      return
c
c  format statements
c
1000  format ('Subroutine iztrsb'/'Input:   lamda  = ',i8,
     1 ' iprofl = ',i8,' pteran = ',f8.4,' pcloud  = ',f8.4,
     2 ' iofset = ',i8)
1100  format ('Intern:  ip     = ',i8,' indx1  = ',i8,' indx2  = ',i8,
     1     ' sb = ',f8.3)
1200  format ('         ip     = ',i8,' ezcor    = ',f8.6,
     1 ' sscor    = ',f8.6,' tcor    = ',f8.6)
1300  format ('Output:  ezgr  = ',f8.5,' ssgr  = ',f8.5,' tgr   = ',
     1   f8.5,' sbgr   = ',f8.5,/,9x,'ezcl  = ',f8.5,' tcl  = ',
     2   f8.5,' sscl  = ',f8.5,' sbcl  = ',f8.5)
c
      end
c
      function pintrp(cofs,ys)
c
      real pintrp, cofs(4), ys(4)
c
      pintrp = 0.0
      do 100 i=1,4
         pintrp= pintrp + ys(i) * cofs(i)
  100 continue
      if(pintrp .lt. 0.0) then
      pintrp = 0.0
      do 101 i=1,4
         pintrp= pintrp + ys(i) * 0.25
  101 continue
      endif
c
      end