program szacr c For itype enter "0" if using PC or LINUX, "1" if SGI or HP parameter (J5MAX=240,ZABINMX=64,TMPBINMX=170,ISEASMX=8) real error(ZABINMX,TMPBINMX), errorlat(J5MAX,ISEASMX) c The example given is for a limb correction for the GOES 8 c satellite (subsat 285E) and a target pixel location of c 41.2 N and 228.0 E with an original brightness temp of 235 Tb = 235 XLAT = 41.2 XLON = 228.0 SUBLON = 285.0 itype = 0 call szamodels (error,errorlat,itype) do 100 jday = 1, 366 call satang (XLAT,XLON,SUBLON,ANGZEN,ANGNDR) xjday = jday call szacorr (error,errorlat,Tb,ANGZEN,xjday,xlat,Tbcorr) write (*,*) ' day ',jday,' sza ', ANGZEN, ' Tb ',Tb, & ' Tbcorr ',Tbcorr 100 continue stop end subroutine szamodels (error,errorlat,itype) c programmer: Bob Joyce c date created: 3/2/2000 c purpose: Subroutine read zenith angle corrections models c c For itype enter "0" if using PC or LINUX, "1" if SGI or HP parameter (J5MAX=240,ZABINMX=64,TMPBINMX=170,ISEASMX=8) real difflat(J5MAX,ISEASMX) real error(ZABINMX,TMPBINMX), errorlat(J5MAX,ISEASMX) character*98 regridin, errorlatin c Read in IR correction table c regridin ="/export/sgi109/irstat/mcidas/data/error_regrid" regridin ="/export-1/lnx57/wd52rj/mcidas/src/sza/error_regrid" open (unit=10, file = regridin, access = 'direct', 1 form='unformatted', recl=64*170*4) read (10,rec=1) error if (itype.eq.1) then call SWAP46(error,10880) endif c Read in latitude adjustment errorlatin = &"/export-1/lnx57/wd52rj/mcidas/src/sza/latitudinal-error_annual-cy &cle-model" open (unit=22, file = errorlatin, & access = 'direct',form='unformatted', recl=240*2*4) do 50 iseas = 1, ISEASMX read (22, rec=iseas) (errorlat(j5,iseas),j5=1,J5MAX), & (difflat(j5,iseas),j5=1,J5MAX) 50 continue if (itype.eq.0) then call SWAP46(errorlat,240*8) endif return end 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 SUBROUTINE SWAP32(A,N) C REVERSE ORDER OF BYTES IN INTEGER*2 WORD INTEGER*2 A(N), ITEMP CHARACTER*1 JTEMP(2) CHARACTER*1 KTEMP SAVE EQUIVALENCE (JTEMP(1),ITEMP) DO 10 I = 1,N ITEMP = A(I) KTEMP = JTEMP(2) JTEMP(2) = JTEMP(1) JTEMP(1) = KTEMP A(I) = ITEMP 10 CONTINUE RETURN END SUBROUTINE SWAP46(A,N) C REVERSE ORDER OF BYTES IN INTEGER*4 WORD, or REAL*4 real*4 A(N),ITEMP CHARACTER*1 JTEMP(4) CHARACTER*1 KTEMP SAVE EQUIVALENCE (JTEMP(1),ITEMP) DO 10 I = 1,N ITEMP = A(I) KTEMP = JTEMP(4) JTEMP(4) = JTEMP(1) JTEMP(1) = KTEMP KTEMP = JTEMP(3) JTEMP(3) = JTEMP(2) JTEMP(2) = KTEMP A(I) = ITEMP 10 CONTINUE RETURN END SUBROUTINE satang (XLAT,XLON,SUBLON,ANGZEN,ANGNDR) DATA TWOPI/6.28318/, R/6371./, H/35680./ DEGRAD=360./TWOPI ZLON=XLON IF(ZLON.LT.0.) ZLON=360.+ZLON ZSUBLN=SUBLON IF(ZSUBLN.LT.0.) ZSUBLN=360.+ZSUBLN DIFLON=ABS(ZLON-ZSUBLN)/DEGRAD DIFLAT=ABS(XLAT)/DEGRAD COSPSI=COS(DIFLON)*COS(DIFLAT) PSI=ACOS(COSPSI) SSQR=R*R + (R+H)*(R+H) - 2*R*(R+H)*COSPSI S=SQRT(SSQR) SINNDR=(R/S)*SIN(PSI) SINZEN=((R+H)/R)*SINNDR ANGZEN=ASIN(SINZEN)*DEGRAD ANGNDR=ASIN(SINNDR)*DEGRAD RETURN END