C Last change: JA 3/26/2015 7:49:05 AM CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C QC_Rad, Version 2.0 C Copyright by Chuck Long C This program applies tests to surface radiation measurement data to assess C data quality. These tests are described in "QC_Rad_qc_tests.doc" and the C QC flags are described in "QC_Rad_flg_key.doc". The program has 5 levels C of testing: physically possible limits, two separate levels of C user configurable climatological limits, cross comparison limits, and user C configurable cross comparison limits. The program reads in configuration C settings from the file "QC_Rad2.cnf", and the list of files to be C processed from the file "QC_Rad2.dir". Output is a file called C "QC_Rad2.asc" which contains daily summaries of the number of data C tests that were possible, and that were failed. (See "QC_Rad_flg_key.doc" C for details.) Other output can be set for individual files corresponding C to the input files (with a "q" added to the extension), and/or one file C with all data called "QC_Rad2.all". The "daily" and "all" output can be set C to include the input measurement values plus only the data QC flags, C or the data QC flags plus all corresponding limits. C C Limit Arrays: C "plim" is physically possible limits, "elim" is 2nd level configurable limits, C "nlim" is the global cross comparison limits, and "clim" is 1st level configurable C limit settings. "nqc", "pqc", "eqc", "cqc" and "tot" are count arrays for C the test result daily output file, "qc" is the QC flag array. C Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C C VERSION 2 (August, 2012): c c Changed DATE AND TIME in output files to YYYYMMDD and hhmm formats. c C Changed to separate allowable QC flag value between max and min to address setting C "-2 W/m^2" to "-9999.0" at "max allowable" setting of "2" value for SW variables. C "qclim" is for all except SW min, "SWlim" is for SW min c c Added "full correction" for IR loss for both GSW and Diffuse. Changed to subroutine c "IRLossCorr" and added extra lines in the config file for including full correction c moist and dry coefficients in addition to the detector-only coefficients. Routine c sends to subroutine if IR loss correction flag set in config, sends back the full c corrected value if possible, else detector-only corrected, else no correction. Also c returns amout of correction value, and flag on whether full dry/moist (1/2), detector c dry/moist (3/4), or no correction (0). Added these variables to output files. c c Moved testing and trying to calculate gsw/dif/dir if missing from other 2 until after c GSW and Dif IR loss correction sections. c c Added "best estimate" SWdn to output. For SZA<90, is sumSW if available, else GSW. c For SZA>90, is greater (least negative) of SumSW or GSW, if both pass QC tests at c level 2 or less, else whichever has passed QCflagging. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC program qcrad2 INTEGER i,yr,mn,dy,hr,minit,date,time,IP5,IP6,IP7,IP8,nh,aflg INTEGER qc(20),nqc(20),pqc(20),eqc(20),tot(20),tflg,timtyp INTEGER a1,a2,a3,a4,a5,a6,a7,a8,a9,t1,t2,d1,d2,d3,dattyp,k(100) INTEGER a10,a11,a12,a13,tdflg,cqc(20),qclim,qflg,e(100),enum,SWlim integer numint,numrel,dirflg,dmin,skip,dyflg,totn,en,sn,j,trflg integer month(12),p,doy,irl1,irl2,rhf,pflg,dnflg,lastim,lasdate INTEGER dflg,gflg REAL*8 x(0:100),bdat,gsw,dif,dir,dirN,swup,lwdn,lwup,ta,clrsw REAL*8 sigma,pi,gx(5),S,ucz,ncdxtim,rh,prs,plim(20),cswa,cswb REAL*8 elim(20),nlim(20),Sa,clim(20),lwdnTc,lwdnTd,ttmp REAL*8 lwupTc,lwupTd,Tslim,cf(24),badval,tmp,tt1,tt2,cf2(24) REAL*8 rtim,Z,AZ,DEC,Au,hrang,slat,slon,xlong,gt,Tu REAL*8 r1,r2,r3,r4,r5,r6,rlim,uprs,dprs,minprs,irc1m,irc1d,irf REAL*8 irc2d,irc2m,difC,gswC,Cgsw,Cdif,rgsw,rdif REAL*8 irt2m,ird2m,irt2d,ird2d,besw REAL*8 irt1m,ird1m,irt1d,ird1d,gtmp CHARACTER*12 dirfil,cnffil,qcfil,allfil, line CHARACTER*80 infil,outfil CHARACTER*309 qcline CHARACTER*360 limlin CHARACTER*220 clmlin CHARACTER*200 path CHARACTER*3 ext CHARACTER*4 qext CHARACTER*1 ans CHARACTER*8 hdr(60),cdate LOGICAL OK ,dodn, endit DATA MONTH /31,28,31,30,31,30,31,31,30,31,30,31/ C Set some values for the code: bdat=-998.0 badval=-9999.0 x(0)=badval sigma=0.0000000567d0 pi=dacos(-1.d0) IP7 = 0 IP6 = 0 IP8 = 0 s=1368.d0 dodn=.false. endit=.false. lastim=-9999 lasdate=-9999 C coefficients for Rayleigh limit calculation: r1=209.3 r2=-708.3 r3=1128.7 r4=-911.2 r5=287.85 r6=0.046725 plim(1)=-4.0 plim(3)=-4.0 plim(5)=-4.0 plim(7)=-4.0 plim(9)=40.0 plim(10)=700.0 plim(11)=40.0 plim(12)=900.0 elim(1)=-2.0 elim(3)=-2.0 elim(5)=-2.0 elim(7)=-2.0 elim(9)=60.0 elim(10)=500.0 elim(11)=60.0 elim(12)=700.0 cnffil='qc_rad2.cnf' dirfil='qc_rad2.dir' qcfil='qc_rad2.asc' allfil='qc_rad2.all' OPEN(UNIT=9,FILE=qcfil,STATUS='unknown') WRITE(9,*)' date pqc1 pqc2 pqc3 pqc4 pqc5 pqc6 pqc7 pqc %8 pqc9 pqc10 pqc11 pqc12 eqc1 eqc2 eqc3 eqc4 eqc5 eqc6 eqc %7 eqc8 eqc9 eqc10 eqc11 eqc12 nqc1 nqc2 nqc3 nqc4 nqc5 nqc %6 nqc7 nqc8 nqc9 nqc10 cqc1 cqc2 cqc3 cqc4 cqc5 cqc6 cqc %7 cqc8 cqc9 cqc10 cqc11 cqc12 cqc13 cqc14 cqc15 cqc16 cqc17 cqc1 %8 cqc19 tq1 tq2 tq3 tq4 tq5 tq6 tn1 tn2 tn4 tn5& &6 tn7&8 tn9&10 tc13 totn' C Read in configuration settings from file: OPEN(UNIT=1,FILE=cnffil,STATUS='old',ERR=98) OPEN(UNIT=19,FILE='QC_Rad2.cnfread',STATUS='unknown') READ(1,101) path 101 format (a200) i=INDEX(path,'*') do p=(i-1),1,-1 IF(path(p:p).eq.' ') then WRITE(6,*) path(1:p)//'<<' else GO TO 102 endif end do 102 WRITE(6,*)' path= <',path(1:p),'>' WRITE(19,*)' path= <',path(1:p),'>' read(1,*) slat WRITE(19,*)' slat= ', slat read(1,*) XLONG WRITE(19,*)' xlong= ', XLONG read(1,*) nh WRITE(19,*)' nh= ', nh read(1,*) numint WRITE(19,*)' numint= ',numint IF(numint.gt.100) then WRITE(6,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' WRITE(6,*)'! Number of integers in line set > max of 100 !' WRITE(6,*)'! Config setting =', numint,' !' WRITE(6,*)'! Exiting Program... !' WRITE(6,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' pause GO TO 100 endif read(1,*) numrel WRITE(19,*)' numrel= ',numrel IF(numrel.gt.100) then WRITE(6,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' WRITE(6,*)'! Number of reals in line set > max of 100 !' WRITE(6,*)'! Config setting =', numrel,' !' WRITE(6,*)'! Exiting Program... !' WRITE(6,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' pause GO TO 100 endif read(1,*) dattyp WRITE(19,*)' dattyp= ',dattyp if(dattyp.eq.1) then read(1,*) d1,d2,d3 WRITE(19,*)' d1,d2,d3= ',d1,d2,d3 elseif(dattyp.eq.4) then read(1,*) d1,d2 WRITE(19,*)' d1,d2= ',d1,d2 else read(1,*) d1 WRITE(19,*)' d1= ', d1 endif read(1,*) timtyp WRITE(19,*)' timtyp= ',timtyp if(timtyp.eq.1) then read(1,*) t1,t2 WRITE(19,*)' t1,t2= ',t1,t2 ELSE read(1,*) t1 WRITE(19,*)' t1= ', t1 endif read(1,*) ip5 WRITE(19,*)' ip5= ', ip5 read(1,*) a1 WRITE(19,*)' a1= ', a1 read(1,*) a2 WRITE(19,*)' a2= ', a2 read(1,*) a3 WRITE(19,*)' a3= ', a3 read(1,*) a4 WRITE(19,*)' a4= ', a4 read(1,*) a5 WRITE(19,*)' a5= ', a5 read(1,*) a6 WRITE(19,*)' a6= ', a6 read(1,*) a7 WRITE(19,*)' a7= ', a7 read(1,*) a8,rhf WRITE(19,*)' a8,rhf= ', a8,rhf read(1,*) a9,pflg WRITE(19,*)' a9,pflg= ', a9,pflg read(1,*) a10 WRITE(19,*)' a10= ', a10 read(1,*) a11 WRITE(19,*)' a11= ', a11 read(1,*) a12 WRITE(19,*)' a12= ', a12 read(1,*) a13 WRITE(19,*)' a13= ', a13 read(1,*) tflg WRITE(19,*)' tflg= ', tflg read(1,*) tdflg WRITE(19,*)' tdflg= ', tdflg read(1,*) TsLim WRITE(19,*)' tslim= ', tslim read(1,*) dirflg WRITE(19,*)' dirflg= ', dirflg read(1,*) skip WRITE(19,*)' skip= ', skip read(1,*) dyflg,dnflg WRITE(19,*)' dyflg,dnflg= ', dyflg,dnflg read(1,*) aflg WRITE(19,*)' aflg= ', aflg read(1,*) bdat WRITE(19,*)' bdat= ', bdat read(1,*) qclim, SWlim WRITE(19,*)' qclim, SWlim= ', qclim, SWlim read(1,*) qflg,trflg WRITE(19,*)' qflg= ', qflg WRITE(19,*)' trflg= ', trflg do i=1,21 IF(i.le.17) then READ(1,*) cf(i),cf2(i) WRITE(19,*)' i,cf(i),cf2(i)= ', i,cf(i),cf2(i) IF(i.eq.5.or.i.eq.7.or.i.eq.11) then IF(cf(i).lt.cf2(i)) then WRITE(6,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' WRITE(6,*)' WARNING: Lim1 < Lim2',i,cf(i),cf2(i) WRITE(6,*)' Terminating this run.....' endit=.true. endif ELSEIF(i.lt.17) then IF(cf(i).gt.cf2(i)) then WRITE(6,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' WRITE(6,*)' WARNING: Lim1 > Lim2',i,cf(i),cf2(i) WRITE(6,*)' Terminating this run.....' endit=.true. endif endif else READ(1,*) cf(i) WRITE(19,*)' i,cf(i)= ', i,cf(i) endif end do IF(endit) then WRITE(6,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' pause GO TO 100 endif READ(1,*) cswa WRITE(19,*)' cswa= ', cswa READ(1,*) cswb WRITE(19,*)' cswb= ', cswb READ(1,*) minprs, dprs WRITE(19,*)' minprs, dprs= ', minprs, dprs IF(minprs.lt.300.0.or.dprs.lt.500.0) then WRITE(6,*)'%%%%%%%%%%%%%%%%% ERROR %%%%%%%%%%%%%%%%%%%%' WRITE(6,*)' minprs < 500 or default prs < 700:' WRITE(6,*)' minprs=',minprs WRITE(6,*)' dprs=',dprs pause GO TO 100 endif READ(1,*) irl1,irc1m,irc1d WRITE(19,103)' irl1,irc1m,irc1d= ', irl1,irc1m,irc1d 103 FORMAT(a, i4, 2f12.6) 104 FORMAT(a, 2f12.6) READ(1,*) irt1m,ird1m,irt1d,ird1d WRITE(19,104) 'Dif ir loss full corr moist: ',irt1m,ird1m WRITE(19,104) 'Dif ir loss full corr dry: ',irt1d,ird1d READ(1,*) irl2,irc2m,irc2d WRITE(19,103)' irl2,irc2m,irc2d= ', irl2,irc2m,irc2d READ(1,*) irt2m,ird2m,irt2d,ird2d WRITE(19,104) 'GSW ir loss full corr moist: ',irt2m,ird2m WRITE(19,104) 'GSW ir loss full corr dry: ',irt2d,ird2d READ(1,*) enum WRITE(19,*)' enum= ', enum IF(enum.gt.0) then do i=1,enum READ(1,1,END=4) line 1 FORMAT(a12) j=INDEX(line,' ') READ(line(1:j),*) e(i) line=line(j:12) 2 continue IF(line(1:1).eq.' ') then line=line(2:12) GO TO 2 endif READ(line,3) hdr(i) 3 FORMAT(a8) WRITE(19,*)' i,e(i),hdr(i)= ', i,e(i),' ',hdr(i) end do endif 4 CLOSE(UNIT=1) CLOSE(UNIT=19) C Set some test limit variables from the config settings: clim(5)=cf(5) clim(6)=cf(6) clim(7)=cf(7) clim(8)=cf(8) clim(17)=cf(17) clim(18)=cf2(17) clim(19)=cf(18) clim(20)=cf(19) IF(TsLim.gt.bdat) TsLim = TsLim + 273.15 IF(TsLim.lt.170.0.or.TsLim.gt.350.0) TsLim=badval cf(20)=cf(20)+273.15 cf(21)=cf(21)+273.15 elim(9)=cf2(5) elim(10)=cf2(6) elim(11)=cf2(7) elim(12)=cf2(8) IF(cf2(1).gt.1.5000) cf2(1)=1.5000 IF(cf2(2).gt.0.95000) cf2(2)=0.9500 IF(cf2(3).gt.1.000) cf2(3)=1.000 IF(cf2(4).gt.1.2000) cf2(4)=1.2000 IF(cf2(5).lt.40.0000) cf2(5)=40.0000 IF(cf2(6).gt.700.0000) cf2(6)=700.0000 IF(cf2(7).lt.40.0000) cf2(7)=40.0000 IF(cf2(8).gt.900.0000) cf2(8)=900.0000 IF(cf2(9).gt.1.0000) cf2(9)=1.0000 IF(cf2(10).gt.1.1000) cf2(10)=1.1000 IF(cf2(11).lt.0.2000) cf2(11)=0.2000 IF(cf2(12).gt.50.0000) cf2(12)=50.0000 IF(cf2(13).gt.25.0000) cf2(13)=25.0000 IF(cf2(14).gt.35.0000) cf2(14)=35.0000 IF(cf2(15).gt.400.0000) cf2(15)=400.0000 IF(cf2(16).gt.50.0000) cf2(16)=50.0000 C Set header lines for output files: qcline=' Date Time BESW GSW DIF DIR % SWup LWdn LWup Ta RH Prs LWd %Tc LWdTd LWuTc LWuTd qc1 qc2 qc3 qc4 qc5 qc6 qc %7 qc8 qc9 qc10 qc11 qc12 qc13 qc14 qc15 qc16 qc17 qc18 qc19 gflg % dflg Z AU ClrSW DifCorr GSWCorr' limlin=' PGSWmin PGSWmax PDifmin PDifmax PDirmin PDirm %ax PSWupmin PSWupmax PLWdnmin PLWdnmax PLWupmin PLWupmax E %GSWmin EGSWmax EDifmin EDifmax EDirmin EDirmax ESWupmin % ESWupmax ELWdnmin ELWdnmax ELWupmin ELWupmax ESWupmxn ESWu %pmxs G/S min G/S max DRmax Swup/Sum LWdnTmin LWdnTmax %LwupTmin LwupTmax Lup-Lm2 Lup+Lm2' clmlin=' CGSWmax CDifmax CDirmax CSWupmax CLWdnmin CLWdnm %ax CLWupmin CLWupmax CSWupmxn CSWupmxs CLdnTmin CLdnTmax CL %upTmin CLupTmax Lup-Clim Lup+Clim TacdLdn TacdLup Tc-TdMin % Tc-TdMax AvgLWT RLim' C Test longitude setting from config file: IF(XLONG .GE. 0.) THEN IF(XLONG .GT. 180.) THEN SLON= 360.-XLONG ELSE SLON= -XLONG ENDIF ELSE IF(XLONG .LT. -180.) THEN SLON= -(360.+XLONG) ELSE SLON= -XLONG ENDIF C Open "all file" if config set and write header: IF(aflg.gt.0) then OPEN(UNIT=8,FILE=allfil,STATUS='unknown') IF(aflg.eq.1) then WRITE(8,8) qcline,(hdr(i),i=1,enum) else WRITE(8,18) qcline//limlin//clmlin,(hdr(i),i=1,enum) endif endif C Open "directory file" containing list of files to be processed: OPEN(UNIT=1,FILE=dirfil,STATUS='old',ERR=99) C Main program loop, read in input file name, construct output file name, and open files: 5 READ(1,6,END=100) infil 6 FORMAT(a80) OPEN(UNIT=2,FILE=path(1:p)//infil,STATUS='old') WRITE(6,*)' Opened Infile: ',path(1:p)//infil IF(nh.gt.0) then do i=1,nh READ(2,*) END do endif IF(dyflg.gt.0) then do i=80,1,-1 IF(infil(i:i).eq.'.') GO TO 7 end do 7 en=i+1 ext=infil(en:en+2) qext='q'//ext IF(dnflg.gt.0) then dodn=.true. else dodn=.false. sn=INDEX(infil,'.') outfil=infil(1:sn)//qext OPEN(UNIT=3,FILE=outfil,STATUS='unknown') IF(dyflg.eq.1) then WRITE(3,8) qcline,(hdr(i),i=1,enum) else WRITE(3,18) qcline//limlin//clmlin,(hdr(i),i=1,enum) endif endif endif 8 FORMAT(a,60(4x,a8)) 18 FORMAT(a,1x,60(4x,a8)) C Clear daily count arrays: do i=1,20 nqc(i)=0 pqc(i)=0 eqc(i)=0 tot(i)=0 cqc(i)=0 end do totn=0 C Main individual input file loop, read data: 10 IF(skip.eq.1) then READ(2,*,END=50,ERR=50) (k(i),i=1,numint),(x(i),i=1,numrel) else READ(2,*,END=50) (k(i),i=1,numint),(x(i),i=1,numrel) endif totn=totn+1 C Initially set all QC flags to "OK": do i=1,20 qc(i)=0 end do C Set date variables per config settings: IF(dattyp.eq.0) then date=k(d1) IF(date.lt.19000000.or.date.gt.21000000) GO TO 96 yr=date/10000 mn=date/100 - yr*100 dy=date-yr*10000-mn*100 ELSEIF(dattyp.eq.1) then yr=k(d1) mn=k(d2) dy=k(d3) date=yr*10000+mn*100+dy IF(date.lt.19000000.or.date.gt.21000000) GO TO 96 ELSEIF(dattyp.eq.2) then date=FLOOR(x(d1)) yr=date/10000 mn=date/100 - yr*100 dy=date-yr*10000-mn*100 IF(date.lt.19000000.or.date.gt.21000000) GO TO 96 ELSEIF(dattyp.eq.3) then date=k(d1) yr=date/1000 doy=date-(yr*1000) call DOY2DATE(yr,doy,date) IF(date.lt.19000000.or.date.gt.21000000) GO TO 96 mn=date/100 - yr*100 dy=date-yr*10000-mn*100 ELSEIF(dattyp.eq.4) then yr=k(d1) doy=k(d2) call DOY2DATE(yr,doy,date) IF(date.lt.19000000.or.date.gt.21000000) GO TO 96 mn=date/100 - yr*100 dy=date-yr*10000-mn*100 else 11 WRITE(6,*)' Problem reading date! date, Dattyp, d1, d2 d3: ', %date,dattyp,d1,d2,d3 IF(dattyp.eq.2) WRITE(6,*)' X(D1):',X(D1) pause GO TO 100 endif c Check for leap year [for "yr" = 4-digit year] IF(MOD(yr,4).ne.0) then month(2)=28 ELSEIF(MOD(yr,400).eq.0) then month(2)=29 ELSEIF(MOD(yr,100).eq.0) then month(2)=28 ELSE month(2)=29 endif C Set time variables per config settings: IF(timtyp.eq.0) then time=k(t1) hr=time/100 minit=time-hr*100 dmin=hr*60+minit IF(dmin.ge.1440) then dmin=dmin-1440 hr=(dmin)/60 minit=dmin-(hr*60) time=(hr*100)+minit dy=dy+1 IF(dy.gt.month(mn)) then dy=1 mn=mn+1 IF(mn.gt.12) then mn=1 yr=yr+1 endif endif date=yr*10000+mn*100+dy endif IF(time.lt.0.or.time.gt.2359) GO TO 97 elseIF(timtyp.eq.4) then time=k(t1) hr=time/10000 minit=(time/100)-(hr*100) time=hr*100+minit dmin=hr*60+minit IF(dmin.ge.1440) then dmin=dmin-1440 hr=(dmin)/60 minit=dmin-(hr*60) time=(hr*100)+minit dy=dy+1 IF(dy.gt.month(mn)) then dy=1 mn=mn+1 IF(mn.gt.12) then mn=1 yr=yr+1 endif endif date=yr*10000+mn*100+dy endif IF(time.lt.0.or.time.gt.2359) GO TO 97 ELSEIF(timtyp.eq.1) then hr=k(t1) minit=k(t2) time=hr*100+minit dmin=hr*60+minit IF(dmin.ge.1440) then dmin=dmin-1440 hr=(dmin)/60 minit=dmin-(hr*60) time=(hr*100)+minit dy=dy+1 IF(dy.gt.month(mn)) then dy=1 mn=mn+1 IF(mn.gt.12) then mn=1 yr=yr+1 endif endif date=yr*10000+mn*100+dy endif IF(time.lt.0.or.time.gt.2359) GO TO 97 ELSEIF(timtyp.eq.2) then dmin=k(t1) hr=(dmin)/60 minit=dmin-(hr*60) time=(hr*100)+minit IF(time.gt.2359) then dmin=dmin-1440 hr=(dmin)/60 minit=dmin-(hr*60) time=(hr*100)+minit dy=dy+1 IF(dy.gt.month(mn)) then dy=1 mn=mn+1 IF(mn.gt.12) then mn=1 yr=yr+1 endif endif date=yr*10000+mn*100+dy endif IF(time.lt.0.or.time.gt.2359) GO TO 97 ELSEIF(timtyp.eq.3) then ncdxtim=x(d1)-date time=NINT(ncdxtim*10000) hr=time/100 minit=time-(hr*100) dmin=hr*60+minit hr=(dmin)/60 minit=dmin-(hr*60) time=(hr*100)+minit IF(dmin.ge.1440.or.time.gt.2359) then dmin=dmin-1440 hr=(dmin)/60 minit=dmin-(hr*60) time=(hr*100)+minit dy=dy+1 IF(dy.gt.month(mn)) then dy=1 mn=mn+1 IF(mn.gt.12) then mn=1 yr=yr+1 endif endif date=yr*10000+mn*100+dy endif IF(time.lt.0.or.time.gt.2400) GO TO 97 else 12 WRITE(6,*)' Problem reading time! date,time,timtyp,t1,t2: ', %date,time,timtyp,t1,t2 IF(dattyp.eq.2) WRITE(6,*)' X(D1):',X(D1) pause GO TO 100 endif IF(date.eq.lasdate.and.time.le.lastim) then GO TO 10 else lasdate=date lastim=time endif IF(dodn) then dodn=.false. OPEN(UNIT=23,FILE='tmp_date.asc',STATUS='unknown') WRITE(23,13) date 13 FORMAT(i8) CLOSE(unit=23) OPEN(UNIT=23,FILE='tmp_date.asc',STATUS='old') read(23,14) cdate 14 FORMAT(a8) CLOSE(unit=23) outfil=cdate//'.'//qext OPEN(UNIT=3,FILE=outfil,STATUS='unknown') WRITE(6,*)' Opened outfile: ',outfil IF(dyflg.eq.1) then WRITE(3,8) qcline,(hdr(i),i=1,enum) else WRITE(3,18) qcline//limlin//clmlin,(hdr(i),i=1,enum) endif endif C Calculate solar zenith angle and Earth-Sun distance: rtim=time*1.d0 IF(rtim.gt.59.59d0.and.rtim.lt.101.00d0.and.ip6.eq.0) then gt=101.01d0 CALL EPHEMS(SLAT,SLON,dy,mn,yr,GT,IP5,IP6,IP7,IP8, + AZ,AU,Z,HRANG,DEC) Gx(1)=az Gx(2)=au Gx(3)=Z Gx(4)=dec gt=59.59d0 CALL EPHEMS(SLAT,SLON,dy,mn,yr,GT,IP5,IP6,IP7,IP8, + AZ,AU,Z,HRANG,DEC) gt=(rtim-100.d0)+.01 gt=gt/0.61d0 az=((Gx(1)-az)*gt)+az au=((Gx(2)-au)*gt)+au Z=((Gx(3)-Z)*gt)+Z dec=((Gx(4)-dec)*gt)+dec else CALL EPHEMS(SLAT,SLON,dy,mn,yr,rtim,IP5,IP6,IP7,IP8, + AZ,AU,Z,HRANG,DEC) endif cosz=dCOS(Z*PI/180.d0) ucz=cosz IF(ucz.lt.0.d0) ucz=0.d0 Sa=S/AU**2 C Set radiation and met variables per config settings: IF(a1.ge.0) then gsw=x(a1) else gsw=badval endif IF(a2.ge.0) then dif=x(a2) else dif=badval endif IF(dirflg.eq.0.and.a3.ge.0) then dirN=x(a3) IF(dirN.gt.bdat.and.ucz.gt.0.0) then dir=dirn*ucz elseIF(dirN.gt.bdat) then dir=dirn else dir=badval endif ELSEIF(a3.ge.0) then dir=x(a3) IF(dir.le.bdat) then dir=badval endif else dir=badval endif swup=x(a4) lwdn=x(a5) lwup=x(a6) ta=x(a7) rh=x(a8) if (rhf.eq.1) rh=rh*100.0 IF(rh.ge.0.0) then IF(rh.gt.110.0) then rh=badval elseIF(rh.gt.100.0) then rh=100.0 endif ELSEIF(rh.GE.-5.0) then rh=0.0 else rh=badval endif prs=x(a9) IF(pflg.eq.1) prs=prs*10.0 IF(prs.gt.minprs.and.prs.lt.1100.0) then uprs=prs else uprs=dprs endif lwdnTc=x(a10) lwdnTd=x(a11) lwupTc=x(a12) lwupTd=x(a13) C Test temperatures: IF(a7.eq.0) then Ta=badval qc(19) = -1 elseIF(ta.lt.bdat) then ta=badval else tot(13)=tot(13)+1 endif IF(tflg.ne.0.and.a7.ne.0.and.ta.gt.bdat) ta = ta + 273.15 IF(ta.lt.cf(20).or.ta.gt.cf(21)) then IF(a7.ne.0.and.ta.gt.0.0) then qc(19) = 3 cqc(19) = cqc(19) + 1 endif ta=badval endif IF(a10.eq.0) lwdnTc=badval IF(a11.eq.0) lwdnTd=badval IF(a12.eq.0) lwupTc=badval IF(a13.eq.0) lwupTd=badval IF(tdflg.ne.0) then IF(a10.ne.0.and.lwdnTc.gt.bdat) lwdnTc = lwdnTc + 273.15 IF(lwdnTc.lt.cf(20).or.lwdnTc.gt.cf(21)) lwdnTc=badval IF(a11.ne.0.and.lwdnTd.gt.bdat) lwdnTd = lwdnTd + 273.15 IF(lwdnTd.lt.cf(20).or.lwdnTd.gt.cf(21)) lwdnTd=badval IF(a12.ne.0.and.lwupTc.gt.bdat) lwupTc = lwupTc + 273.15 IF(lwupTc.lt.cf(20).or.lwupTc.gt.cf(21)) lwupTc=badval IF(a13.ne.0.and.lwupTd.gt.bdat) lwupTd = lwupTd + 273.15 IF(lwupTd.lt.cf(20).or.lwupTd.gt.cf(21)) lwupTd=badval else IF(a10.ne.0.and.lwdnTc.lt.cf(20).or.lwdnTc.gt.cf(21)) %lwdnTc=badval IF(a11.ne.0.and.lwdnTd.lt.cf(20).or.lwdnTd.gt.cf(21)) %lwdnTd=badval IF(a12.ne.0.and.lwupTc.lt.cf(20).or.lwupTc.gt.cf(21)) %lwupTc=badval IF(a13.ne.0.and.lwupTd.lt.cf(20).or.lwupTd.gt.cf(21)) %lwupTd=badval endif tt1=100.0 tt2=100.0 IF(lwdnTc.gt.0.0.and.lwdnTd.gt.0.0) tt1=ABS(lwdntc-lwdntd) IF(lwupTc.gt.0.0.and.lwupTd.gt.0.0) tt2=ABS(lwuptc-lwuptd) cc calculate TTmp ttmp=0.0 i=0 IF(lwdnTc.gt.0.0.and.tt1.lt.10.0) then ttmp=ttmp+lwdnTc i=i+1 endif IF(lwdnTd.gt.0.0.and.tt1.lt.10.0) then ttmp=ttmp+lwdnTd i=i+1 endif IF(lwupTc.gt.0.0.and.tt2.lt.10.0) then ttmp=ttmp+lwupTc i=i+1 endif IF(lwupTd.gt.0.0.and.tt2.lt.10.0) then ttmp=ttmp+lwupTd i=i+1 endif IF(i.gt.0) then ttmp=ttmp/(1.0*i) else ttmp=-9.0 endif cc test Ta vs TTmp IF(ttmp.gt.0.0.and.ta.gt.0.0) then IF(ABS(ta-ttmp).gt.20.0) then ta=badval qc(19)=4 cqc(19) = cqc(19) + 1 endif endif cc test LWdn & LWup Tc & Td vs TTmp ok=.true. IF(ttmp.gt.0.0.and.lwdntc.gt.0.0) then IF(lwdntc.lt.ttmp-15.0) then lwdntc=badval qc(13)=5 IF(ok) cqc(17)=cqc(17)+1 IF(ok) OK=.false. elseIF(lwdntc.gt.ttmp+15.0) then lwdntc=badval qc(13)=6 IF(ok) cqc(17)=cqc(17)+1 IF(ok) OK=.false. endif endif IF(ttmp.gt.0.0.and.lwdntd.gt.0.0) then IF(lwdntd.lt.ttmp-15.0) then lwdntd=badval qc(14)=5 IF(ok) cqc(17)=cqc(17)+1 IF(ok) OK=.false. elseIF(lwdntd.gt.ttmp+15.0) then lwdntd=badval qc(14)=6 IF(ok) cqc(17)=cqc(17)+1 IF(ok) OK=.false. endif endif IF(ttmp.gt.0.0.and.lwuptc.gt.0.0) then IF(lwuptc.lt.ttmp-15.0) then lwuptc=badval qc(15)=5 IF(ok) cqc(17)=cqc(17)+1 IF(ok) OK=.false. elseIF(lwuptc.gt.ttmp+15.0) then lwuptc=badval qc(15)=6 IF(ok) cqc(17)=cqc(17)+1 IF(ok) OK=.false. endif endif IF(ttmp.gt.0.0.and.lwuptd.gt.0.0) then IF(lwuptd.lt.ttmp-15.0) then lwuptd=badval qc(16)=5 IF(ok) cqc(17)=cqc(17)+1 IF(ok) OK=.false. elseIF(lwuptd.gt.ttmp+15.0) then lwuptd=badval qc(16)=6 IF(ok) cqc(17)=cqc(17)+1 IF(ok) OK=.false. endif endif cc calculate LWdn detector flux if possible IF(lwdntc.gt.0.0.and.lwdntd.gt.0.0.and.lwdn.gt.0.0) then irf=lwdn-(sigma*lwdntc**4)+(4.0*sigma*(lwdntd**4-lwdntc**4)) else irf=badval endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cc adjust global SW for generic IR loss if flag set IF(lwdnTc.gt.0.0) then Tu=lwdnTc ELSEIF(Ta.gt.0.0) then Tu=Ta else Tu=-9999.0 endif rgsw=gsw gswC=0.0 Cgsw=0.0 gflg=0 IF(irl2.eq.1.and.gsw.GT.bdat) then call IRLossCorr(rgsw,LWdn,RH,Tu,lwdnTc,lwdnTd,irf,z,irt2m,ird2m, %irt2d,ird2d,irc2m,irc2d,Cgsw,gswC,gflg) gsw=Cgsw endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cc adjust diffuse SW for generic IR loss if flag set difC=0.0 Cdif=0.0 dflg=0 rdif=dif IF(irl1.eq.1.and.dif.gt.bdat) then call IRLossCorr(rdif,LWdn,RH,Tu,lwdnTc,lwdnTd,irf,z,irt1m,ird1m, %irt1d,ird1d,irc1m,irc1d,Cdif,difC,dflg) dif=Cdif endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc IF(gsw.lt.bdat.and.dir.GT.-2.0.and.dif.GT.-2.0) gsw=dir+dif IF(dif.lt.bdat.and.dir.GT.-2.0.and.gsw.GT.-2.0) dif=gsw-dir IF(dir.lt.bdat.and.dif.GT.-2.0.and.gsw.GT.-2.0) dir=gsw-dif C Limit calculations and Measurements tests: plim(2)=(Sa*1.5*ucz**1.2+100.0) cc elim(2)=(Sa*1.2*ucz**1.2+50.0) elim(2)=(Sa*cf2(1)*ucz**1.2+55.0) clim(1)=(Sa*cf(1)*ucz**1.2+50.0) IF(gsw.gt.bdat) then tot(1)=tot(1)+1 IF(gsw.LT.plim(1)) then qc(1)=5 pqc(1)=pqc(1)+1 elseIF(gsw.LT.elim(1)) then qc(1)=3 eqc(1)=eqc(1)+1 endif IF(gsw.GT.plim(2)) then qc(1)=6 pqc(2)=pqc(2)+1 elseIF(gsw.GT.elim(2)) then qc(1)=4 eqc(2)=eqc(2)+1 elseIF(gsw.GT.clim(1)) then qc(1)=2 cqc(1)=cqc(1)+1 endif else qc(1)=-1 gsw=badval endif plim(4)=(Sa*0.95*ucz**1.2+50.0) cc elim(4)=(Sa*0.75*ucz**1.2+30.0) elim(4)=(Sa*cf2(2)*ucz**1.2+35.0) clim(2)=(Sa*cf(2)*ucz**1.2+30.0) rlim=r1*ucz+r2*ucz**2+r3*ucz**3+r4*ucz**4+r5*ucz**5+r6*ucz*uprs IF(dif.gt.bdat) then tot(2)=tot(2)+1 IF(gsw.gt.50.0.and.(dif/gsw).lt.0.8.and.dif.LT.(rlim-1.0)) then qc(2)=8 pqc(3)=pqc(3)+1 elseIF(dif.LT.plim(3)) then qc(2)=5 pqc(3)=pqc(3)+1 elseIF(dif.LT.elim(3)) then qc(2)=3 eqc(3)=eqc(3)+1 endif IF(dif.GT.plim(4)) then qc(2)=6 pqc(4)=pqc(4)+1 elseIF(dif.GT.elim(4)) then qc(2)=4 eqc(4)=eqc(4)+1 elseIF(dif.GT.clim(2)) then qc(2)=2 cqc(2)=cqc(2)+1 endif else qc(2)=-1 dif=badval endif plim(6)=(Sa) cc elim(6)=(Sa*0.95*ucz**1.2+10.0) elim(6)=(Sa*cf2(3)*ucz**1.2+15.0) clim(3)=(Sa*cf(3)*ucz**1.2+10.0) IF(dir.gt.bdat) then tot(3)=tot(3)+1 IF(dir.LT.plim(5)) then qc(3)=5 pqc(5)=pqc(5)+1 elseIF(dir.LT.elim(5)) then qc(3)=3 eqc(5)=eqc(5)+1 endif IF(dirN.GT.plim(6)) then qc(3)=6 pqc(6)=pqc(6)+1 ELSEIF(dir.GT.elim(6)) then qc(3)=4 eqc(6)=eqc(6)+1 ELSEIF(dir.GT.clim(3)) then qc(3)=2 cqc(3)=cqc(3)+1 endif else qc(3)=-1 dirn=badval dir=badval endif plim(6)=(Sa*ucz+10.0) plim(8)=(Sa*1.2*ucz**1.2+50.0) cc elim(8)=(Sa*ucz**1.2+50.0) elim(8)=(Sa*cf2(4)*ucz**1.2+55.0) clim(4)=(Sa*cf(4)*ucz**1.2+50.0) IF(swup.gt.bdat) then tot(4)=tot(4)+1 IF(swup.LT.plim(7)) then qc(4)=5 pqc(7)=pqc(7)+1 elseIF(swup.LT.elim(7)) then qc(4)=3 eqc(7)=eqc(7)+1 endif IF(swup.GT.plim(8)) then qc(4)=6 pqc(8)=pqc(8)+1 elseIF(swup.GT.elim(8)) then qc(4)=4 eqc(8)=eqc(8)+1 elseIF(swup.GT.clim(4)) then qc(4)=2 cqc(4)=cqc(4)+1 endif else qc(4)=-1 swup=badval endif IF(lwdn.gt.bdat) then tot(5)=tot(5)+1 IF(lwdn.LT.plim(9)) then qc(5)=5 pqc(9)=pqc(9)+1 elseIF(lwdn.LT.elim(9)) then qc(5)=3 eqc(9)=eqc(9)+1 elseIF(lwdn.LT.clim(5)) then qc(5)=1 cqc(5)=cqc(5)+1 endif IF(lwdn.GT.plim(10)) then qc(5)=6 pqc(10)=pqc(10)+1 elseIF(lwdn.GT.elim(10)) then qc(5)=4 eqc(10)=eqc(10)+1 elseIF(lwdn.GT.clim(6)) then qc(5)=2 cqc(6)=cqc(6)+1 endif else qc(5)=-1 lwdn=badval endif IF(lwup.gt.bdat) then tot(6)=tot(6)+1 IF(lwup.LT.plim(11)) then qc(6)=5 pqc(11)=pqc(11)+1 elseIF(lwup.LT.elim(11)) then qc(6)=3 eqc(11)=eqc(11)+1 elseIF(lwup.LT.clim(7)) then qc(6)=1 cqc(7)=cqc(7)+1 endif IF(lwup.GT.plim(12)) then qc(6)=6 pqc(12)=pqc(12)+1 elseIF(lwup.GT.elim(12)) then qc(6)=4 eqc(12)=eqc(12)+1 elseIF(lwup.GT.clim(8)) then qc(6)=2 cqc(8)=cqc(8)+1 endif else qc(6)=-1 lwup=badval endif IF(z.le.75.0) then nlim(1)=0.92 nlim(2)=1.08 elseIF(z.le.93.0) then nlim(1)=0.85 nlim(2)=1.15 else nlim(1)=badval nlim(2)=badval endif IF(qc(1).eq.0.and.qc(2).eq.0.and.qc(3).eq.0.and.z.lt.93.0) then tot(7)=tot(7)+1 IF(z.le.75.0.and.(dif+dir).gt.50.0) then IF(gsw/(dif+dir).lt.nlim(1).or.gsw/(dif+dir).gt.nlim(2)) then qc(7)=1 nqc(1)=nqc(1)+1 endif elseIF((dif+dir).gt.50.0) then IF(gsw/(dif+dir).lt.nlim(1).or.gsw/(dif+dir).gt.nlim(2)) then qc(7)=2 nqc(1)=nqc(1)+1 endif else qc(7)=-1 endif else qc(7)=-1 endif IF(z.le.75.0) then nlim(3)=1.05 elseIF(z.le.93.0) then nlim(3)=1.10 else nlim(3)=badval endif IF(qc(1).eq.0.and.qc(2).eq.0.and.z.lt.93.0) then IF(gsw.gt.50.0) tot(8)=tot(8)+1 IF(z.le.75.0.and.gsw.gt.50.0) then IF((dif/gsw).gt.nlim(3)) then qc(8)=1 nqc(2)=nqc(2)+1 endif ELSEIF(gsw.gt.50.0) then IF((dif/gsw).gt.nlim(3)) then qc(8)=2 nqc(2)=nqc(2)+1 endif else qc(8)=-1 endif else qc(8)=-1 endif IF((dir+dif).gt.50.0) then nlim(4)=(dir+dif) IF(gsw.gt.50.0) then tmp=gsw else tmp=badval endif elseIF(gsw.gt.50.0) then nlim(4)=(gsw) IF(dir+dif.gt.50.0) then tmp=dir+dif else tmp=badval endif else nlim(4)=badval endif IF(ta.gt.bdat.and.tslim.gt.bdat.and.nlim(4).gt.bdat) then clim(9)=cf(9)*nlim(4)+25.0 clim(10)=cf(10)*nlim(4)+25.0 IF(clim(10).gt.nlim(4)) clim(10)=nlim(4) elim(13)=cf2(9)*nlim(4)+30.0 elim(14)=cf2(10)*nlim(4)+30.0 IF(elim(14).gt.nlim(4)+25.0) elim(14)=nlim(4)+25.0 else clim(9)=badval clim(10)=badval elim(13)=badval elim(14)=badval endif IF(qc(2).eq.0.and.qc(3).eq.0.and.qc(4).eq.0) then IF((dir+dif).gt.50.0) then tot(9)=tot(9)+1 IF(swup.GT.nlim(4)+25.0) then qc(9)=5 nqc(4)=nqc(4)+1 IF(tmp.gt.0.0.and.swup.gt.tmp+25.0) qc(9)=6 elseIF(ta.ge.tslim.and.clim(9).gt.bdat) then IF(swup.gt.elim(13)) then qc(9)=3 cqc(9)=cqc(9)+1 elseIF(swup.gt.clim(9)) then qc(9)=1 cqc(9)=cqc(9)+1 endif elseIF(ta.lt.tslim.and.clim(10).gt.bdat) then IF(swup.gt.elim(14)) then qc(9)=4 cqc(10)=cqc(10)+1 elseIF(swup.gt.clim(10)) then qc(9)=2 cqc(10)=cqc(10)+1 endif elseIF(clim(10).lt.bdat.and.nlim(4).gt.bdat) then clim(10)=cf(10)*nlim(4)+25.0 elim(14)=cf2(10)*nlim(4)+30.0 IF(swup.GT.elim(14)) then qc(9)=4 cqc(10)=cqc(10)+1 elseIF(swup.GT.clim(10)) then qc(9)=2 cqc(10)=cqc(10)+1 endif endif else qc(9)=-1 endif elseif(qc(1).eq.0.and.qc(4).eq.0) then IF(gsw.gt.50.0) then tot(9)=tot(9)+1 IF(swup.GT.nlim(4)+25.0) then qc(9)=5 nqc(4)=nqc(4)+1 IF(tmp.gt.0.0.and.swup.gt.tmp+25.0) qc(9)=6 elseIF(ta.ge.tslim.and.clim(9).gt.bdat) then IF(swup.gt.elim(13)) then qc(9)=3 cqc(9)=cqc(9)+1 elseIF(swup.gt.clim(9)) then qc(9)=1 cqc(9)=cqc(9)+1 endif elseIF(ta.lt.tslim.and.clim(10).gt.bdat) then IF(swup.gt.elim(14)) then qc(9)=4 cqc(10)=cqc(10)+1 elseIF(swup.gt.clim(10)) then qc(9)=2 cqc(10)=cqc(10)+1 endif elseIF(clim(10).lt.bdat.and.nlim(4).gt.bdat) then clim(10)=cf(10)*nlim(4)+25.0 elim(14)=cf2(10)*nlim(4)+30.0 IF(swup.gt.elim(14)) then qc(9)=4 cqc(10)=cqc(10)+1 elseIF(swup.gt.clim(10)) then qc(9)=2 cqc(10)=cqc(10)+1 endif endif else qc(9)=-1 endif else qc(9)=-1 endif IF(ta.GT.0.0) then nlim(5)=(cf2(11)*sigma*ta**4) nlim(6)=(sigma*ta**4+cf2(12)) clim(11)=(cf(11)*sigma*ta**4) clim(12)=(sigma*ta**4+cf(12)) else nlim(5)=badval nlim(6)=badval clim(11)=badval clim(12)=badval endif IF(ta.GT.0.0.and.qc(5).eq.0) then tot(10)=tot(10)+1 IF(lwdn.LT.nlim(5)) then qc(10)=3 nqc(5)=nqc(5)+1 elseIF(lwdn.LT.clim(11)) then qc(10)=1 cqc(11)=cqc(11)+1 ELSEIF(lwdn.gT.nlim(6)) then qc(10)=4 nqc(6)=nqc(6)+1 ELSEIF(lwdn.gT.clim(12)) then qc(10)=2 cqc(12)=cqc(12)+1 endif else qc(10)=-1 endif IF(ta.GT.0.0) then nlim(7)=(sigma*(ta-cf2(13))**4) nlim(8)=(sigma*(ta+cf2(14))**4) clim(13)=(sigma*(ta-cf(13))**4) clim(14)=(sigma*(ta+cf(14))**4) else nlim(7)=badval nlim(8)=badval clim(13)=badval clim(14)=badval endif IF(ta.GT.0.0.and.qc(6).eq.0) then tot(11)=tot(11)+1 IF(lwup.LT.nlim(7)) then qc(11)=3 nqc(7)=nqc(7)+1 elseIF(lwup.LT.clim(13)) then qc(11)=1 cqc(13)=cqc(13)+1 ELSEIF(lwup.gT.nlim(8)) then qc(11)=4 nqc(8)=nqc(8)+1 ELSEIF(lwup.gT.clim(14)) then qc(11)=2 cqc(14)=cqc(14)+1 endif else qc(11)=-1 endif IF(LWup.GT.-1.0) then nlim(9)=(lwup-cf2(15)) nlim(10)=(lwup+cf2(16)) clim(15)=(lwup-cf(15)) clim(16)=(lwup+cf(16)) else nlim(9)=badval nlim(10)=badval clim(15)=badval clim(16)=badval endif IF(qc(5).eq.0.and.qc(6).eq.0) then tot(12)=tot(12)+1 IF(lwdn.lt.nlim(9)) then qc(12)=3 nqc(9)=nqc(9)+1 elseIF(lwdn.lt.clim(15)) then qc(12)=1 cqc(15)=cqc(15)+1 elseIF(lwdn.gt.nlim(10)) then qc(12)=4 nqc(10)=nqc(10)+1 elseIF(lwdn.gt.clim(16)) then qc(12)=2 cqc(16)=cqc(16)+1 endif else qc(12)=-1 endif ccccccccccccccccccccccccccccccccccccccccccccccccccccccc OK=.true. IF(lwdnTc.gt.bdat.and.Ta.gt.bdat) then IF(lwdnTc.gt.Ta+cf(17)) then qc(13)=4 IF(ok) cqc(17)=cqc(17)+1 IF(ok) OK=.false. elseIF(lwdnTc.lt.Ta-cf(17)) then qc(13)=3 IF(ok) cqc(17)=cqc(17)+1 IF(ok) OK=.false. endif else qc(13)=-1 endif IF(lwdnTd.gt.bdat.and.Ta.gt.bdat) then IF(lwdnTd.gt.Ta+cf(17)) then qc(14)=4 IF(ok) cqc(17)=cqc(17)+1 IF(ok) OK=.false. elseIF(lwdnTd.lt.Ta-cf(17)) then qc(14)=3 IF(ok) cqc(17)=cqc(17)+1 IF(ok) OK=.false. endif else qc(14)=-1 endif IF(lwupTc.gt.bdat.and.Ta.gt.bdat) then IF(lwupTc.gt.Ta+cf2(17)) then qc(15)=4 IF(ok) cqc(17)=cqc(17)+1 IF(ok) OK=.false. elseIF(lwupTc.lt.Ta-cf2(17)) then qc(15)=3 IF(ok) cqc(17)=cqc(17)+1 IF(ok) OK=.false. endif else qc(15)=-1 endif IF(lwupTd.gt.bdat.and.Ta.gt.bdat) then IF(lwupTd.gt.Ta+cf2(17)) then qc(16)=4 IF(ok) cqc(17)=cqc(17)+1 IF(ok) OK=.false. elseIF(lwupTd.lt.Ta-cf2(17)) then qc(16)=3 IF(ok) cqc(17)=cqc(17)+1 IF(ok) OK=.false. endif else qc(16)=-1 endif IF(qflg.gt.1) then do i=9,18 IF(qc(i).gt.qclim) then IF(i.eq.9) swup=badval IF(i.eq.10) lwdn=badval IF(i.eq.13) lwdn=badval IF(i.eq.14) lwdn=badval IF(i.eq.17) lwdn=badval IF(i.eq.11) lwup=badval IF(i.eq.15) lwup=badval IF(i.eq.16) lwup=badval IF(i.eq.18) lwup=badval endif end do endif ccccccccccccccccccccccccccccccccccccccccccccccccccccccc OK=.true. IF(lwdnTc.gt.bdat.and.lwdnTd.gt.bdat) then IF((lwdnTc-lwdnTd).gt.cf(19)) then qc(17)=4 IF(ok) cqc(18)=cqc(18)+1 IF(ok) OK=.false. elseIF((lwdnTc-lwdnTd).lt.cf(18)) then qc(17)=3 IF(ok) cqc(18)=cqc(18)+1 IF(ok) OK=.false. endif else qc(17)=-1 endif IF(lwupTc.gt.bdat.and.lwupTd.gt.bdat) then IF((lwupTc-lwupTd).gt.cf(19)) then qc(18)=4 IF(ok) cqc(18)=cqc(18)+1 IF(ok) OK=.false. elseIF((lwupTc-lwupTd).lt.cf(18)) then qc(18)=3 IF(ok) cqc(18)=cqc(18)+1 IF(ok) OK=.false. endif else qc(18)=-1 endif ccccccccccccccccccccccccccccccccccccccccccccccccccccccc IF(dir.GT.-90.0.and.dif.GT.-90.0) then tmp=dir+dif ELSEIF(gsw.GT.-90.0) then tmp=gsw else tmp=-999.0 endif IF(ucz.gt.0.0.and.cswa.gt.0.0) then clrsw=(cswa/au**2)*ucz**cswb IF(dif.GT.50.0.and.tmp.GT.-90.0.and.trflg.gt.0) then cc put in a counter to ratio with this test for plotting % of failures cc IF(tmp/clrsw.ge.0.9.and.dif/tmp.ge.0.9) then IF(tmp/clrsw.ge.0.85.and.dif/tmp.ge.0.85) then qc(2)=9 qc(3)=9 nqc(3)=nqc(3)+1 IF(qc(2).gt.qclim) dif=badval IF(qc(3).gt.qclim) dir=badval endif endif ELSE clrsw=badval endif c Set values to "-999.0" if set in config file for definitive tests: c check for min setting versus "SWlim" and all else versus "QClim" for SW variables IF(qflg.gt.0) then do i=1,6 IF(qc(i).gt.qclim) then IF(i.eq.1.and.MOD(i,2).eq.0) gsw=badval IF(i.eq.2.and.MOD(i,2).eq.0) dif=badval IF(i.eq.3.and.MOD(i,2).eq.0) dir=badval IF(i.eq.3.and.MOD(i,2).eq.0) dirN=badval IF(i.eq.4.and.MOD(i,2).eq.0) swup=badval IF(i.eq.5) lwdn=badval IF(i.eq.6) lwup=badval endif IF(qc(i).gt.SWlim) then IF(i.eq.1.and.MOD(i,2).eq.1) gsw=badval IF(i.eq.2.and.MOD(i,2).eq.1) dif=badval IF(i.eq.3.and.MOD(i,2).eq.1) dir=badval IF(i.eq.3.and.MOD(i,2).eq.1) dirN=badval IF(i.eq.4.and.MOD(i,2).eq.1) swup=badval endif end do endif IF(qflg.gt.1) then do i=9,18 IF(qc(i).gt.qclim) then IF(i.eq.9) swup=badval IF(i.eq.10) lwdn=badval IF(i.eq.13) lwdn=badval IF(i.eq.14) lwdn=badval IF(i.eq.17) lwdn=badval IF(i.eq.11) lwup=badval IF(i.eq.15) lwup=badval IF(i.eq.16) lwup=badval IF(i.eq.18) lwup=badval endif end do endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Determine best estimate SWdn c IF(dif.Ge.-2.0.and.dir.Ge.-2.0.and.qc(2).le.2.and.qc(3).le.2) then gtmp=dif+dir else gtmp=-9999.0 endif besw=-9999.0 IF(z.le.90.0) then IF(gtmp.GE.-2.0) then besw= gtmp else besw=gsw endif ELSE IF(gtmp.GE.-2.0) then IF(qc(1).le.2) then IF(gtmp.gt.gsw) then besw=gtmp else besw=gsw endif else besw=gtmp endif else besw=gsw endif endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccc C write data to daily output file if set in config: 21 FORMAT(2i8,14f10.1,21i5,f9.1,f9.5,3f10.1,61f12.4) 22 FORMAT(2i8,14f10.1,21i5,f9.1,f9.5,29f10.1,3f10.3,29f10.1, %61f12.4) IF(dyflg.eq.1) then WRITE(3,21)date,time,besw,GSW,DIF,DIR,SWup,LWdn,LWup,Ta,RH, %Prs,lwdnTc,lwdnTd,lwupTc,lwupTd,(qc(i),i=1,19),gflg,dflg,z,au, %clrsw,difC,gswC,(x(e(i)),i=1,enum) ELSEIF(dyflg.eq.2) then WRITE(3,22)date,time,besw,GSW,DIF,DIR,SWup,LWdn,LWup,Ta,RH, %Prs,lwdnTc,lwdnTd,lwupTc,lwupTd,(qc(i),i=1,19),gflg,dflg,z,au, %clrsw,difC,gswC,(plim(i),i=1,12),(elim(i),i=1,14),(nlim(i),i=1,10) %,(clim(i),i=1,20),ttmp,rlim,(x(e(i)),i=1,enum) endif C write data to "all" output file if set in config: IF(aflg.eq.1) then WRITE(8,21)date,time,besw,GSW,DIF,DIR,SWup,LWdn,LWup,Ta,RH, %Prs,lwdnTc,lwdnTd,lwupTc,lwupTd,(qc(i),i=1,19),gflg,dflg,z,au, %clrsw,difC,gswC,(x(e(i)),i=1,enum) ELSEIF(aflg.eq.2) then WRITE(8,22)date,time,besw,GSW,DIF,DIR,SWup,LWdn,LWup,Ta,RH, %Prs,lwdnTc,lwdnTd,lwupTc,lwupTd,(qc(i),i=1,19),gflg,dflg,z,au, %clrsw,difC,gswC,(plim(i),i=1,12),(elim(i),i=1,14),(nlim(i),i=1,10) %,(clim(i),i=1,20),ttmp,rlim,(x(e(i)),i=1,enum) endif GO TO 10 50 CLOSE(unit=2) CLOSE(unit=3) C write data to daily count output file: WRITE(9,51) date,(pqc(i),i=1,12),(eqc(i),i=1,12),(nqc(i),i=1,10), %(cqc(i),i=1,19),(tot(i),i=1,13),totn 51 FORMAT(i8,100i6) GO TO 5 C ERROR messages to screen: 95 FORMAT(a1) 96 WRITE(6,*)' Problem reading date! date, Dattyp, d1, d2 d3: ', %date,dattyp,d1,d2,d3 IF(dattyp.eq.2) WRITE(6,*)' X(D1):',X(D1) WRITE(6,*)' continue (c), skip file (s), or end program (e)? ' READ(5,95) ans IF(ans.eq.'c'.or.ans.eq.'C') GO TO 10 IF(ans.eq.'s'.or.ans.eq.'S') GO TO 50 GO TO 100 97 WRITE(6,*)' Problem reading time! date,time,timtyp,t1,t2: ', %date,time,timtyp,t1,t2 IF(timtyp.eq.4) WRITE(6,*)' in time: ',k(t1) IF(timtyp.eq.4) WRITE(6,*)' hr,minit: ',hr,minit IF(dattyp.eq.2) WRITE(6,*)' X(D1),ncdxtim:',X(D1),ncdxtim WRITE(6,*)' continue (c), skip file (s), or end program (e)? ' READ(5,95) ans IF(ans.eq.'c'.or.ans.eq.'C') GO TO 10 IF(ans.eq.'s'.or.ans.eq.'S') GO TO 50 GO TO 100 98 WRITE(6,*)'======================================================' WRITE(6,*)'======= Error opening "',cnffil,'" ==========' WRITE(6,*)' >>>>> hit ENTER to end program <<<<<' READ(5,*) GO TO 100 99 WRITE(6,*)'======================================================' WRITE(6,*)'======= Error opening "',dirfil,'" ==========' WRITE(6,*)' >>>>> hit ENTER to end program <<<<<' READ(5,*) GO TO 100 100 CLOSE(UNIT=1) CLOSE(UNIT=9) CLOSE(UNIT=8) end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE DOY2DATE(year,doy,date) C subroutine for converting DOY to YYYYMMDD format integer year,doy,date,MONTHS(12),i,num DATA MONTHS /31,28,31,30,31,30,31,31,30,31,30,31/ c Check for leap year [for "yr" = 4-digit year] IF(MOD(year,4).ne.0) then months(2)=28 ELSEIF(MOD(year,400).eq.0) then months(2)=29 ELSEIF(MOD(year,100).eq.0) then months(2)=28 ELSE months(2)=29 endif c Determine month and day from doy num=doy do i=1,12 num=num-months(i) if(num.le.0) then num=num+months(i) go to 10 endif enddo 10 date=year*10000+i*100+num RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE IRLossCorr(SWdn,LWdn,RH,Ta,Tc,Td,det,z,irt1m,ird1m, %irt1d,ird1d,irc2m,irc2d,SWcorr,corrfac,cflg) c Subroutine to apply IR loss correction to downwelling SW measurements given c pre-determined coefficients. integer cflg Real*8 SWdn,LWdn,RH,Ta,Tc,Td,det,z,irt1m,ird1m,irt1d,ird1d REAL*8 sigma,Te,cgswd,cgswf,dfac,tfac,corrfac,SWcorr,Tflux REAL*8 irc2m,irc2d cc Input variables cc SWdn, The SW variable to be corrected for IR loss cc LWdn, LW downwelling irradiance cc RH, Relative Humidity in % cc Ta, ambient air temperature in K cc Tc, LW pyrgeometer case temperature in K cc Td, LW pyrgeometer dome temperature in K cc det, LW pyrgeometer detector flux in W/m^2 cc z solar zenith angle cc [irt1m,ird1m,irt1d,ird1d] are moist (Tflux, det) and dry (Tflux, det) full corr coeffs (neg, pos) cc [irc2m,irc2d] are the moist and dry detector-only coeffs (positive) cc cc Returned variables: cc SWcorr: corrected SW value cc corrfac: correction factor applied cc cflg: type of correction applied, none/full moist/full dry/detector moist/det dry (0/1/2/3/4) sigma=0.0000000567d0 corrfac=0.00000 SWcorr=0.000000 cc calculate LWdn sky brightness temperature if possible IF(lwdn.gt.0.0) then te=(lwdn/sigma)**0.25 else te=-9999.0 endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cc adjust SW for IR loss cgswd=0.0 cgswf=0.0 cflg=0 IF(tc.gt.10.0.and.td.gt.10.0.and.(tc-td).ge.-1.5.and. %(tc-td).lt.3.0) then Tflux=sigma*(tc**4.d0-td**4.d0) else Tflux=-9999.0 endif IF(z.le.80.0) then dfac=1.4d0 tfac=2.0d0 ELSEIF(z.gt.80.0.and.z.lt.90.0) then dfac=1.d0+((90.d0-z)/10.d0)*(0.4d0) tfac=1.d0+((90.d0-z)/10.d0)*(1.d0) else dfac=1.d0 tfac=1.d0 endif IF(det.GT.-250.0.and.te.gt.0.0.and.ta.gt.0.0) then IF(RH.ge.0.0) then IF(rh.ge.80.0.AND.(ta-te).le.6.0) then cflg=2 cgswd= dfac*(irc2m*det) IF(tflux.GT.-9990.0) then cgswf= tfac*(ird1m*det) + irt1m*tflux else cgswf= -9999.0 endif else cflg=1 cgswd= dfac*(irc2d*det) IF(tflux.GT.-9990.0) then cgswf= tfac*(ird1d*det) + irt1d*tflux else cgswf= -9999.0 endif endif else IF((ta-te).le.6.0) then cflg=2 cgswd= dfac*(irc2m*det) IF(tflux.GT.-9990.0) then cgswf= tfac*(ird1m*det) + irt1m*tflux else cgswf= -9999.0 endif else cflg=1 cgswd= dfac*(irc2d*det) IF(tflux.GT.-9990.0) then cgswf= tfac*(ird1d*det) + irt1d*tflux else cgswf= -9999.0 endif endif endif cccccccccc if no Ta available, do dry correction for all data ccccccccccccccccccc elseIF(det.GT.-250.0) then cflg=1 cgswd= dfac*(irc2d*det) IF(tflux.GT.-9990.0) then cgswf= tfac*(ird1d*det) + irt1d*tflux else cgswf= -9999.0 endif endif cgswd=-cgswd IF(cgswf.GT.-900.0) then cgswf=-cgswf SWcorr=SWdn+cgswf corrfac=cgswf ELSEIF(cgswd.GT.-900.0) then SWcorr=SWdn+cgswd cflg=cflg+2 corrfac=cgswd else SWcorr=-9999.0 cflg=-1 endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C====================================================================== C====================================================================== SUBROUTINE EPHEMS (LAT,LONG,DAY,MONTH,YEAR,TIME,P5,P6,P7,P8, * A,R,Z,LHA,DEC) C C SUBROUTINE TO COMPUTE SOLAR POSITION C C WRITTEN BY WAYNE H WILSON, JR. C VISIBILITY LABORATORY C SCRIPPS INSTITUTION OF OCEANOGRAPHY C UNIVERSITY OF CALIFORNIA, SAN DIEGO C LA JOLLA, CALIFORNIA, 92093 C PH. (714)-294-5534 C 25 FEB 1980 C C C C %%%%%%%%%%%%%%%%% NOTE: KNOWN BUG PATCH %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C Wilson's code produces the incorrect output in the TIME is set C to between 100.00 and 100.59. The following patch loop in the C MAIN program will avoid this problem. Declare REAL*8 gt,x(5) C CCCCCCCCCCCCCCCCCCCCCC PATCH LOOP CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C IF(time.gt.59.59d0.and.time.lt.101.00d0.and.ip6.eq.0) then C gt=101.01d0 C CALL EPHEMS(SLAT,SLON,IDAY,IMONTH,IYEAR,gt,IP5,IP6,IP7,IP8, C + AZANG,RADAU,SOLZA,HRANG,DECANG) C x(1)=azang C x(2)=radau C x(3)=solza C x(4)=decang C gt=59.59d0 C CALL EPHEMS(SLAT,SLON,IDAY,IMONTH,IYEAR,gt,IP5,IP6,IP7,IP8, C + AZANG,RADAU,SOLZA,HRANG,DECANG) C gt=(time-100.d0)+.01 C gt=gt/0.61d0 C azang=((x(1)-azang)*gt)+azang C radau=((x(2)-radau)*gt)+radau C solza=((x(3)-solza)*gt)+solza C decang=((x(4)-decang)*gt)+decang C else C CALL EPHEMS(SLAT,SLON,IDAY,IMONTH,IYEAR,TIME,IP5,IP6,IP7,IP8, C + AZANG,RADAU,SOLZA,HRANG,DECANG) C endif C CCCCCCCCCCCCCCCCCCCC End Patch CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C IMPLICIT REAL*8(P) IMPLICIT INTEGER*4(I-N) REAL*8 LAT,LONG,TIME REAL*8 DARSIN,DARCOS,DTAN REAL*8 A, R, Z, HA, DEC, LHA REAL*8 DEGRD,PI,TOHMS,ONE INTEGER*4 P5,P6,P7,P8,P10 INTEGER*4 DAY,MONTH,YEAR,g C C LAT - LATITUDE (IN DEGREES AND FRACTION OF DEGREES), P19 C LONG - LONGITUDE (IN DEGREES AND FRACTIONS OF DEGREES),P20 C DAY - DAY OF THE MONTH,P16 C MONTH - MONTH OF YEAR,P17 C YEAR - YEAR,P18 C TIME - TIME (HHMM.SS) LOCAL STANDARD LST OR GMT, P4 C P5 - TIME TYPE (0 - LST, 1 - GMT) C P6 - APPARENT NOON CALCULATIONS (1-YES, 0-NO) C P7 - PRINT (0-NONE, 1-DATE, 2-DATA, 3 - NA COMPAR) C P8 - APPROXIMATION (0-FULL, 1-NO PLANETS, 2-NO LUNAR) C C RETURNED VARIABLES C C A - AZIMUTH ANGLE(IN DEGREES AND FRACTION OF DEGREE) C R - RADIUS VECTOR C Z - ZENITH ANGLE(IN DEGREES AND FRACTION OF DEGREE) C LHA - LOCAL HOUR ANGLE (IN DEGREES) C TIME - APPARENT NOON IF P6 = 1 (LST OR GMT SEE P5) C C P21 - TIME (IN FRACTIONS OF DAY) C P51 - GMT IN DEGREES C P53 - ZONE LONGITUDE C C DATA ONE/1.0D0/ C C DTAN(P)=DSIN(P)/DCOS(P) DARSIN(P)=DATAN2(P,DSQRT(1.0D0-P**2)) DARCOS(P) = PID2 - DARSIN(P) C C C----------------------------------------------------------------------- g=1 C----------------------------------------------------------------------- C C PID2 = DATAN(1.0D0)*2.0D0 PI = acos(-1.d0) DEGRD = PI/180.D0 C C++++++++++++++++++++++++++++++++++++++++++++ C P57 - (ET-UT), DIFFERENCE OF EPHEMERIS & UNIVERSAL TIME P57 = 0 C C P16 = DAY P17 = MONTH P18 = YEAR P19 = LAT P20 = LONG P54 = 10 P55 = 14 P10 = P7 1 CONTINUE P53 = 15.*AINT((7.5+ABS(LONG))/15.0)*DSIGN(1.d0,LONG) P11 = P53 IF(P5.EQ.1) P11=0D0 IF(P6.NE.1) GO TO 10 P56 = (P54+P55)/2D0 P51 = P56/24D0*360D0+P53 g=g+1 if(g.ge.100) go to 620 GO TO 20 10 CONTINUE P21 = (DMOD(TIME,1.d0)*100/3600+DMOD(AINT(TIME)/100.d0,1.d0) * *100./60. + AINT(AINT(TIME)/100.0))/24.0 P51 = P21*360D0+P11 20 CONTINUE I = P18 J = P17+1 IF(J.GT.3) GO TO 40 J = J + 12 I = I - 1 40 CONTINUE P13 = INT(I*365.25) + INT(J*30.6) + P16 P22 = P13 - 694038.5D0 P11 = P18 + P17/100.D0 + P16/1.0D4 IF(P11.LT.1900.0228D0) P22 = P22 + 1 IF(P11.LT.1800.0228D0) P22 = P22+1 P23 = (P51/360D0+P22+P57)/36525D0 P22 = P23*36525D0 C C MEAN LONGITUDE - P24 C P11 = 279.69668D0 P12 = 0.9856473354D0 P13 = 3.03D-4 P24 = P11 + DMOD(P12*P22,360D0) + P13*P23**2 P24 = DMOD(P24,360D0) C C MEAN ANOMALY - P25 C P11 = 358.47583D0 P12 = 0.985600267D0 P13 = -1.5D-4 P14 = -3.D-6 P25 = P11 + DMOD(P12*P22,360D0) + P13*P23**2 + P14*P23**3 P25 = P11 + DMOD(P12*P22,360D0) + P13*P23**2 + P14*P23**3 P25 = DMOD(P25,360D0) C C ECCENTRICITY - P26 C P11 = 0.01675104D0 P12 = -4.18D-5 P13 = -1.26D-7 P26 = P11 + P12*P23 + P13*P23**2 P11 = P25*DEGRD P12 = P11 C C ECCENTRIC ANOMALY - P13 (TEMP) C 100 CONTINUE P13 = P12 P12 = P11 + P26*DSIN(P13) IF(DABS((P12-P13)/P12).GT.1.0D-8) GO TO 100 P13 = P12/DEGRD C C TRUE ANOMALY - P27 C P27 =2.0D0*DATAN(DSQRT((1.0D0+P26)/(1.0D0-P26))*DTAN(P13/2.0D0 * *DEGRD))/DEGRD IF(DSIGN(1.0D0,P27).NE.DSIGN(1.0D0,DSIN(P13*DEGRD))) P27 = P27 + 1 IF(DSIGN(1.0D0,P27).NE.DSIGN(1.0D0,DSIN(P13*DEGRD))) * P27 = P27 + 180.0D0 IF(P27.LT. 0.0D0) P27 = P27 + 360D0 C C RADIUS VECTOR - R C R = 1.0D0 - P26*DCOS(P13*DEGRD) C C ABERRATION - P29 C P29 = -20.47/R/3600.0 C C MEAN OBLIQUITY - P43 C P11 = 23.452294D0 P12 = -0.0130125D0 P13 = -1.64D-6 P14 = 5.03D-7 P43 = P11 + P12*P23 + P13*P23**2 + P14*P23**3 C C MEAN ASCENSION - P45 C P11 = 279.6909832D0 P12 = 0.98564734D0 P13 = 3.8708D-4 P45 = P11 + DMOD(P12*P22,360D0) + P13*P23**2 P45 = DMOD(P45,360D0) C C----------------------------------------------------------------------- C CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA C C NUTATION AND LONGITUDE PERT C IF(P8.GT.1) GO TO 300 C----------------------------------------------------------------------- C C MOON'S MEAN ANOMALY - P28 P11 = 296.104608D0 P12 = 1325D0*360D0 P13 = 198.8491083D0 P14 = .00919167D0 P15 = 1.4388D-5 P28 = P11 + DMOD(P12*P23,360D0) + P13*P23 + P14*P23**2 + * P15*P23**3 P28 = DMOD(P28,360D0) C C MEAN ELONGATION OF MOON - P30: C P11 = 350.737486 P12 = 1236D0 * 360D0 P13 = 307.1142167D0 P14 = -1.436D-3 P30 = P11 + DMOD(P12*P23,360D0) + P13*P23 + P14*P23**2 P30 = DMOD(P30,360D0) C C MOON LONG OF ASCENDING NODE - P31: C P11 = 259.183275D0 P12 = -5D0 * 360D0 P13 = -134.142008D0 P14 = 2.0778D-3 P31 = P11 + DMOD(P12*P23,360D0) + P13*P23 + P14*P23**2 P31 = DMOD(P31,360D0) C C MEAN LONG OF MOON - P32: C P11 = 270.434164D0 P12 = 1336D0 * 360D0 P13 = 307.8831417D0 P14 = -1.1333D-3 P32 = P11 + DMOD(P12*P23,360D0) + P13*P23 + P14*P23**2 P32 = DMOD(P32,360D0) C C MOON PERTURBATION OF SUN LONG - P74: C P74 = 6.454D0*DSIN(P30*DEGRD) + .013D0*DSIN(3D0*P30*DEGRD) + * .177D0*DSIN((P30+P28)*DEGRD) - .424D0*DSIN((P30-P28)*DEGRD) P74 = P74 + 0.039D0*DSIN((3D0*P30-P28)*DEGRD) * - 0.064D0*DSIN((P30+P25)*DEGRD) + 0.172D0*DSIN((P30-P25)*DEGRD) P74 = P74/3600D0 P33 = P74 C C NUTATION OF LONG - P34 C P34 = -(17.2327D0-.017D0*P23) * DSIN(P31*DEGRD) + * .2088D0*DSIN(2D0*P31*DEGRD) - .2037D0*DSIN(2D0*P32*DEGRD) P34 = P34 - 1.2729D0*DSIN(2D0*P24*DEGRD) + .1261D0*DSIN(P28*DEGRD) P34 = P34/3600D0 C C NUTATION IN OBLIQUITY - P34: C P35 = 9.210D0*DCOS(P31*DEGRD) + .5522D0*DCOS(2D0*P24*DEGRD) - * .09D0*DCOS(2D0*P31*DEGRD) + .0889D0*DCOS(2D0*P32*DEGRD) P35 = P35/3600D0 C C INEQUALITIES OF LONG PERIOD - P36: C P36 = .266D0*DSIN((31.8D0+119D0*P23)*DEGRD) + *((1.882D0-.016D0*P23)*DEGRD) * DSIN((57.24D0+150.27D0*P23)*DEGRD) P36 = P36 + .202D0*DSIN((315.6D0+893.3D0*P23)*DEGRD) + * 1.089D0*P23**2 + 6.4D0*DSIN((231.19D0+20.2D0*P23)*DEGRD) P36 = P36/3600D0 C C MOON MEAN ARGUMENT OF LATITUDE - P63: C P11 = 11.250889D0 P12 = 1342D0*360D0 P13 = 82.02515D0 P14 = .003211D0 P63 = P11 + DMOD(P12*P23,360D0) + P13*P23 + P14*P23**2 C C----------------------------------------------------------------------- C CBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB C C PERTURBATIONS DUE TO PLANETS - P33: C IF(P8.GT.0) GO TO 300 C C VENUS: C C MEAN ANOMALY OF VENUS - P37: C P11 = 212.603222D0 P12 = 162D0*360D0 P13 = 197.803875D0 P14 = 1.286D-3 P37 = P11 + DMOD(P12*P23,360D0) + P13*P23 + P14*P23**2 P37 = DMOD(P37,360D0) P11 = 4.838D0*DCOS((299.171D0+P37-P25)*DEGRD) + * 5.526D0*DCOS((148.5259D0+2D0*P37-2D0*P25)*DEGRD) P11 = P11 + 2.497D0*DCOS((316.5759D0+2D0*P37-3D0*P25)*DEGRD) + * .666D0*DCOS((177.71D0+3D0*P37-3D0*P25)*DEGRD) P11 = P11 + 1.559D0*DCOS((345.4559D0+3D0*P37-4D0*P25)*DEGRD) + * 1.024D0*DCOS((318.15D0+3D0*P37-5D0*P25)*DEGRD) P70 = P11 P33 = P33+P11/3600D0 C C PERT. OF LATITUDE OF SUN BY VENUS - P61: C P61 = .21D0*DCOS((151.8D0+3D0*P37-4D0*P25)*DEGRD) + * .092D0*DCOS((93.7D0+P37-2D0*P25)*DEGRD) + * .067D0*DCOS((123D0+2D0*P37-3D0*P25)*DEGRD) C C MARS: C C MEAN ANOMALY OF MARS - P38: C P11 = 319.529022D0 P12 = 53D0*360D0 P13 = 59.8592194D0 P14 = 1.8083D-4 P38 = P11 + DMOD(P12*P23,360D0) + P13*P23 + P14*P23**2 P38 = DMOD(P38,360D0) P11 = .273D0*DCOS((217.7D0-P38+P25)*DEGRD) + * 2.043D0*DCOS((343.888D0-2D0*P38+2D0*P25)*DEGRD) P11 = P11 + 1.77D0*DCOS((200.4017D0-2D0*P38+P25)*DEGRD)+ * .425D0*DCOS((338.88D0-3D0*P38+2D0*P25)*DEGRD) P11 = P11 + .5D0*DCOS((105.18D0-4D0*P38+3D0*P25)*DEGRD) + * .585D0*DCOS((334.06D0-4D0*P38+2D0*P25)*DEGRD) P71=P11 P33 = P33+P11/3600D0 C C JUPITER: C C MEAN ANOMALY OF JUPITER - P39: C P11 = 225.3225D0 P12 = 8D0*360D0 P13 = 154.583D0 P39 = P11 + DMOD(P12*P23,360D0) + P13*P23 P39 = DMOD(P39,360D0) P11 = 7.208D0*DCOS((179.5317D0-P39+P25)*DEGRD) + * 2.6D0*DCOS((263.2167D0-P39)*DEGRD) P11 = P11 + 2.731*DCOS((87.1450D0-2D0*P39+2D0*P25)*DEGRD) P11 = P11 + 1.61D0*DCOS((109.4933D0-2D0*P39+P25)*DEGRD) + * .556D0*DCOS((82.65D0-3D0*P39+2D0*P25)*DEGRD) P72=P11 P33 = P33+P11/3600D0 C C PERT OF SOLAR LATITUDE BY JUPITER - P62: C P62 = .166D0*DCOS((265.5D0-2D0*P39+P25)*DEGRD) C C SATURN: C C MEAN ANOMALY OF SATURN - P40: C P11 = 175.613D0 P12 = 3D0*360D0 P13 = 141.794D0 P40 = P11 + DMOD(P12*P23,360D0) + P13*P23 P40 = DMOD(P40,360D0) P11 = .419*DCOS((100.58D0-P40+P25)*DEGRD) + * .32D0*DCOS((269.46D0-P40)*DEGRD) P73=P11 P33 = P33+P11/3600D0 C C----------------------------------------------------------------------- C CBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB C CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA C C PARAMETERS C 300 CONTINUE C C PRECESSION: C P42 = (50.2564D0+.0222D0*(P18-1900D0)/100D0) * * (P23-(P18-1900D0)/100D0) * 100D0/3600D0 C C APPARENT LONGITUDE - P41: C P41 = P27-P25+P24+P29+P33+P36+P34 C C SOLAR LATITUDE - P60: C P60 = (.576D0*DSIN(P63*DEGRD) + P61+P62)/3600D0 C C OBLIQUITY - P75 C P75 = P35 + P43 C C APPARENT RIGHT ASCENSION - P44: C P44 = DATAN(DTAN(P41*DEGRD)*DCOS(P75*DEGRD))/DEGRD IF(DSIGN(ONE,P44).NE.DSIGN(ONE,DSIN(P41*DEGRD))) P44 = P44+180D0 IF(P44.LT.0.D0) P44 = 360D0+P44 C C EQUATION OF TIME - P46: C P46 = P45-P44 IF(P46.GT.180D0) P46 = P46-360D0 C C HOUR ANGLE - P48: C P48 = P51+P46-180D0 HA = P48 C C LOCAL HOUR ANGLE - P49 C P49 = P48-P20 P49 = DMOD(P49,360D0) IF(P49.LT.0) P49 = P49 + 360D0 LHA=P49 C C C C DECLINATION - P47: C P47 = DARSIN(DSIN(P41*DEGRD)*DSIN(P75*DEGRD)*DCOS(P60 *DEGRD) + * DSIN(P60*DEGRD) * DCOS(P75*DEGRD))/DEGRD DEC = P47 C C ZENITH ANGLE - Z: C Z = DARCOS(DSIN(P19*DEGRD)*DSIN(P47*DEGRD) + DCOS(P19*DEGRD) * * DCOS(P47*DEGRD) * DCOS(P49*DEGRD))/DEGRD C C AZIMUTH ANGLE - A: C P11 = (-DSIN(P19*DEGRD) * DCOS(P49*DEGRD) * DCOS(P47*DEGRD) + * DSIN(P47*DEGRD) *DCOS(P19*DEGRD))/DSIN(Z*DEGRD) IF(DABS(P11).GT.1.0D0) P11 = DSIGN(1.0D0,P11) A= DARCOS(P11)/DEGRD IF(DSIGN(1.D0,-DCOS(P47*DEGRD) * DSIN(P49*DEGRD)/DSIN(Z*DEGRD)) * .NE.DSIGN(1D0,DSIN(A*DEGRD))) A = 360-A C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF(P10.EQ.0) GO TO 30 CALL NUPAGE PRINT 2000,MONTH,DAY,YEAR 2000 FORMAT(' MONTH',I5,/,' DAY',I5,/,' YEAR',I6) IF(P8.EQ.0) PRINT 2001 IF(P8.EQ.1) PRINT 2002 IF(P8.EQ.2) PRINT 2003 2001 FORMAT(' FULL COMPUTATIONS') 2002 FORMAT(' NO PLANETS') 2003 FORMAT(' NO MOON') PRINT 2004,LAT,LONG 2004 FORMAT(' LATITUDE ', F10.4,/,' LONGITUDE',F10.4) IF(P6.NE.1.AND.P5.EQ.0) PRINT 2050,TIME 2050 FORMAT(' TIME (LST)',F10.2) IF(P6.NE.1.AND.P5.EQ.1) PRINT 2051,TIME 2051 FORMAT(' TIME (GMT)',F10.2) C IF(P10.LT.2) GO TO 30 C CDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD C PRINT 2005,P23,P22 2005 FORMAT(' T',F15.6,' D',F15.6) P0 = TOHMS(P24) PRINT 2006,P0 2006 FORMAT(' MEAN LONG',F15.5) P0 = TOHMS(P25) PRINT 2007,P0 2007 FORMAT(' MEAN ANOMALY',F15.5) PRINT 2008,P26 2008 FORMAT(' ECCENTRICITY',F15.5) P0 = TOHMS(P27-P25) PRINT 2009,P0 2009 FORMAT(' EQUATION OF CENTER',F15.5) PRINT 2010,R 2010 FORMAT(' RADIUS VECTOR',F15.5) P0 = TOHMS(P29) PRINT 2011,P0 2011 FORMAT (' ABERRATION ',F15.5) P0 = TOHMS(P43) PRINT 2012,P0 2012 FORMAT(' MEAN OB',F15.5) P0 = TOHMS(P45/360D0*24D0) PRINT 2013,P0 2013 FORMAT(' MEAN ASCENSION',F15.5) IF(P10.GT.1) GO TO 35 P0 = TOHMS(P74) PRINT 2014, P0 2014 FORMAT (' P MOON',F15.5) IF(P10.GT.0) GO TO 35 P0 = P70/3600D0 P0 = TOHMS(P0) PRINT 2015, P0 2015 FORMAT (' P VENUS',F15.5) P0 = P61/3600D0 P0 = TOHMS(P0) PRINT 2016, P0 2016 FORMAT (' P LAT BY VENUS ',F15.5) P0 = P71/3600D0 P0 = TOHMS(P0) PRINT 2017, P0 2017 FORMAT (' P MARS ',F15.5) P0 = P72/3600D0 P0 = TOHMS(P0) PRINT 2018, P0 2018 FORMAT (' P JUPITER ',F15.5) P0 = P62/3600D0 P0 = TOHMS(P0) PRINT 2019, P0 2019 FORMAT (' P OF LAT BY JUP ',F15.5) P0 = P73/3600D0 P0 = TOHMS(P0) PRINT 2020, P0 2020 FORMAT (' P SATURN ',F15.5) C 35 CONTINUE C P0 = TOHMS(P33) PRINT 2021, P0 2021 FORMAT (' PERTURBATIONS ',F15.5) P0 = TOHMS(P36) PRINT 2022, P0 2022 FORMAT (' LONG PERIOD ', F15.5) P0 = TOHMS(P34) PRINT 2023, P0 2023 FORMAT (' NUTATION OF LONG ', F15.5) P0 = TOHMS(P35) PRINT 2024, P0 2024 FORMAT (' NUT OBLIQUITY ', F15.5) P0 = TOHMS(P42) PRINT 2025, P0 2025 FORMAT (' PRECESSION ', F15.5) P0 = TOHMS(P41) PRINT 2026, P0 2026 FORMAT (' APPAR. LONG. ',F15.5) P0 = TOHMS(P60) PRINT 2027, P0 2027 FORMAT (' SOLAR LATITUDE ',F15.5) P0 = TOHMS(P75) PRINT 2028, P0 2028 FORMAT (' OBLIQUITY ',F15.5) P0 = P44/360D0*24D0 P0 = TOHMS(P0) PRINT 2029, P0 2029 FORMAT (' APPAR. ASCENSION ',F15.5) P0 = P46/360D0*24D0 P0 = TOHMS(P0) PRINT 2030, P0 2030 FORMAT (' EQUATION OF TIME ', F15.5) P0 = TOHMS(P49) PRINT 2032, P0 2032 FORMAT (' LHA = ',F15.5) P0 = TOHMS(P47) PRINT 2033, P0 2033 FORMAT (' DECLIN. ',F15.5) C CDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C 30 CONTINUE C CEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE C C IF(P6.EQ.0) GO TO 500 P12 = P53/360D0*24D0 IF(P5.EQ.0) P12=0D0 P10 = 0 IF(DMOD(A,180.d0).GT.0.01) GO TO 400 P11 = P56+P12 IF(P11.LT.0D0) P11 = P11+24D0 P13 = DMOD(P11,1D0)*60.0D0 TIME = DINT(P11)*100.0D0+DINT(P13) * + DMOD(P13,1.D0)*60.0D0/100.D0 IF(P56+P12.LT.0D0) TIME = -TIME RETURN C C 400 CONTINUE IF(A.LT.90.0.OR.A.GT.270.0) GO TO 410 IF(A.LE.180.) GO TO 405 402 CONTINUE P55 = P56 GO TO 1 405 CONTINUE P54 = P56 GO TO 1 410 CONTINUE IF(A.LT.90.0) GO TO 405 GO TO 402 C C 500 CONTINUE C CEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE C IF(P10.EQ.0) GO TO 550 P0 = A P0 = TOHMS(P0) PRINT 2034, P0 2034 FORMAT (' A DDD.MMSSS ',F15.5) P0 = Z P0 = TOHMS(P0) PRINT 2035, P0 2035 FORMAT (' Z DD.MMSSS ',F15.5) C 550 CONTINUE IF(P7.LT.3) RETURN C CFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF C IF(P10.EQ.0) GO TO 600 P52 = P46 P0 = P41-(P29+P34) P0 = TOHMS(P0) PRINT 2036, P0 2036 FORMAT (' LONGITUDE ',F15.6) P0 = (P29+P34)*10000 P0 = TOHMS(P0) PRINT 2037, P0 2037 FORMAT (' REDN TO LONG ',F15.6) P0 = P60*10000 P0 = TOHMS(P0) PRINT 2038, P0 2038 FORMAT (' LATITUDE ',F15.6) P0 = P42*10000 P0 = TOHMS(P0) PRINT 2039, P0 2039 FORMAT (' PRECESSION ',F15.6) P0 = P34*10000 P0 = TOHMS(P0) PRINT 2040, P0 2040 FORMAT (' NUT IN LONG ',F15.6) P0 = P35*10000 P0 = TOHMS(P0) PRINT 2041, P0 2041 FORMAT (' NUT IN OBLI ',F15.6) P0 = TOHMS(P43) PRINT 2042, P0 2042 FORMAT (' OBLIQUITY ',F15.6) P0 = P44/360D0*24D0 P0 = TOHMS(P0) PRINT 2043, P0 2043 FORMAT (' APP RT. ASCEN ',F15.6) P0 = TOHMS(P47) PRINT 2044, P0 2044 FORMAT (' APP DECLINTION ',F15.6) C 600 CONTINUE P16 = P16+1.D0 IF(P10.EQ.0) GO TO 610 P10 = 0D0 GO TO 1 610 CONTINUE P11 = 180.D0-P52-(P46-P52)*((-P52+180.D0)/360D0) P0 = TOHMS(P11/360.D0*24D0) PRINT 2045,P0 2045 FORMAT (' ET ',F15.6) C C 620 continue RETURN C CFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF C C *** END SUBROUTINE EPHEMS *** END C C SUBROUTINE NUPAGE C DUMMY NEW PAGE ROUTINE PRINT 2000 2000 FORMAT(1H1) RETURN END C C REAL*8 FUNCTION TOHMS(P1) C FUNCTION RETURNS P IN HH MM SS FORMAT IMPLICIT REAL*8(A-H,O-Z) DATA ONE/1.0D0/ P3 = DABS(P1) B = DMOD(P3,ONE)*60. C = DINT(B) D = C/100. P2 = DINT(P3) + D + DMOD(B,ONE)*60./10000. P2 = P2*DSIGN(ONE,P1) TOHMS = P2 RETURN END