c-------------------------------------------------------------------------- C Merges ISCCP "B1" 0.10 x 0.10 deg lat/lon data into global merged file, c applies zenith angle correction & intercalibrates the TBs. 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.0, & 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.......... Add parallax correction c write(*,*) 'here1' 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 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 c........... Perform satellite zenith angle correction c write(*,*) 'here3' call zangcorr(tb,szapix,nx,ny,sza,i5max,j5max,xmissing, & resx,resy,ZAPRLMAX,errorr,errorlat,seasjday, & ZABINMX,ITMPBINMX,ISEASMX,xjday,Tbcorr, & xscale,yscale,difflat) c write(6,'('' after zang corr'')') do j = 900,920 write(6,'(i8,10(1x,f8.1))') j,(tb(i,j),i=301,310) enddo irec=irec+1 write(12,rec=irec) tb stop 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 c--------------------------------------------------------------------------- c ccc subroutine outfile(high_level_path,len,iyear,mon,nday,hour,imin, c x fileout) c CHARACTER*98 high_level_path c CHARACTER*90 fileout c fileout = high_level_path(1:len) // '/global_ir/gef/ge_fd1_YYYY0M0 c &D0H00.global-sza' c szaout = high_level_path(1:len) // '/global_ir/gef/ge_fd1_YYYY0M0D c &0H00.5sza' c write(fileout(len+23:len+26),'(i4)') IYEAR c write(szaout(len+23:len+26),'(i4)') IYEAR c if (MON.gt.9) then c write(fileout(len+27:len+28),'(i2)') MON c write(szaout(len+27:len+28),'(i2)') MON c else c write(fileout(len+28:len+28),'(i1)') MON c write(szaout(len+28:len+28),'(i1)') MON c endif c if (NDAY.gt.9) then c write(fileout(len+29:len+30),'(i2)') NDAY c write(szaout(len+29:len+30),'(i2)') NDAY c else c write(fileout(len+30:len+30),'(i1)') NDAY c write(szaout(len+30:len+30),'(i1)') NDAY c endif c if (Hour.gt.9) then c write(fileout(len+31:len+32),'(i2)') Hour c write(szaout(len+31:len+32),'(i2)') Hour c else c write(fileout(len+32:len+32),'(i1)') Hour c write(szaout(len+32:len+32),'(i1)') Hour c endif c if (imin.gt.9) then c write(fileout(len+33:len+34),'(i2)') imin c write(szaout(len+33:len+34),'(i2)') imin c else c write(fileout(len+34:len+34),'(i1)') imin c write(szaout(len+34:len+34),'(i1)') imin c endif c c write(6,*) 'output data set =',fileout c c return c end c c-------------------------------------------------------------------------------------- c subroutine parllx (satlat,satlon,alat,alon, + zht,zang,clat,clon) c rewritten in fortran by clay davenport 8-7-97 c from c code parallax.c by russ? irvin? c the program converts lat/lon c-os into Cartesian c-os; c solves for the line joining the satellite and apparent surface point; c finds point on that line a distance zht over the surface; c finds the surface point connecting that with center of earth c converts back to lat/lon implicit real*8 (a-f) implicit real*8 (r-z) dpi=3.1415927 aeqr=6378.077 aplr=6356.577 aratio=aeqr/aplr aorbit=42165.39 r=(aplr+aeqr)/2.0 c covert angles from degrees to radians and height from ft to m satlat=satlat*dpi/180.0 satlon=-satlon*dpi/180.0 alat=alat*dpi/180.0 alon=-alon*dpi/180.0 c zht=zht*.000304806 c get x-y-z coordinates for satellite c slat1 is the geometric latitude for geodetic lat satlat slat1=atan(tan(satlat)*aratio**2) xs=aorbit*cos(slat1)*sin(satlon) ys=aorbit*sin(slat1) zs=aorbit*cos(slat1)*cos(satlon) c get x-y-z coordinates for surface point alat1=atan(tan(alat)*aratio**2) ri=aeqr/sqrt(cos(alat1)**2+aratio**2*sin(alat1)**2) xi=ri*cos(alat1)*sin(alon) yi=ri*sin(alat1) zi=ri*cos(alat1)*cos(alon) c b is the new aratio b=((aeqr+zht)/(aplr+zht))**2 xsmxi=xs-xi ysmyi=ys-yi zsmzi=zs-zi c get zenith angle xfact=sqrt(xsmxi**2+ysmyi**2+zsmzi**2) zang=((xsmxi*xi + ysmyi*yi + zsmzi*zi)/(r*xfact)) zang=acos(zang) zang=zang*180.0/dpi c if (zang.gt.70.0) goto 10 c quadratic equation to solve for the line values e=xsmxi**2 + b*ysmyi**2 + zsmzi**2 ef=2.0*(xi*xsmxi+b*yi*ysmyi+zi*zsmzi) eg=xi**2 + zi**2 + b*yi**2 - (aeqr+zht)**2 a=(sqrt(ef**2-4.0*e*eg)-ef)/2.0/e c corrected xyz xc=xi + a*xsmxi yc=yi + a*ysmyi zc=zi + a*zsmzi c convert back to lat/lon clat1=atan(yc/sqrt(xc**2+zc**2)) clat=atan(tan(clat1)/aratio**2)*180/dpi clon=-atan2(xc,zc)*180/dpi 10 return end c c-------------------------------------------------------------------------------------- c 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 c c-------------------------------------------------------------------------------------- c 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 c c-------------------------------------------------------------------------------------- c 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,ITMPBINMX), 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=ZABINMX*ITMPBINMX*4) read (10,rec=1) errorr ccc if (itype.eq.1) then ccc call SWAP46(errorr,10880) ccc endif 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 ccc if (itype.eq.0) then call SWAP46(buferr,240*8) call SWAP46(bufdif,240*8) ccc endif 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 c c-------------------------------------------------------------------------------------- c subroutine swap17(sza,N) real*4 sza(259200), itemp character*1 jtemp(4) character*1 ktemp save equivalence (jtemp(1), itemp) do 10 i=1, 259200 itemp = sza(i) ktemp = jtemp(4) jtemp(4) = jtemp(1) jtemp(1) = ktemp ktemp = jtemp(3) jtemp(3) = jtemp(2) jtemp(2) = ktemp sza(i) = itemp 10 continue return end c c-------------------------------------------------------------------------------------- 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 init(nx,ny,szapix,xmissing,i5max,j5max,sza,szacnt,tb) real*4 szapix(nx,ny), tb(nx,ny), sza(i5max,j5max), & szacnt(i5max,j5max) do j = 1, ny do i = 1, nx szapix(i,j) = xmissing tb(i,j)=xmissing enddo enddo c c.........Initialize sza 0.5 deg lat/lon arrays c do j5 = 1, J5MAX do i5 = 1, I5MAX sza(i5,j5) = 0. szacnt(i5,j5) = 0. 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 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 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))/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 zangcorr(tb,szapix,nx,ny,sza,i5max,j5max,xmissing, & resx,resy,ZAPRLMAX,errorr,errorlat,seasjday, & ZABINMX,ITMPBINMX,ISEASMX,xjday,Tbcorr, & xscale,yscale,difflat) c c Perform satellite zenith angle correction c real*4 tb(nx,ny), szapix(nx,ny), errorr(ZABINMX,ITMPBINMX), & errorlat(J5MAX,ISEASMX), sza(i5max,j5max), & difflat(J5MAX,ISEASMX), seasjday(ISEASMX) write(*,*) '3.1' 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 (angzen.eq.XMISSING) then j5 = ((((j-1)*1.0*resy))/YSCALE) + 1 if(j5.lt.1 .or. j5.gt.j5max) + write(*,*) j,j5,resy,XLAT1STPIX,YSCALE i5 = ((((i-1)*1.0*resx)+XLON1STPIX)/XSCALE)+1 if(i5.lt.1 .or. i5.gt.i5max) + write(*,*) i,i5,resx,XLON1STPIX,XSCALE angzen = sza(i5,j5) endif c Set any location with sza greater than 83.0 to missing c to prevent warm bias induced from parallax correction 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 return end