c--- 05/04/01 - turn off distance calculation; setup dist=1.0 c--- 05/03/01 - calculate and output satellite zenith angle c--- 04/16/01 - output surface flag category (int) real lat_s,lon_s parameter(nsa=13,nod=5) integer n_site real sdepth parameter(n_site=0,sdepth=0.0) integer ibuf(525) common/frame/orbit,time,seq,sync,year,day,alt,nadir, 1 ref(37),sza(37),lat(37),long(37),phi(37), 2 nwl1(37),nwl2(37),nwl3(37),nwl4(37),nwl5(37),nwl6(37), 3 isnow(37),o3bst(37),clfrac(37),pres(37), 4 senswl1(37),senswl2(37),senswl3(37),senswl4(37),senswl5(37), 5 resnwl1(37),resnwl2(37),resnwl3(37),resnwl4(37),resnwl5(37), 6 dndrwl1(37),dndrwl2(37),dndrwl3(37),dndrwl4(37),dndrwl5(37), 7 dndrwl6(37),algflg(37),errflg(37),ozcld(37),soi(37),pcloud(37), 8 prfrac(37),pthir(37),surf(37) integer orbit, time,seq,sync,year,day,alt,algflg,errflg,surf real lat,long,nwl1,nwl2,nwl3,nwl4,nwl5,nwl6 real*8 sat_toms,rtod,dtor real*8 dpi /3.14159265/ real nadir logical*1 ex,jn7,jm3,jep,jad,jall character dsnamei*40,dirn7t*36,local*29 character dsnameo*17, fname*5 character yy*4,ss*2,dd*3,suf*4,rootdir*25,rootn7t*21,yy2*2 dimension qa1smo(nsa,nod,9,4,7,36),qa2smo(nsa,nod,9,4,7,36), 1 qa1des(nsa,nod,9,4,7,36),qa2des(nsa,nod,9,4,7,36) open(unit=15,file='sat_input.dat',status='unknown') ss='/s' read(15,*)jn7 read(15,*)jm3 read(15,*)jep read(15,*)jad pi = 3.141592654 fact = pi/180. nadpos = 18 npos = 35 if (jn7) then rootn7t='/misc/isdat1/n7/v2d/y' suf ='.n7t' endif if (jm3) then rootn7t='/archive/m3toms/v2rep/y' zfact = 1.188 suf = '.m3t' endif if (jep) then c rootdir='/misc/mhdat6/ep/nrt/v2d/y' rootdir='/misc/isdat2/ep/nrt/v2d/y' suf='.ept' endif if (jad) then rootdir='/misc/mhdat7/a1/nrt/v2d/y' suf = '.a1t' nadpos = 19 npos = 37 endif read(15,*)refco read(15,*)aico read(15,*)scaco read(15,*)xlatlo,xlathi read(15,*)xlonlo,xlonhi read(15,*)kref read(15,*)refnom read(15,*)iyr read(15,*)idayf,idayl read(15,'(a5)') fname read(15,*) lat_s read(15,*) lon_s c write(6,*) lat_s,lon_s call convtc(iyr,4,yy) iyr2=iyr-1900 if(iyr.ge.2000)iyr2=iyr-2000 call convtc(iyr2,2,yy2) do 2100 jday = idayf,idayl call convtc(jday,3,dd) dsnamei=rootdir//yy//ss//yy2//dd//suf dirn7t=rootn7t//yy//ss//yy2//dd//suf dsnameo=fname//yy2//dd//suf//'lv2' if (jn7.or.jm3) then inquire(file=dirn7t,exist=ex) if (ex) then open(unit=30,file=dirn7t,status='old',form='unformatted') c open(unit=11,name=dsnameo,status='unknown') c write(6,'(a40,a20)')dirn7t,dsnameo else write(6,*)'File ',dirn7t,' does not exist' go to 2100 endif else inquire(file=dsnamei,exist=ex) if (ex) then open(unit=30,file=dsnamei,status='old',form='unformatted') c open(unit=11,name=dsnameo,status='unknown') c write(6,'(a40,a20)')dsnamei,dsnameo else write(6,*)'File ',dsnamei,' does not exist' go to 2100 endif endif irec=0 100 read(30,end=500)ibuf call rdrec(ibuf) if(seq.le.-1)go to 100 do 200 isp = 1,npos c.. Check for geographical boundaries if((lat(isp).ge.xlatlo.and.lat(isp).le.xlathi).and. 1 (long(isp).ge.xlonlo.and.long(isp).le.xlonhi))then c.. Check quality flags if((errflg(isp).ge.0.and.errflg(isp).le.5) . or . 1 (errflg(isp).ge.10.and.errflg(isp).le.15))then irec=irec+1 if(irec.eq.1)then open(unit=11,name=dsnameo,status='unknown') c23456789012345678901234567890123456789012345678901234567890123456789012 write(11,*) 'yr,day,time, satza, lat, long, ptr, pcl, sza, phi, 1 o3bst, ref, ai, soi, nwl6, algf, errflg, surf, snow, sun_gli 2nt_ang' c write(6,'(a40,a20)')dsnamei,dsnameo endif c-- Compute distance from the center of fov to the site radeg=57.2957795 ! rad to deg conversion dlat= (lat(isp) - lat_s)/radeg dlon= (long(isp) - lon_s)/radeg c write(6,*) lat_s,lon_s c write(6,*) lat(isp),long(isp) c write(6,*) dlat,dlon c dist= 2.*radeg*asin(sqrt( c & sin(dlat/2.)**2 + c & sin(dlon/2.)**2 * cos(lat_s/radeg) * cos(lat_s/radeg))) c idist= NINT(dist*111.2) !Convert to km c NK remove distance calculation to get a global map dist=1.0 ! set up dist=1.0 c dist=6370.*sqrt( c & dlat**2 + dlon**2 * cos(lat_s/radeg) * cos(lat_s/radeg) ) c idist= NINT(dist) !Convert to integer c write(6,*) 'dist= ',dist c-- calculate satellite zenith angle from scan position: sat_toms, deg c-- This is only valid for N7 or EP TOMS since they both have 35 scan positions c-- and nadir poition =18 (ADEOS has 37 positions!! nadir: isp=19) c-- !! This calculation is not valid for EP/TOMS till Dec 97 due to low orbit rtod=180./dpi dtor=dpi/180. sat_toms = rtod*asin(1.15*sin(dtor*3.*abs(isp-nadpos) ) ) c write(6,*) 'isp=',isp,' sat_toms=',sat_toms cc write(11,900)year,day,time,isp,lat(isp),long(isp), c write(11,900)year,day,time,sat_toms,lat(isp),long(isp), c 1 pres(isp),pcloud(isp),sza(isp),phi(isp),o3bst(isp), c 2 ref(isp),resnwl3(isp),resnwl4(isp),resnwl5(isp), c 3 soi(isp),nwl6(isp),dist,algflg(isp),errflg(isp), c 4 surf(isp) c 900 format(2i4,i6,1x,f8.5,1x,f7.2,1x,f7.2,1x,f4.2,1x,f4.2,1x,f5.2, cc 900 format(2i4,i6,i3,1x,f7.2,1x,f7.2,1x,f4.2,1x,f4.2,1x,f5.2,1x, c 1 1x,f5.1,1x,f6.1,1x,f5.1,1x,3(f5.1,1x),f5.1,1x,f6.2,1x,f6.1,1x, c 2 i3,1x,i11,1x,i11) c--- ai calculation IF (year .LT. 1994) THEN ai=resnwl4(isp) ELSE ai=resnwl5(isp) END IF c--- sun glint angle calculation sat_r =sat_toms/radeg sza_r = sza(isp)/radeg phi_r = phi(isp)/radeg myu=cos(sat_r) myu0=cos(sza_r) cglint=myu*myu0 + 1 sqrt(1. - myu*myu)*sqrt(1.-myu0*myu0)*cos(phi_r) sun_glint_ang=radeg*acos(cglint) c--- write test_global.input if (o3bst(isp) .ge. 1000.0 .or. o3bst(isp) .lt. 0.0) go to 77 write(11,901)year,day,time,isp,sat_toms,lat(isp),long(isp), 1 pres(isp),pcloud(isp),sza(isp),phi(isp),o3bst(isp),ref(isp), 2 ai,soi(isp), 3 nwl6(isp), algflg(isp),errflg(isp),surf(isp), 4 sdepth,sun_glint_ang 901 format(2i4,i6,1x,i2,1x,f8.5,1x,f7.2,f8.2,1x, 1 2(f4.2,1x),4(f6.2,1x),2(f5.1,1x), 2 f6.2,1x,i3,1x,i3,1x,i2,1x,f5.1,1x,f5.1) 77 continue endif endif 200 continue go to 100 500 close(30) close(11) 2100 continue stop end subroutine rdrec(ibuf) c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c author c colin seftor c c purpose c read and unpack 1 record of version 7 level 2 data c if record is a data record store in common block /frame/ c otherwise don't do anything c c variables c name type i/o description c ---- ---- --- ------------ c orbit int o orbit number c time int o gmt time (in seconds) c seq int o logical sequence number (> 0 for data c records, < -1 for orbital summary records, c = -1 for trailer record) c sync int o chopper sync flag c year int o year c day int o julian day of year c lat real o latitude c long real o longitude c sza real o solar zenith angle c phi real o azimuthal angle c nwl1 real o wl1 nm n value c nwl2 real o wl2 nm n value c nwl3 real o wl3 nm n value c nwl4 real o wl4 nm n value c nwl5 real o wl5 nm n value c nwl6 real o wl6 nm n value c senswl1 real o wl1 nm sensitivity c senswl2 real o wl2 nm sensitivity c senswl3 real o wl3 nm sensitivity c senswl4 real o wl4 nm sensitivity c senswl5 real o 360 nm sensitivity c ref real o reflectivity c o3bst real o ozone c dndrwl1 real o wl1 nm sensitivity to ref c dndrwl2 real o wl2 nm sensitivity to ref c dndrwl3 real o wl3 nm sensitivity to ref c dndrwl4 real o wl4 nm sensitivity to ref c dndrwl5 real o wl5 nm sensitivity to ref c dndrwl6 real o wl6 nm sensitivity to ref c pthir real o thir cloud pressure (when available) c pres real o terrain pressure (when available) c resnwl1 real o wl1 nm residue c resnwl2 real o wl2 nm residue c resnwl3 real o wl3 nm residue c resnwl4 real o wl4 nm residue c resnwl5 real o wl5 nm residue c ozcld real o ozone below cloud c soi real o so2 index c pcloud real o cloud pressure c algflg int o algorithm flag (0 if a triplet used, 1 if b, c 2 if b with profile mixing, and 3 if c) c errflg int o error flag c prfrac int o profile mixing fraction c surf int o surface category c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c integer ibuf(525) common /frame/ orbit,time,seq,sync,year,day,alt,nadir, 1 ref(37),sza(37),lat(37),long(37),phi(37),nwl1(37),nwl2(37), 2 nwl3(37),nwl4(37),nwl5(37),nwl6(37),isnow(37),o3bst(37), 3 clfrac(37),pres(37),senswl1(37),senswl2(37), 4 senswl3(37),senswl4(37),senswl5(37),resnwl1(37),resnwl2(37), 5 resnwl3(37),resnwl4(37),resnwl5(37),dndrwl1(37),dndrwl2(37), 6 dndrwl3(37),dndrwl4(37),dndrwl5(37),dndrwl6(37),algflg(37), 7 errflg(37),ozcld(37),soi(37),pcloud(37),prfrac(37),pthir(37), 8 surf(37) integer orbit,time,seq,sync,year,day,alt,algflg,errflg,surf real lat, long, nwl1, nwl2, nwl3, nwl4, nwl5, nwl6 real nadir c -- read record into buffer c read (lun) ibuf c -- unpack values that pertain to the whole scan orb = ibuf(1) time = ibuf(2) seq = iunp2b(1,ibuf(3)) sync = iunp2b(2,ibuf(3)) day = iunp2b(1,ibuf(4)) year = iunp2b(2,ibuf(4)) alt = iunp2b(1,ibuf(5)) nadir = iunp2b(2,ibuf(5))/100. c -- if not a data record return if (seq.eq.-1) return c -- unpack values for each sample in the scan do 100 i=1,37 kk = 6 + (i-1)*14 lat(i) = iunp2b(1,ibuf(kk))/100. long(i) = iunp2b(2,ibuf(kk))/100. sza(i) = iunp2b(1,ibuf(kk+1))/100. phi(i) = iunp2b(2,ibuf(kk+1))/100. nwl1(i) = iunp2b(1,ibuf(kk+2))/50. nwl2(i) = iunp2b(2,ibuf(kk+2))/50. nwl3(i) = iunp2b(1,ibuf(kk+3))/50. nwl4(i) = iunp2b(2,ibuf(kk+3))/50. nwl5(i) = iunp2b(1,ibuf(kk+4))/50. nwl6(i) = iunp2b(2,ibuf(kk+4))/50. senswl1(i) = iunp2b(1,ibuf(kk+5))/10000. senswl2(i) = iunp2b(2,ibuf(kk+5))/10000. senswl3(i) = iunp2b(1,ibuf(kk+6))/10000. senswl4(i) = iunp2b(2,ibuf(kk+6))/10000. senswl5(i) = iunp2b(1,ibuf(kk+7))/10000. ref(i) = iunp2b(2,ibuf(kk+7))/100. o3bst(i) = iunp2b(1,ibuf(kk+8))/10. errflg(i) = iunp2b(2,ibuf(kk+8)) dndrwl1(i) = iunp1b(1,ibuf(kk+9))/-50. dndrwl2(i) = iunp1b(2,ibuf(kk+9))/-50. dndrwl3(i) = iunp1b(3,ibuf(kk+9))/-50. dndrwl4(i) = iunp1b(4,ibuf(kk+9))/-50. dndrwl5(i) = iunp1b(1,ibuf(kk+10))/-50. dndrwl6(i) = iunp1b(2,ibuf(kk+10))/-50. pthir(i) = iunp1b(3,ibuf(kk+10))/100. pres(i) = iunp1b(4,ibuf(kk+10))/100. resnwl1(i) = (iunp1b(1,ibuf(kk+11))-127.)/10. resnwl2(i) = (iunp1b(2,ibuf(kk+11))-127.)/10. resnwl3(i) = (iunp1b(3,ibuf(kk+11))-127.)/10. resnwl4(i) = (iunp1b(4,ibuf(kk+11))-127.)/10. resnwl5(i) = (iunp1b(1,ibuf(kk+12))-127.)/10. ozcld(i) = iunp1b(2,ibuf(kk+12)) soi(i) = iunp1b(3,ibuf(kk+12))-50. pcloud(i) = iunp1b(4,ibuf(kk+12))/100. algflg(i) = iunp1b(1,ibuf(kk+13)) clfrac(i) = iunp1b(2,ibuf(kk+13))/100. prfrac(i) = iunp1b(3,ibuf(kk+13))/10. surf(i) = iunp1b(4,ibuf(kk+13)) 100 continue return end function iunp1b (lbyte,ipackd) c c*********************************************************************** c* c* purpose: c* return the value of the specified byte as a fullword integer c* c* arguments: c* name i/o description c* ---- --- ----------- c* lbyte i byte location (1 to 4) of packed value c* ipackd i fullword containing selected byte c* c* date: jan 1991 programmer: b. kelley c* c*********************************************************************** c c** right-shift input value to move selected byte to rightmost byte c** field, then use logical-and to get value of that byte c nbits = (4-lbyte) * 8 iunp1b = iand(ishft(ipackd,-nbits),255) c end function iunp2b (ihalf,ipackd) c c*********************************************************************** c* c* purpose: return the value of specified two-byte half of input c* integer as a fullword integer c* c* arguments: c* name i/o description c* ---- --- ----------- c* ihalf i specifies 2-byte half of word to be extracted c* (1 for first half, 2 for second half) c* ipackd i packed word c* c* local variables: c* rmask - mask to clear left half of word c* lmask - mask to clear right half of word c* c* date: jan 1991 programmer: b. kelley c* c*********************************************************************** c integer rmask data rmask/65535/,lmask/-65536/ c if (ihalf.eq.1) then c c** iunpak = ipackd right-shifted 16 bits. if ipackd is negative, c** logical-or iunpak with lmask to set leftmost 16 bits. c iunpak=ishft(ipackd,-16) if (ipackd.lt.0) iunpak=ior(lmask,iunpak) else c c** clear left half of ipackd. if iunpak is greater than maximum c** value of two-byte signed integer, then packed value is a c** negative number and is obtained by subtracting 2**16. c iunpak=iand(ipackd,rmask) if (iunpak.gt.32767) iunpak=iunpak-65536 end if c iunp2b = iunpak c end SUBROUTINE CONVTC (NUM,NCHAR, LOC) C*********************************************************************** C LIBRARY SUBROUTINE C C MODULE NAME: CONVTC C C VERSION : 1.0 C C PROGRAMMER : LUCY LIU AND MICHAEL PENG, STX, AUGUST 1990 C C PURPOSE: TO CONVERT FIRST NCHAR DIGITS OF AN INTEGER INTO A CHARACTER C OR A STRING. (NOTE: NCHAR DOESN'T INCLUDE THE SIGN IF NUM C IS NEGATIVE) C C CALLING SEQUENCE: CONVTC (NUM,NCHAR, LOC) C C SUBROUTINES CALLED: NONE C C INTRINSIC FUNCTIONS USED: CHAR, ICHAR, IABS C C COMMON BLOCKS USED: NONE C C ARGUMENTS AND LOCAL VARIABLES: C C NAME TYPE I/O DESCRIPTIONS C --------- ---- --- ---------------------------------------------- C NUM I*4 I INTEGER NUMBER FOR CONVERSION C NCHAR I*4 I NUMBER OF CHARS TO BE CONVERTED C LOC C*X O CHAR VARIABLE OR ARRAY (STRING) C ISTART I*4 START POSITION OF INTEGRAL CHAR STRING: C 1 FOR POSITIVE INTEGER; 2 FOR NEGATIVE INTEGER C IREM I*4 ABSOLUTE VALUE OR REMAINDER BY 10 C ITEMP I*4 TEMPORARY BUFFER C*********************************************************************** C C --- ARGUMENTS C INTEGER NUM, NCHAR CHARACTER LOC(1) C C --- LOCAL VARIABLES C INTEGER ISTART, IREM, ITEMP C IREM = IABS (NUM ) ISTART = 1 C C --- KEEP THE SIGN IF THE INTEGER IS NEGATIVE C IF (NUM.LT.0) THEN NCHAR = NCHAR + 1 LOC(1) = '-' ISTART = ISTART + 1 ENDIF C C --- CONVERT DIGIT BY DIGIT STARTING WITH THE LEAST SIGNIFICANT DIGIT C DO 100 I = NCHAR, ISTART, -1 ITEMP = IREM / 10 LOC(I) = CHAR (IREM - ITEMP * 10 + ICHAR('0')) IREM = ITEMP 100 CONTINUE C C --- PUT AN ASTERISK AT THE FIRST POSITION IF THE ACTUAL NUMBER OF C --- DIGITS IN THE INTEGER IS FOUND GREATER THAN NCHAR C IF (IREM.NE.0) LOC(1) = '*' C RETURN END