!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! 2|---------------------------------------------> satellite: (see below) ! F:\DVD_SBUV| ----------------------------------> drive letter for DVD (or data source) ! -----| 19.53 |---| 2.50|-------------------> center latitude & Dlat (-90 to +90) ! -----| -155.58|---| 20.00|-------------------> center longitude & Dlon (-180 to +180) ! T|---| 02/02/1985|---| 03/01/1985|-------------> (option 1) select time period to analyze: begin time & end time (MM DD YYYY) ! F|---| 0.00|---| + 60.00|-------------------> (option 2) screen on solar zenith angle: min ZA & max ZA (0 to 90) ! F|---------------------------------------------> (option 3) use only flag 0 and 1 data (and 10 and 11) ! F|---| 0.20|--------------------------------> (option 4) screen on Quality factor ! T|---------------------------------------------> (option 5) do distance weighted average in addition to close scans only ! MaunaLoa_tst.ov1 |------> output file name for individual overpass scans near site each day ! MaunaLoa_tst.ovp |------> output file name for distance weighted overpass average ! =========================================================================================================== ! for logical options 1 - 5 set T(rue) or F(alse) | ! | ! satellite: | ! 1: Nimbus 7 10/31/1978 - 6/21/1990 | ! 2: NOAA 9a 2/02/1985 - 12/31/1989 | ! 3: NOAA 11 12/01/1988 - 3/27/2001 | ! 4: NOAA 9d 1/01/1993 - 2/19/1998 | ! 5: NOAA 16 10/03/2000 - 12/31/2003 | ! ================================================== ! ! Program to read SBUV version 8 ASCII data set and find overpasses ! (38 values saved in outarr for each of up to 24 FOVs in the lat/lon window each day) ! INTEGER y1, y2, II(4) REAL*4 pmr(15), MR(15), plyr(13), Ozlyr(13), outarr(38,24), wt(24) LOGICAL CTL(5) CHARACTER*8 DATE,dummy CHARACTER*12 drive CHARACTER*17 NAMlist CHARACTER*40 OUTscan, OUTovp_av CHARACTER*41 nam2 CHARACTER*43 FNAM CHARACTER*53 FNAM2 Data pmr/0.5,0.7,1.0,1.5,2.0,3.0,4.0,5.0,7.0,10.0,15.0,20.0,30.0,40.0,50.0/ !hPa !======================================= read input from CONTROL.DAT ============================================== open(1,FILE='CONTROL_OVP.DAT') READ(1,'(1I1)') IS ! select satellite READ(1,'(1A12)') drive ! set drive letter for DVD (or data source) READ(1,'(6x,F8.2,5x,F8.2))') alat0,Dlat ! set center latitude and range alatmin=alat0-Dlat ; alatmax=alat0+Dlat if(alatmin.lt.-90.0) alatmin=-90.0 if(alatmax.gt.90.0) alatmax=90.0 READ(1,'(6x,F8.2,5x,F8.2))') alon0,Dlon ! set center longitude and range alonmin=alon0-Dlon ; alonmax=alon0+Dlon READ(1,'(1L1,2(5x,I3,1x,I2,1x,I4))') CTL(1),imb,idb,iyb,ime,ide,iye ! select time period to read READ(1,'(1L1,5x,F8.2,5x,F8.2))') CTL(2),ZAmin,ZAmax ! screen by solar zenith angle READ(1,'(1L1)') CTL(3) ! use error flag 0 and 1 (and 10 and 11) only READ(1,'(1L1,5x,F8.2)') CTL(4),QClimit ! screen on residue QC limit READ(1,'(1L1)') CTL(5) ! do distance weighted average rather than closest READ(1,'(1A40)') OUTscan ! name of output data set for scans near ground site READ(1,'(1A40)') OUTovp_av ! name of output data set - overpass average close(1) !================================================================================================================== if(IS.eq.1) then ! Nimbus 7 y1=1978 ; y2=1990 NAMlist='N07_filenames.txt' write(*,*) ' Nimbus 7: data coverage from 1978 10 31 to 1990 06 21' elseif(IS.eq.2) then ! NOAA 9a y1=1985 ; y2=1989 NAMlist='N9a_filenames.txt' write(*,*) ' NOAA 9a: data coverage from 1985 02 02 to 1989 12 31' elseif(IS.eq.3) then ! NOAA 11 y1=1988 ; y2=2001 NAMlist='N11_filenames.txt' write(*,*) ' NOAA 11: data coverage from 1988 12 01 to 2001 03 27' elseif(IS.eq.4) then ! NOAA 9d y1=1993 ; y2=1998 NAMlist='N9d_filenames.txt' write(*,*) ' NOAA 9d: data coverage from 1993 01 01 to 1998 02 19' elseif(IS.eq.5) then ! NOAA 16 y1=2000 ; y2=2003 NAMlist='N16_filenames.txt' write(*,*) ' NOAA 16: data coverage from 2000 10 03 to 2003 12 31' else write(*,*) ' Did not enter a valid satellite number' stop endif write(*,'('' control parameters: '',5L1)') CTL if(CTL(1)) write(*,*) ' CTL(1): select time period to process' if(CTL(2)) write(*,*) ' CTL(2): select solar zenith angle range to process' if(CTL(3)) write(*,*) ' CTL(3): select good total ozone scans only (0,1,10,11)' if(CTL(4)) write(*,*) ' CTL(4): screen out scans with QC factor over limit' if(CTL(5)) write(*,*) ' CTL(5): output distance weighted average overpass ozone also' write(*,*) ' ' ! ======================= do a simple check on input parameters ===================================== if(CTL(1)) then if(iyb.lt.y1 .or. iyb.gt.y2) then write(*,'('' >>> Bad year number'',i5)') iyb stop endif if(imb.lt.1 .or. imb.gt.12) then write(*,'(''>>> Bad month'',I4)') imb stop endif if(idb.lt.1 .or. idb.gt.31) then write(*,'(''>>> Bad day'',I4)') idb stop endif if(iye.lt.y1 .or. iye.gt.y2) then write(*,'('' >>> Bad year number'',i5)') iyb stop endif if(ime.lt.1 .or. ime.gt.12) then write(*,'(''>>> Bad month'',I4)') imb stop endif if(ide.lt.1 .or. ide.gt.31) then write(*,'(''>>> Bad day'',I4)') idb stop endif endif if(alatmin.lt.-90.1 .or. alatmax.gt.90.1) then write(*,'(''>>> Bad latitude limit'',2F8.1)') alatmin,alatmax go to 999 endif if(alonmin.lt.-180.1 .or. alonmax.gt.180.1) then write(*,'(''>>> Bad longitude limit'',2F8.1)') alonmin,alonmax go to 999 endif if(CTL(2)) then if(ZAmin.lt.0.0 .or. ZAmax.gt.90.1) then write(*,'(''>>> Bad solar zenith angle limit'',2F8.1)') ZAmin,ZAmax go to 999 endif endif ! ============================================================================ get ready to read data ! Read series of scans Ngood=0; Nscan=0; NN=0; IDS=0 if(CTL(5)) then write(*,'('' opening output files: '',1A40,/,23x,1A40)') OUTscan,OUTovp_av open(22,file=OUTovp_av) WRITE(22,'('' date year DOY Sec pts Lat Lon ZA TOZ Ref QC dist'',13I7,4x,15F7.1)') (j,j=1,13),(pmr(j),j=1,15) else write(*,'('' opening output file: '',1A40)') OUTscan endif open(21,file=OUTscan) pause ! ---------------------------------------------------------------------------- loop through daily files open(11,file=NAMlist) ! NAMlist contains path and name of every file for selected satellite dataset 50 continue read(11,'(1A43)',end=500) FNAM if(CTL(1)) then call extractYMD(FNAM,iy,im,id) ! determine expected year, month, and day from file name if(iy.lt.iyb) go to 50 if(iy.eq.iyb .and. im.lt.imb) go to 50 if((iy.eq.iyb .and. im.eq.imb) .and. id.lt.idb) go to 50 if(iy.gt.iye) go to 500 if(iy.eq.iye .and. im.gt.ime) go to 500 if((iy.eq.iye .and. im.eq.ime) .and. id.gt.ide) go to 500 endif ! write(*,'('' opening input file '',1A43)') FNAM ! write(21,'('' input file '',1A43)') FNAM ! ---------------------------------------------------------------------------- loop through data records each day nam2=FNAM(3:43) FNAM2 = drive // nam2 open (2,file=FNAM2) ! there are 8 header lines in each file read(2,'(36x,I3)') iday0 ! first line in each file contains day number read(2,'(1I7)') nrecords read(2,'(1A8)') dummy ! these informational lines are skipped read(2,'(1A8)') dummy read(2,'(1x,8F7.4,2F8.4,F7.4,2F8.4)') plyr read(2,'(1A8)') dummy read(2,'(2x,9F4.1,6F5.1)') pmr read(2,'(1A8)') dummy ! read next scan 100 CONTINUE read(2,2000,end=400) iyear,iday,isec,alat,alon,ZA,TOZ,ref,AI,QC,iflag 2000 format(1x,I4,I4,I6,F7.2,F8.2,F6.2,F6.1,F7.3,F6.1,F6.3,I4) read(2,'(3F7.2,7F7.3,3F7.4)') Ozlyr read(2,'(15F7.3)') MR ! write(*,'(I6,I4,I6,F6.1F7.1,F6.1)') iyear,iday,isec,alat,alon,TOZ Nscan=Nscan+1 ! ********* Do screening here ***************************** if(alat.lt. alatmin .or. alat.gt.alatmax) GO TO 100 ! latitude screen if(alon.lt. alonmin .or. alon.gt.alonmax) GO TO 100 ! longitude screen if(TOZ.le.10.0 .or. TOZ.gt.900.0) go to 100 ! reject scans without valid total ozone ! optional: solar zenith angle screen if(CTL(2) .and. (ZA.LT.zamin .or. ZA.gt.zamax)) GO TO 100 ! optional: reject scans with bad error flags if(CTL(3) .and. (iflag.ne.0 .and. iflag.ne.1 .and. iflag.ne.10 .and. iflag.ne.11)) go to 100 ! optional: reject scans where the residue QC value is too big if(CTL(4) .and. (QC.gt.QClimit)) go to 100 Ngood=Ngood+1 ! ********* passed screens ***************************** IF(IDAY.ne.IDS) THEN CALL SETDAT(IYEAR,IDAY,DATE,ID) write(*,'('' reading data for: '',a8,i3,'','',i5)') DATE,ID,IYEAR write(21,'('' Data for: '',a8,i3,'','',i5)') DATE,ID,IYEAR WRITE(21,*) 'Yr Day Sec Lat Lon ZA TOZ Ref QC dist' IDS=IDAY ENDIF NN=NN+1 outarr(1,NN)=iyear outarr(2,NN)=iday outarr(3,NN)=isec outarr(4,NN)=alat outarr(5,NN)=alon outarr(6,NN)=ZA outarr(7,NN)=TOZ outarr(8,NN)=ref outarr(9,NN)=QC call Distcal(alat0,alon0,alat,alon,dist) outarr(10,NN)=dist DO j=1,13 outarr(j+10,NN)=OZlyr(j) END DO DO j=1,15 outarr(j+23,NN)=MR(j) END DO ! write(*,'(I6,I8,3f8.1)') iday,isec,alat,alon,toz write(21,'(2F5.0,F7.0,3F8.2,F8.1,F8.2,F8.3,F7.0,13F8.3,3x,15F7.3)') (outarr(j,NN),j=1,38) GO TO 100 ! ----------------------- finished reading data for one day 400 CONTINUE close(2) if(.not.CTL(5)) GO TO 450 ! do distance weighted averages if(NN.eq.0) GO TO 450 WTT=0 if(NN.gt.22) then write(*,'('' NN='',I4)') NN NN=22 endif ! get weights for average do j=1,NN dist=outarr(10,j) if(dist.lt.0.1) dist=0.1 ! avoid singularity when distance = 0 wt(j)=1.0/dist WTT=WTT+wt(j) end do wt=wt/WTT do j=1,NN ! calculate weighted average do k=1,38 outarr(k,NN+1)=outarr(k,NN+1)+outarr(k,j)*wt(j) ! put weighted average into the N+1 position in outarr end do end do do k=1,3 ! use integer array for outputting: year day seconds II(k)=outarr(k,NN+1) + 0.5 end do II(4)=outarr(10,NN+1) + 0.5 iyear=outarr(1,NN+1) iday=outarr(2,NN+1) CALL SETDAT(iyear,iday,DATE,ID) write(22,2100) DATE,ID,II(1),II(2),II(3),NN,(outarr(j,NN+1),j=4,9),II(4),(outarr(j,NN+1),j=11,38) 2100 format(1x,A8,I3,2x,2I4,I6,I4,F7.2,F8.2,F7.2,F6.1,2F6.2,I5,2x,6F7.2,7F7.3,3x,15F7.3) 450 continue ! zero the counters and output array NN=0 outarr=0.0 wr=0.0 GO TO 50 500 CONTINUE close(11) close(21) close(22) 900 CONTINUE WRITE(*,'('' Selected '',i6,'' scans out of '',I9,'' scans read'')') Ngood, Nscan WRITE(21,'('' Selected '',i6,'' scans out of '',I9,'' scans read'')') Ngood, Nscan WRITE(*,*) ' Job completion normal' STOP 999 CONTINUE WRITE(*,*) ' Job stopped for errors' WRITE(21,*) ' Job stopped for errors' close(21) close(22) STOP END ! ---------- subroutine to get month and day from day of year ------------------------------------------------- SUBROUTINE SETDAT(IYEAR,IDAY,DATE,ID) CHARACTER*8 DATE,MON(12) DIMENSION IMON(13),IMONL(13) DATA IMON/1,32,60,91,121,152,182,213,244,274,305,335,366/ DATA IMONL/1,32,61,92,122,153,183,214,245,275,306,336,367/ DATA MON/' January','February',' March',' April',' May',' June',' July',' August','Septembr',' October','November','December'/ ! IF(MOD(IYEAR,4).EQ.0) GO TO 200 DO 100 J=1,12 IF(IDAY.GE.IMON(J+1)) GO TO 100 ID=IDAY-IMON(J)+1 GO TO 400 100 CONTINUE ! 200 CONTINUE DO 300 J=1,12 IF(IDAY.GE.IMONL(J+1)) GO TO 300 ID=IDAY-IMONL(J)+1 GO TO 400 300 CONTINUE ! 400 CONTINUE k=1 DATE=MON(J) RETURN END ! ---------- subroutine to extract year, month, and day from file name ---------------------------------------- SUBROUTINE extractYMD(name,iy,im,id) CHARACTER*1 c1,c2,c3,c4,c5,c6,c7,c8 character*43 name ! n:\Data\Nimbus7\Y1978\sbuv_n07_19781031.txt ! 1234567890123456789012345678901234567890123 ! 1 2 3 4 ! c1 = name(32:32) ; c2=name(33:33) ; c3=name(34:34) ; c4=name(35:35) ! year iy = 1000*(ichar(c1)-48) + 100*(ichar(c2)-48) + 10*(ichar(c3)-48) + ichar(c4)-48 c5 = name(36:36) ; c6=name(37:37) ! month im = 10*(ichar(c5)-48) + ichar(c6)-48 c7=name(38:38) ; c8=name(39:39) ! day id = 10*(ichar(c7)-48) + ichar(c8)-48 return end ! ---------- subroutine to calculate distance from ground site in kilometers ----------------------------------- SUBROUTINE Distcal(alat0,alon0,alat,alon,dist) dx=abs(alat-alat0)*111.0 dy=abs(alon-alon0)*111.0*cos(alat0*0.0174533) dist=sqrt(dy*dy + dx*dx) return end