subroutine szacorr (error,errorlat,Tb,sza,xjday,xlat,Tbcorr)

c     programmer:      Bob Joyce
c     date created:    3/2/2000
c     purpose:         Subroutine takes the following inputs
c          1. uncorrected Tb
c          2. satellite zenith angle
c          3. Julian day
c          4. latitude of IR Tb pixel
c
c                      Subroutine outputs
c          1. Seasonal, Latitudinal, Tb dependent
c                            Zenith angle IR Tb correction

      parameter (ZAMAX=87.0, ISTYR=1999, J5MAX=240, XMISSING=-9999,
     2 ZABINMX=64,TMPBINMX=170,SZAANGMN=26.,SZATMPMN=160.,
     3 YEARDYS=365, ISEASMX=8, ZERO=0, ONE=1, XLATMX=60, XLATMN=-60)

      real XMISSING, seasjday(ISEASMX)
      real error(ZABINMX,TMPBINMX), errorlat(J5MAX,ISEASMX), xjday

      data seasjday/57.5,104.5,148.5,206.5,268.5,326.5,395,422.5/

      XLATBININT=0.5

c     add yearly day total if after feb 26 (when model began)

      if (xjday.lt.seasjday(ONE)) xjday = xjday + YEARDYS

c     Find latitude bin for target

      if (xlat.gt.XLATMX) then
       j5x = ONE
      elseif (xlat.gt.XLATMN.and.xlat.lt.XLATMX) then
       j5x = (XLATMX - xlat)/XLATBININT + ONE
      else
       j5x = J5MAX
      endif

c   Check to see if current day is exactly seasonal model

      do 35 icday = 1, ISEASMX-1
       if (xjday.eq.seasjday(icday)) then
        xerrorlat = errorlat(j5x,icday)
        goto 47
       else
        if (xjday.gt.seasjday(icday).and.xjday.lt.seasjday(icday+1))then

c   Find day difference from current day and early and late seasonal tables

         diff1 = xjday - seasjday(icday)
         diff2 = seasjday(icday+1) - xjday
         xlnth = seasjday(icday+1) - seasjday(icday)
         w1    = diff2/xlnth
         w2    = diff1/xlnth
         xerrorlat = errorlat(j5x,icday)*w1 + errorlat(j5x,icday+1)*w2
         goto 47
        endif
       endif

 35   continue

 47   continue

c     Add zenith angle correction to IR

      iz = (sza-SZAANGMN) + 1
      it = (Tb-SZATMPMN) + 1

      if (iz.ge.1)then
       if (iz.gt.ZABINMX)iz=ZABINMX
       if (it.lt.1)it=1
       if (it.gt.TMPBINMX)it=TMPBINMX
       Tbcorr = Tb - xerrorlat*error(iz,it)
      else
       Tbcorr = Tb
      endif

c     if sza > 87 of if Tb is non zero set to missing

      if (sza.gt.ZAMAX.or.sza.lt.ZERO.or.Tb.lt.ZERO) then
       Tb     = XMISSING
       Tbcorr = XMISSING
      endif

      return

      end