c Merges ISCCP "B1" 0.10 x 0.10 deg lat/lon data into global merged file, c applies zenith angle correction c c The input "B1" data are separate files for each geo satellite separately, c but contained in global arrays (3600 x 1800) for each 3 hr synoptic period c parameter (resx=0.10, resy=0.10, nx=3600, ny=1800, jhmax=ny/2, & i5max=720, j5max=360,CTMPBINMX=34,INMX=6, IPMX=3, & ZABINMX=64,ITMPBINMX=170, XLAT1STPIX=0.018192, & XLON1STPIX=0.018189, ZABINMX=64,ITMPBINMX=170, c & ZABINMX=64,ITMPBINMX=170, XLAT1STPIX=0.0, c & XLON1STPIX=0.0,YSCALE=0.5, XSCALE=0.5, & XMISSING=-9999., ZAPRLMAX=83.0, ISEASMX=8,tbmin=130.) real*4 tb(nx,ny), szapix(nx,ny), sza(i5max,j5max), & szacnt(i5max,j5max), calbin(CTMPBINMX,INMX,IPMX), & errorr(ZABINMX,ITMPBINMX), errorlat(J5MAX,ISEASMX), & difflat(J5MAX,ISEASMX), seasjday(ISEASMX), & work(nx,ny) integer*2 ibuf(nx,ny) c Zenith angle variables real*4 xlat, xlon, sublon, angzen, angndr, xmerge, tbmin, & xmergec, Tbcorr c Parallax variables real*8 satlat,satlon,alat,alon,zht,zang,clat,clon real*4 Tsfcmid, Tsfcequ, Tsfc, Tcld, alatx, rad CHARACTER*98 high_level_path, fileout integer*4 kbuf(CTMPBINMX*INMX*IPMX) data seasjday/57.5,104.5,148.5,206.5,268.5,326.5,395,422.5/ equivalence (kbuf(1), calbin(1,1,1)) fileout='../satellite_files/gef/gef_za.2005073100' open(12,file=fileout,access='direct',recl=nx*ny*4) irec=0 c c c...... Get high-level pathnames c call getarg(1,high_level_path) do i=1,98 if(high_level_path(i:i).eq.' ') go to 110 enddo write(*,*)'WARNING: No blank found in high-level path' 110 len=i-1 c c................................................................................... c CC CC Begin MAIN PROGRAM .... call init(nx,ny,szapix,xmissing,i5max,j5max,sza,szacnt,tb) itype=0 ! running on Linux machine c c....... Set standard atmosphere parameters c Tsfcmid = 288.15 !Surface temperature for 45 degrees latitude and poleward Tsfcequ = 305.0 !Surface temperature for equator rad = 57.2958 c isat = 1 ! GMS c sublon = 140.0 c isat = 1 ! GOES-9 (GMS replacement @ 155E) c sublon = 155.0 c isat = 2 ! G-W c sublon = 225.0 isat = 3 ! G-E sublon = 285.0 c isat = 4 ! Met-7 c sublon = 0.0 c isat = 5 ! Met-5 c sublon = 63.0 iyear=2005 imon=7 iday=31 ihour=0 imin=0 call date2doy(iyear,imon,iday,jday) xjday=jday c c............. Read in the ISCCP B1 data & unscale c call readb1(iyear,imon,iday,ihr,ibuf,nx,ny,tbmin,tb) c c....... ISCCP B1 files ordered 180-180 & south > north; reorder to 0-360 & north > south c call reordertb(nx,ny,work,tb) c itype = 1 is to read correct binary byte type for sza models for SGI or HP itype = 0 c Read in Tb, Seasonal, Latitudinal, sza dependent IR correction models call szamodels(errorr,errorlat,difflat,itype, & high_level_path,len,J5MAX,ZABINMX,ITMPBINMX,ISEASMX) c set standard atmosphere parameters c Surface temperature for 45 degrees latitude and poleward Tsfcmid = 288.15 c Surface temperature for equator Tsfcequ = 305.0 rad = 57.2958 c c........ compute the zenith angles for each location c call zangget(tb,nx,ny,resx,resy,sublon,xmissing,szapix) c c........... Average satellite zenith angle into 0.5 lat/lon bins c write(*,*) 'here2' call zangavg(sza,szacnt,i5max,j5max,szapix,nx,ny,resx,resy, & xscale,yscale,XLON1STPIX,XLAT1STPIX,xmissing) open (77,file='sza_half_deg.data',access='direct', & recl=i5max*j5max*4) write(77,rec=1) sza c Perform satellite zenith angle correction do j = 1,ny xlat = (j-1)*resy - 90. do 250 i = 1,nx xlon=(i-1)*resx if(tb(i,j).ne.XMISSING) then angzen = szapix(i,j) if (izang.eq.IMISSING) then j5 = ((((j-1)*1.0*resy)+XLAT1STPIX)/YSCALE)+1 i5 = ((((i-1)*1.0*resx)+XLON1STPIX)/XSCALE)+1 angzen = sza(i5,j5) endif if (angzen.gt.ZAPRLMAX) tb(i,j)=XMISSING tir=tb(i,j) call szacorr(errorr,errorlat,difflat,seasjday,J5MAX,ISEASMX, & ZABINMX,ITMPBINMX,XMISSING,tir,angzen,xjday,xlat, & Tbcorr,xlon) if (Tbcorr.gt.XMISSING) then tb(i,j)=Tbcorr else tb(i,j)=XMISSING endif endif 250 continue enddo -------------------------------------------------------------------------------------- subroutine szacorr(errorr,errorlat,difflat,seasjday,J5MAX,ISEASMX, & ZABINMX,ITMPBINMX,XMISSING,Tb,sza,xjday,xlat, & Tbcorr,xlon) 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 5. enter "0" if using PC or LINUX, "1" if SGI or HP c c Subroutine outputs c 1. Seasonal, Latitudinal, Tb dependent c Zenith angle IR Tb correction ccc parameter (ZAMAX=85.0, ISTYR=1999, J5MAX=240, XMISSING=-9999, ccc 2 ZABINMX=64,TMPBINMX=170,SZAANGMN=26.,SZATMPMN=160., ccc 3 YEARDYS=365, ISEASMX=8, ZERO=0, ONE=1, XLATMX=60, XLATMN=-60) parameter (ZAMAX=85.0, ISTYR=1999, SZAANGMN=26.,SZATMPMN=160., 3 YEARDYS=365, ZERO=0, ONE=1, XLATMX=90, XLATMN=-90) real errorr(ZABINMX,ITMPBINMX),errorlat(J5MAX,ISEASMX),xjday, & difflat(J5MAX,ISEASMX), seasjday(ISEASMX), LATBININT character*98 regridin, errorlatin character*3 mnth, month(12) real*8 zht,zang,clat,clon ccc 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.ITMPBINMX)it=ITMPBINMX Tbcorr = Tb - xerrorlat*errorr(iz,it) else Tbcorr = Tb endif c if sza > 85 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 zangavg(sza,szacnt,i5max,j5max,szapix,nx,ny,resx,resy, & xscale,yscale,XLON1STPIX,XLAT1STPIX,xmissing) c c Average satellite zenith angle into 0.5 lat/lon bins c real*4 sza(i5max,j5max), szacnt(i5max,j5max), szapix(nx,ny) c do j = 1,ny j5 = ((((j-1)*1.0*resy)+XLAT1STPIX)/YSCALE) + 1 do i = 1,nx i5 = ((((i-1)*1.0*resx)+XLON1STPIX)/XSCALE)+1 if(i5.lt.1 .or. i5.gt.i5max) write(*,*) 'i5=',i5 if(szapix(i,j).ne.xmissing) then sza(i5,j5) = sza(i5,j5) + szapix(i,j) szacnt(i5,j5) = szacnt(i5,j5) + 1. endif enddo enddo do j5 = 1, J5MAX do i5 = 1, I5MAX if (szacnt(i5,j5).gt.0.) then sza(i5,j5) = sza(i5,j5)/szacnt(i5,j5) else sza(i5,j5) = xmissing endif enddo enddo return end subroutine zangget(tb,nx,ny,resx,resy,sublon,xmissing,szapix) real*4 szapix(nx,ny), tb(nx,ny) c c.........Computes zeith angles for each (i,j); "sublon" is the longitude of c at the sub-satellite point of the satellite c c c....... Set standard atmosphere parameters c Tsfcmid = 288.15 !Surface temperature for 45 degrees latitude and poleward Tsfcequ = 305.0 !Surface temperature for equator rad = 57.2958 do j =1,ny do i = 1, nx if(tb(i,j).gt. 0.) then xlat = (j-1)*resy - 90. xlon = (i-1)*resx call SATANG(XLAT, XLON, SUBLON, szapix(i,j), ANGNDR) if(szapix(i,j).lt.0.) write(*,*) i,j,xlat,xlon,szapix(i,j) else szapix(i,j)=xmissing endif c enddo enddo return end subroutine szamodels(errorr,errorlat,difflat,itype, & high_level_path,len,J5MAX,ZABINMX,ITMPBINMX,ISEASMX) c programmer: Bob Joyce c date created: 3/2/2000 c purpose: Subroutine read zenith angle corrections models c Modifications: Janowiak (2/8/2007) - zenith angle model inputs into c "errorlat" array are only 60N-60S ... this routine was c modified to use the values at 60N/60S at all latitudes c poleward of 60N/60S parameter (J5MAX=240,ZABINMX=64,TMPBINMX=170,ISEASMX=8) real difflat(J5MAX,ISEASMX) real errorr(ZABINMX,TMPBINMX), errorlat(J5MAX,ISEASMX) real bufdif(240,8), buferr(240,8) character*98 regridin, errorlatin, high_level_path c Read in IR correction table regridin = high_level_path(1:len) // '/error_regrid_1deg-sza_resiz &e_met+goes_0524-06.plx' open (unit=10, file = regridin, access = 'direct', 1 form='unformatted', recl=64*170*4) read (10,rec=1) errorr c Read in latitude adjustment errorlatin = 1high_level_path(1:len) // '/latitudnal-error_annual-cycle-model_22 26-226.met+goes.plx' open (unit=22, file = errorlatin, & access = 'direct',form='unformatted', recl=240*2*4) do 50 iseas = 1, ISEASMX read (22, rec=iseas) (buferr(j5,iseas),j5=1,240), & (bufdif(j5,iseas),j5=1,240) 50 continue c if (itype.eq.0) then c call SWAP46(errorlat,240*8) c call SWAP46(difflat,240*8) c endif call SWAP46(buferr,240*8) call SWAP46(bufdif,240*8) c c....... Set values poleward of 60N/S to the 60N/S value c do j= 1,60 do k=1,ISEASMX errorlat(j,k)=buferr(1,k) difflat(j,k)=bufdif(1,k) enddo enddo do j=61,300 do k=1,ISEASMX errorlat(j,k)=buferr(j-60,k) difflat(j,k)=bufdif(j-60,k) enddo enddo do j= 301,360 do k=1,ISEASMX errorlat(j,k)=buferr(240,k) difflat(j,k)=bufdif(240,k) enddo enddo return end subroutine readb1(iyear,imon,iday,ihr,ibuf,nx,ny,tbmin,tb) real*4 tb(nx,ny) integer*2 ibuf(nx,ny) c open(1,file='../data/GPCP.B1.a.GOE-9.2005.07.31.0000.GRID' c open(1,file='../data/GPCP.B1.a.GOE-10.2005.07.31.0000.GRID' open(1,file='../data/GPCP.B1.a.GOE-12.2005.07.31.0000.GRID' c open(1,file='../data/GPCP.B1.a.MET-7.2005.07.31.0000.GRID' c open(1,file='../data/GPCP.B1.a.MET-5.2005.07.31.0000.GRID' & ,access='direct',recl=nx*ny*2) read(1,rec=1) ibuf do j=1,ny do i=1,nx tb(i,j)=ibuf(i,j) tb(i,j)=tb(i,j)*0.01 + 100. if(tb(i,j).lt.tbmin) tb(i,j)=-9999. enddo enddo return end subroutine reordertb(nx,ny,work,tb) real*4 tb(nx,ny), work(nx,ny) c c Data enters this subroutine ordered in longitude from 180-180 & exits c ordered from 0-360; also reorders data so that they are c arranged from north to south c do j=1,ny do i=1,nx ii=i+nx/2 if(ii.gt.nx) ii=ii-nx work(ii,j)=tb(i,j) enddo enddo c do j=1,ny jj=ny+1-j do i=1,nx tb(i,j)=work(i,jj) enddo enddo c return end c c-------------------------------------------------------------------------------------- c c SUBROUTINE date2doy(NYY,NMM,NDD,NCEN) C c computes day-of-year from year(NYY), month(NMM), and day(NDD) c INTEGER*2 NTAB(12) C DATA NTAB/0,31,59,90,120,151,181,212,243,273,304,334/ C LP = 0 IF (MOD(NYY,4).EQ.0 ) LP = 1 NLP = NYY / 4 IF (LP.EQ.1 .AND. NLP.GT.0) NLP = NLP - 1 IF (LP.EQ.1 .AND. NMM.LT.3) LP = 0 JDAY = NTAB(NMM) + NDD + LP NCEN = JDAY RETURN END c------------------------------------------------------------- SUBROUTINE SWAP46(A,N) C C REVERSE ORDER OF BYTES IN INTEGER*4 WORD, or REAL*4 C real*4 A(N),ITEMP C CHARACTER*1 JTEMP(4) CHARACTER*1 KTEMP C SAVE C EQUIVALENCE (JTEMP(1),ITEMP) C 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