c Computes statistics for comparison of precip. estimates from various c algorithms with CPC real-time daily rain gauge analysis c c NOTE: look at output_statistics/INFO for info about order of scores & c statistics in the data sets c cccc 01/21/2005: added "scampr" algo. c cccc 04/01/2005: Changes made to Main & "rdvald" and to "scores" and c "stats" subroutines c c************************************************************************** c PARAMETER (nx=241, ny=121, xmsng=-9999., nalgo=24, x ix1=1, ix2=241, jy1=1,jy2=121, iday1=92, iyr1=2003) c real*4 algo(nx,ny,nalgo), gag(nx,ny), xmask(nx,ny), x ctabl(2,2,nalgo), a(nx*ny), b(nx*ny), gagn(nx,ny), x rad(nx,ny), radn(nx,ny) c integer*4 iexist(nalgo), istag(nalgo) logical there c character*26 infile character*27 infile2 character*72 out_ascii character*72 out_grads character*8 jdate character*6 calgo (nalgo) character*8 title(nalgo) c data infile /'data_regrid/ '/ data infile2 /'data_regrid/ '/ data io_algo /1/ data io_date /9/ data io_out /20/ data iout_gr /50/ data jxistg /0/, jxistr /0/ c c------------------------------------------------------------------------------------------ c data calgo /'nrlgeo', 'nrlbld', 'gsfc40', 'gsfc41', x 'gsfc42' ,'cstxxx', x 'cmorph', 'gfsprc', 'radar', 'mwcomb', 'pramsu', x 'gpiprc', 'gmsrax', 'irfreq', 'hydroe', 'nogaps', x 'gmsrad', 'emorph', 'imorph', 'persin', 'ukpmir', x 'gmorph','scampr', 'qprop1'/ data title /'NRL-GEO ', 'NRL-PMW ', '3B40RT ', '3B41RT ', x '3B42RT ', 'CST ', x 'CMORPH ', 'NCEP_MDL', 'RADAR ','CPC-PMW ', 'AMSU-B ', x 'GPI ', 'GMSRA ', 'IR_FREQ ', 'HYDRO-E ', 'NOGAPS ', x 'GMSRAD ', 'CMORPHEV', 'CMORPHIR', 'PERSIANN', 'UKPMIR ', x 'GMORPH ', 'SCAMPR ', 'QMORPH '/ c c "If "istag=1" then the estimate is on the "staggered grid" (i.e. lattice-centered) c data istag / 1, 1, 1, 1, x 1, 1, x 0, 0, 0, 0, 0, x 0, 0, 0, 0, 0, x 0, 0, 0, 0, 0, x 0, 0, 0/ c c------------------------------------------------------------------------------------------ c c..... Read in land-sea mask c open (89,file='data_regrid/mask.data',access='direct',recl=nx*ny) read(89,rec=1) xmask c c call gtdate(io_date,iyr1,iday1,jdate,ndays) call rdvald(infile,infile2,nx,ny,jdate,gag,gagn,rad,radn, x jxistg,jxistr) call init(algo,iexist,nx,ny,iy1,iy2,jy1,jy2,nalgo) c do ialgo=1,nalgo call rdalgo(io_algo,infile2,algo,iexist,nx,ny,ialgo,calgo,nalgo, x jdate) enddo c call chkdat(algo,nx,ny,nalgo,iexist) out_ascii = 'output_statistics/----- ' out_grads = 'output_statistics_grads/----- ' call mask(xmask,nx,ny,xmsng,gag,gagn) ! Mask out gage data over ocean call fair(xmsng,nx,ny,nalgo,algo,gag,gagn,iexist) c if(jxistg.eq.1) then call vgauge(xmsng,nx,ny,ix1,ix2,jy1,jy2,nalgo,algo,gag,gagn,istag, x ctabl,calgo,jdate,iexist,out_ascii,out_grads,iout, x iout_gr,a,b,ndays,title) endif out_ascii = 'output_statistics/------/------_vs_radar.yyyymmdd' out_grads = 'output_statistics_grads/radar.----- ' call mask(xmask,nx,ny,xmsng,rad,radn) ! Mask out gage data over ocean call fair(xmsng,nx,ny,nalgo,algo,rad,radn,iexist) c if(jxistr.eq.1) then call vradar(xmsng,nx,ny,ix1,ix2,jy1,jy2,nalgo,algo,rad,radn,istag, x ctabl,calgo,jdate,iexist,out_ascii,out_grads,iout, x iout_gr,a,b,ndays,title) endif stop end c c............................................................................. c subroutine chkdat(algo,nx,ny,nalgo,iexist) real*4 algo(nx,ny,nalgo) integer*4 iexist(nalgo) c c computes the average of each algorithm, and if teh average is < -9000 then the c element of "iexist" that corresponds tho that algorithm is set to "0". This c is done because some fields may have a file for a particular data, but all the c data may be missing, then the "fair" subroutine sets ALL valaues to missing for c ALL algorithms which is not desirable. c c do ialg=1,nalgo s=0. xn=0. write(*,*) 'ialg=',ialg do j=1,ny do i=1,nx s=s+algo(i,j,ialg) xn=xn+1. enddo enddo avg=s/xn if(avg.lt.-9000.) iexist(ialg)=0 write(*,*) 'ialg,avg,iexist(ialg)',ialg,avg,iexist(ialg) enddo return end c c............................................................................. subroutine chkgag(gag,nx,ny,iyesno) c real*4 gag(nx,ny) c c computes the average of gauge data & sets "iyesno" to "0" if avg < 0. c c iyesno=1 s=0. xn=0. do j=1,ny do i=1,nx s=s+gag(i,j) xn=xn+1. enddo enddo avg=s/xn if(avg.lt.-9000.) then iyesno=0 write(*,*) 'Gauge avg. precip. = ',avg endif return end c............................................................................. c subroutine contin(xmsng,nx,ny,ix1,ix2,jy1,jy2,nalgo,algo,gag, x gagn,istag,ctabl) real*4 algo(nx,ny,nalgo), gag(nx,ny), ctabl(2,2,nalgo), x gagn(nx,ny) integer*4 istag(nalgo) c c Assembles a contingency table for rain/no-rain (< 1mm/day vs. >= 1mm/day) c for each technique to be validated in array "ctabl". The "i" & "j" c elements represent the rows and columns of the contingency table. c NOTE that prior to entering this subroutine, a "call" should have been c made in "Main" to subroutine "fair" which, when a value is missing at an c (i,j), set the value to missing at that (i,j) in ALL of the algorithms c (i.e. algo(i,j,-). do k=1,nalgo do j=1,2 do i=1,2 ctabl(i,j,k)=0. enddo enddo enddo c c do j=jy1,jy2 do 100 i=ix1,ix2 if(gagn(i,j).eq.xmsng .or. gag(i,j).eq.xmsng) go to 100 commented out 3/28/2005 if(algo(i,j,k).eq.xmsng) go to 100 c do k=1,nalgo if(istag(k).eq.1) then c if(gagn(i,j).ge.1. .and. algo(i,j,k).ge. 1.) x ctabl(1,1,k)=ctabl(1,1,k)+1. c if(gagn(i,j).lt. 1. .and. algo(i,j,k).ge. 1.) x ctabl(2,1,k)=ctabl(2,1,k)+1. c if(gagn(i,j).ge. 1. .and. algo(i,j,k).lt. 1.) x ctabl(1,2,k)=ctabl(1,2,k)+1. c if(gagn(i,j).lt. 1. .and. algo(i,j,k).lt. 1.) x ctabl(2,2,k)=ctabl(2,2,k)+1. c else c if(gag(i,j).ge.1. .and. algo(i,j,k).ge. 1.) x ctabl(1,1,k)=ctabl(1,1,k)+1. c if(gag(i,j).lt. 1. .and. algo(i,j,k).ge. 1.) x ctabl(2,1,k)=ctabl(2,1,k)+1. c if(gag(i,j).ge. 1. .and. algo(i,j,k).lt. 1.) x ctabl(1,2,k)=ctabl(1,2,k)+1. c if(gag(i,j).lt. 1. .and. algo(i,j,k).lt. 1.) x ctabl(2,2,k)=ctabl(2,2,k)+1. endif enddo 100 continue enddo return end c c............................................................................. c SUBROUTINE CORR(A,B,NUM,R,SLOPE,YINT,IPRINT) C C JANOWIAK - OCTOBER, 1988 C C C************************************************************* C * C CORRELATION FORMULATION FROM PANOFSKY & BRIER * C SLOPE & INTERCEPT FORMULATION FROM SHAUM'S OUTLINE * C * C************************************************************* C C DIMENSION A(NUM), B(NUM) C REAL*8 ABAR,BBAR,ABBAR,ASQ,BSQ,ASIG,BSIG,ASUM,BSUM C IF(IPRINT.NE.1) GO TO 150 PRINT 105,(A(I),I=1,NUM) PRINT 115,(B(I),I=1,NUM) 105 FORMAT(///' A:',(10F8.2)) 115 FORMAT(///' B:',(10F8.2)) C 150 ASUM=0. BSUM=0. ABBAR=0. ASQ=0. BSQ=0. ASQSUM=0. BSQSUM=0. ABSUM=0. C DO 300 I=1,NUM ASUM=ASUM+A(I) BSUM=BSUM+B(I) ABSUM=ABSUM+A(I)*B(I) ASQSUM=ASQSUM+A(I)*A(I) BSQSUM=BSQSUM+B(I)*B(I) 300 CONTINUE C ABAR=ASUM/NUM BBAR=BSUM/NUM ABBAR=ABSUM/NUM ASQ=ASQSUM/NUM BSQ=BSQSUM/NUM ASIG=SQRT(ASQ-ABAR*ABAR) BSIG=SQRT(BSQ-BBAR*BBAR) C IF(ASIG.GT.0. .AND. BSIG.GT.0.) GO TO 310 R=-999. SLOPE=-999. CCC WRITE(*,305) 305 FORMAT(' FROM "CORR": SIGMAS OF AT LEAST 1 TIME SERIES', 1 ' ARE ZERO - SETTING "R" & "SLOPE" TO "-9999."') GO TO 350 310 R=(ABBAR-ABAR*BBAR)/(ASIG*BSIG) SLOPE=(NUM*ABSUM-ASUM*BSUM)/(NUM*ASQSUM-ASUM*ASUM) YINT = (BSUM*ASQSUM-ASUM*ABSUM)/(NUM*ASQSUM-ASUM*ASUM) C 350 RETURN END 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 c c--------------------------------------------------------------------- c 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 subroutine fair(xmsng,nx,ny,nalgo,algo,gag,gagn,iexist) real*4 algo(nx,ny,nalgo), gag(nx,ny), gagn(nx,ny) integer*4 iexist(nalgo) c c Sets value to missing in all arrays wherever data are missing in c any array. This is done to ensure that verification is uniform among c all estimates to be verified. c do j=1,ny do 100 i=1,nx if(gag(i,j).eq.xmsng .or. gagn(i,j).eq.xmsng) then do k=1,nalgo algo(i,j,k)=xmsng enddo go to 100 endif c do 50 n=1,nalgo if(iexist(n).eq.0) go to 50 cjul23,20: do not set values to missing in other arrays based on "3b40rt" data if(n.eq.3) go to 50 if(algo(i,j,n).eq.xmsng) then do k=1,nalgo gag(i,j)=xmsng gagn(i,j)=xmsng algo(i,j,k)=xmsng enddo go to 100 endif 50 continue 100 continue enddo open (77,file='test',access='direct',recl=nx*ny) write(77,rec=1) gag return end c c............................................................................. c subroutine gtdate(io_date,iyr1,iday1,jdate,ndays) integer*4 ndys(2003:2020) character*8 jdate character*4 lyr character*2 lmn,ldy c c ..... following is the # of days in each year during 2003 thru 2020 data ndys /365,366,365,365,365,366,365,365,365,366,365,365,365, x 366,365,365,365,366/ c open (io_date,file='date_file') read(io_date,'(a4)') lyr read(io_date,'(a2)') lmn read(io_date,'(a2)') ldy jdate=lyr//lmn//ldy c write(*,*) 'lyr,lmn,ldy,jdate:',lyr,lmn,ldy,jdate c c c---------------------------------------------------------- c c "iday1" is used to compute the direct access record to which c statistics are to be written in a data set for time series c plots. It is the day-of-the-year (DOY) in 2003 that the statistics c archive began (2 apr 2003). To compute the record number, c need to get date of data being processed, then compute the c DOY for that date & determine the # of days since "iday1" c read(lyr,'(i4)') iyr read(lmn,'(i2)') imon read(ldy,'(i2)') idy call date2doy(iyr,imon,idy,idoy) if(iyr.eq.iyr1) then ndays=(idoy-iday1) + 1 else nd=(iyr-iyr1)*ndys(iyr-1) + idoy ndays=(nd-iday1)+1 endif return end c c............................................................................. c subroutine init(algo,iexist,nx,ny,iy1,iy2,jy1,jy2,nalgo) real*4 algo(nx,ny,nalgo) integer*4 iexist(nalgo) do n=1,nalgo iexist(n)=0 enddo c do k=1,nalgo do j=jy1,jy2 do i=iy1,iy2 algo(i,j,k)=-9999. enddo enddo enddo return end c c............................................................................. c subroutine mask(xmask,nx,ny,xmsng,gag,gagn) real*4 xmask(nx,ny), gag(nx,ny),gagn(nx,ny) c c Uses a land-sea mask to set rain gauge values in "gag" array to missing c over coastal regions. "xmask"=+1 means land; "xmask"=-1 means ocean c do j=1,ny do i=1,nx if(xmask(i,j).lt.0.) gag(i,j)=xmsng if(xmask(i,j).lt.0.) gagn(i,j)=xmsng enddo enddo return end c c............................................................................. c subroutine precor(algo,gag,gagn,istag,nx,ny,ix1,ix2,jy1,jy2,nalgo, x ialgo,xmsng,a,b,npts) real*4 algo(nx,ny,nalgo), gag(nx,ny), a(nx*ny), b(nx*ny ), x gagn(nx,ny) integer*4 istag(nalgo) npts=0 do j=jy1,jy2 do 100 i=ix1,ix2 if(gag(i,j).eq.xmsng .or. algo(i,j,ialgo).eq.xmsng) go to 100 if(gagn(i,j).eq.xmsng) go to 100 ccc if(gag(i,j).lt.1.) go to 100 npts=npts+1 a(npts)=algo(i,j,ialgo) c c.......... chg following if test when algo added that uses "gaugez" for validation c if(istag(ialgo).eq.1) then b(npts)=gagn(i,j) else b(npts)=gag(i,j) endif 100 continue enddo return end c c............................................................................. c subroutine rdalgo(io_algo,infile2,algo,iexist,nx,ny,ialgo,calgo, x nalgo,jdate) character*27 infile2 character*8 jdate character*6 calgo(nalgo) integer*4 iexist(nalgo) logical*4 there real*4 algo(nx,ny,nalgo) c if(calgo(ialgo).eq.'radar') then infile2(13:17)=calgo(ialgo) infile2(18:18)='/' infile2(19:26)=jdate infile2(27:27)=' ' else infile2(13:18)=calgo(ialgo) infile2(19:19)='/' infile2(20:27)=jdate endif write(*,*) infile2 inquire(file=infile2,exist=there) if(there) then write(*,*) 'opening file:',infile2 open (io_algo,file=infile2,access='direct',recl=nx*ny) read(io_algo,rec=1) ((algo(i,j,ialgo),i=1,nx),j=1,ny) iexist(ialgo)=1 close (io_algo) endif return end c c............................................................................. c subroutine rdvald(infile,infile2,nx,ny,jdate,gag,gagn,rad,radn, x jxistg,jxistr) cc real*4 gag(nx,ny), gagn(nx,ny), rad(nx,ny), radn(nx,ny) logical there character*26 infile character*27 infile2 character*8 jdate do j=1,ny do i=1,nx gag(i,j)=-9999. gagn(i,j)=-9999. rad(i,j)=-9999. radn(i,j)=-9999. enddo enddo c c....... Read in gauge data on cell-centered grid c infile(13:18)='gauge/' infile(19:26)=jdate inquire(file=infile,exist=there) if(there) then jxistg=1 write(*,*) 'opening file:',infile open (1,file=infile,access='direct',recl=nx*ny) read(1,rec=1) gag close (1) call chkgag(gag,nx,ny,iyesno) if(iyesno.eq.0) then write(*,*) 'No Gauge data for: ',jdate write(*,*) 'Program aborting' stop endif else write(*,*) 'No Gauge file for: ',jdate write(*,*) 'Program aborting' stop endif c c....... Read in gauge data on lattice-centered grid c infile2(13:19)='gaugez/' infile2(20:27)=jdate inquire(file=infile2,exist=there) if(there) then write(*,*) 'opening file:',infile2 open (1,file=infile2,access='direct',recl=nx*ny) read(1,rec=1) gagn call chkgag(gagn,nx,ny,iyesno) if(iyesno.eq.0) then write(*,*) 'No Gauge (staggered) data for: ',jdate write(*,*) 'Program aborting' stop endif else write(*,*) 'No Gauge file (staggered) for: ',jdate write(*,*) 'Program aborting' stop endif c c....... Read in radar data on cell-centered grid c infile(13:18)='radar/' infile(19:26)=jdate inquire(file=infile,exist=there) if(there) then jxistr=1 write(*,*) 'opening file:',infile open (1,file=infile,access='direct',recl=nx*ny) read(1,rec=1) rad close (1) call chkgag(rad,nx,ny,iyesno) if(iyesno.eq.0) then write(*,*) 'No Radar data for: ',jdate ccc write(*,*) 'Program aborting' ccc stop endif else write(*,*) 'No Radar file for: ',jdate ccc write(*,*) 'Program aborting' ccc stop endif c c....... Read in gauge data on lattice-centered grid c infile2(13:19)='radarz/' infile2(20:27)=jdate inquire(file=infile2,exist=there) if(there) then write(*,*) 'opening file:',infile2 open (1,file=infile2,access='direct',recl=nx*ny) read(1,rec=1) radn call chkgag(radn,nx,ny,iyesno) if(iyesno.eq.0) then write(*,*) 'No Radar (staggered) data for: ',jdate ccc write(*,*) 'Program aborting' ccc stop endif else write(*,*) 'No Radar file (staggered) for: ',jdate ccc write(*,*) 'Program aborting' ccc stop endif c return end c c............................................................................. c subroutine scores(nalgo,ialgo,ctabl,calgo,iout,iout_gr,ndays, x iexist) real*4 ctabl(2,2,nalgo) character*6 calgo (nalgo) integer*4 iexist(nalgo) c c hits = Proportion Correct c threat = Threat Score c pod = Probability of Detection c far = False Alarm Ratio c bias = bias score c hss = Heidke Skill Score c xkss = Kuiper's Score c c--------------------------------------------- c if(iexist(ialgo).eq.0) then xn=-999. hits=-999. ethret=-999. pod=-999. far=-999. bias=-999. hss=-999. hks=-999. a=-999. b=-999. c=-999. d=-999. xnr1=-999. xnr2=-999. go to 550 endif a = ctabl(1,1,ialgo) b = ctabl(2,1,ialgo) c = ctabl(1,2,ialgo) d = ctabl(2,2,ialgo) xn = a+b+c+d hits = (a+d)/xn c c.....Compute "threat" c if(a+b+c .gt.0.) then threat = a / (a+b+c) else threat=-999. endif c c..... compute "equitable" threat score (Ebert) c hitran = ((a+c)*(a+b)) / xn if (a+b+c - hitran .gt.0.) then ethret = (a - hitran) / (a+b+c - hitran) else ethret=-999. endif c c.............................................. c if(a+c .gt.0.) then pod = a / (a+c) bias = (a+b) / (a+c) hks = (a*d-b*c) / ( (a+c)*(b+d) ) else pod=-999. bias=-999. hks=-999. endif if(a+b .gt.0) then far = b / (a+b) else far=-999. endif if( (a+c)*(c+d) + (a+b)*(b+d) .gt.0.) then hss = 2*(a*d-b*c) / ( (a+c)*(c+d) + (a+b)*(b+d) ) else hss=-999. endif write(*,'(''n :'',f8.0)') xn write(*,'(''Proportion Correct :'',f8.3)') hits write(*,'(''Equitable Threat Score :'',f8.3)') ethret write(*,'(''Prob. of Detection :'',f8.3)') pod write(*,'(''False Alarm Ratio :'',f8.3)') far write(*,'(''Bias Score :'',f8.3)') bias write(*,'(''Heidke Skill Score :'',f8.3)') hss write(*,'(''Hanssen-Kuiper Score :'',f8.3)') hks write(*,*) write(*,'('' VALID'',3x,a6)') x calgo(ialgo) xnr1=a+c xnr2=a+b write(*,'(''# gridpnts raining :'',2f8.0)') xnr1,xnr2 550 write(iout,'(f8.0)') xn write(iout,'(f8.3)') hits write(iout,'(f8.3)') ethret write(iout,'(f8.3)') pod write(iout,'(f8.3)') far write(iout,'(f8.3)') bias write(iout,'(f8.3)') hss write(iout,'(f8.3)') hks write(iout,'(2f8.0)') xnr1 write(iout,'(2f8.0)') xnr2 write(iout_gr,rec=(ndays-1)*10+1) ethret write(iout_gr,rec=(ndays-1)*10+2) pod write(iout_gr,rec=(ndays-1)*10+3) far write(iout_gr,rec=(ndays-1)*10+4) bias write(iout_gr,rec=(ndays-1)*10+5) hss write(iout_gr,rec=(ndays-1)*10+6) hks return end c c............................................................................. c subroutine stats(algo,gag,gagn,istag,nx,ny,ix1,ix2,jy1,jy2,nalgo, x ialgo,xmsng,xmae,rmse,cgag,calg,rmaxg,rmaxa, x algomn,gagemn) real*4 algo(nx,ny,nalgo), gag(nx,ny), gagn(nx,ny) integer*4 istag(nalgo) c c..... Compute MAE & RMSE ................................. c err=0. errsq=0. algomn=0. gagemn=0. xn=0. do j=jy1,jy2 do 100 i=ix1,ix2 if(istag(ialgo).eq.1) then if(gagn(i,j).eq.xmsng .or. algo(i,j,ialgo).eq.xmsng)go to 100 xn=xn+1. algomn=algomn+algo(i,j,ialgo) gagemn=gagemn+gagn(i,j) err=err + abs(algo(i,j,ialgo)-gagn(i,j)) errsq=errsq + (algo(i,j,ialgo)-gagn(i,j)) * x (algo(i,j,ialgo)-gagn(i,j)) else if(gag(i,j).eq.xmsng .or. algo(i,j,ialgo).eq.xmsng) go to 100 xn=xn+1. algomn=algomn+algo(i,j,ialgo) gagemn=gagemn+gag(i,j) err=err + abs(algo(i,j,ialgo)-gag(i,j)) errsq=errsq + (algo(i,j,ialgo)-gag(i,j)) * x (algo(i,j,ialgo)-gag(i,j)) endif 100 continue enddo xmae=err/xn rmse=sqrt(errsq/xn) algomn=algomn/xn gagemn=gagemn/xn c c..... Compute Conditional Rain Rate & determine max. rain value c rmaxg=0. rmaxa=0. cgag=0. calg=0. xng=0. xna=0. do j=jy1,jy2 do 200 i=ix1,ix2 c if(istag(ialgo).eq.1) then if(gagn(i,j).eq.xmsng) go to 200 if(gagn(i,j).ge.1.) then cgag=cgag+gagn(i,j) xng=xng+1. if(gagn(i,j).gt.rmaxg) rmaxg=gagn(i,j) endif c if(algo(i,j,ialgo).ge.1.) then calg=calg+algo(i,j,ialgo) xna=xna+1. if(algo(i,j,ialgo).gt.rmaxa) rmaxa=algo(i,j,ialgo) endif else if(gag(i,j).eq.xmsng) go to 200 if(gag(i,j).ge.1.) then cgag=cgag+gag(i,j) xng=xng+1. if(gag(i,j).gt.rmaxg) rmaxg=gag(i,j) endif c if(algo(i,j,ialgo).ge.1.) then calg=calg+algo(i,j,ialgo) xna=xna+1. if(algo(i,j,ialgo).gt.rmaxa) rmaxa=algo(i,j,ialgo) endif endif 200 continue enddo if(xng.gt.0.) then cgag=cgag/xng else cgag=-999. endif if(xna.gt.0.) then calg=calg/xna else calg=-999. endif c return end c c------------------------------------------------------------------------------------- c subroutine vgauge(xmsng,nx,ny,ix1,ix2,jy1,jy2,nalgo,algo,gag,gagn, x istag,ctabl,calgo,jdate,iexist,out_ascii,out_grads, x iout,iout_gr,a,b,ndays,title) c c------- Validates estimates using RAIN GAUGE data as the reference standard c real*4 algo(nx,ny,nalgo), gag(nx,ny), xmask(nx,ny), x ctabl(2,2,nalgo), a(nx*ny), b(nx*ny), gagn(nx,ny) c integer*4 iexist(nalgo), istag(nalgo) logical there c character*72 out_ascii character*72 out_grads character*8 jdate character*6 calgo (nalgo) character*8 title(nalgo) call contin(xmsng,nx,ny,ix1,ix2,jy1,jy2,nalgo,algo,gag,gagn,istag, x ctabl) do 500 ialgo=1,nalgo ccc if(iexist(ialgo).eq.0) go to 500 write(iout,*) calgo(ialgo) if(calgo(ialgo).eq.'radar') then out_ascii(19:23)=calgo(ialgo) out_ascii(24:24)='/' out_ascii(25:29)=calgo(ialgo) out_ascii(30:30)='.' out_ascii(31:38)=jdate out_ascii(39:40)=' ' out_grads(25:29)=calgo(ialgo) out_grads(30:35)='.grads' out_grads(36:36)=' ' else out_ascii(19:24)=calgo(ialgo) out_ascii(25:25)='/' out_ascii(26:31)=calgo(ialgo) out_ascii(32:32)='.' out_ascii(33:40)=jdate out_grads(25:30)=calgo(ialgo) out_grads(31:36)='.grads' endif open (iout,file=out_ascii) open (iout_gr,file=out_grads,access='direct',recl=1) write(iout,*) calgo(ialgo) if(iexist(ialgo).eq.0) go to 550 cc if(iexist(ialgo).eq.0) then cc do istat=1,10 cc write(iout_gr,rec=(ndays-1)*10+istat) xmsng cc enddo cc do istat=2,25 cc write(iout,'(f8.0)') xmsng cc enddo cc go to 500 cc endif c c write(*,*) write(*,*) '********************************************' write(*,*) write(*,*) write(*,*) ' GAUGE' write(*,*) write(*,*) ' >=1 <1' write(*,*) write(*,'(8x,''>=1'',2f8.0)') (ctabl(i,1,ialgo),i=1,2) write(*,*) calgo(ialgo) write(*,'(9x,''<1'',2f8.0)') (ctabl(i,2,ialgo),i=1,2) write(*,*) c call precor(algo,gag,gagn,istag,nx,ny,ix1,ix2,jy1,jy2,nalgo,ialgo, x xmsng,a,b,npts) call corr(a,b,npts,r,slope,yint,0) write(*,'(''Correlation :'',f8.3)') r c call stats(algo,gag,gagn,istag,nx,ny,ix1,ix2,jy1,jy2,nalgo,ialgo, x xmsng,xmae,rmse,cgag,calg, rmaxg,rmaxa,algomn,gagemn) write(*,'(''Mean Abs Error :'',f8.3)') xmae write(*,'(''RMSE :'',f8.3)') rmse write(*,'(''Normalized (gauge)RMSE :'',f8.3)') rmse/gagemn c call scores(nalgo,ialgo,ctabl,calgo,iout,iout_gr,ndays,iexist) write(*,'(''Conditional rain rate :'',2f8.1)') cgag,calg write(*,'(''Maximum rain value :'',2f8.1)') rmaxg,rmaxa write(*,'(''Mean rain :'',2f8.1)') gagemn,algomn c c.... Write statistics to file for Graphics c 550 rmsen=rmse/gagemn diff=algomn-gagemn if(iexist(ialgo).eq.0) then ctabl(1,1,ialgo)=-999. ctabl(2,1,ialgo)=-999. ctabl(1,2,ialgo)=-999. ctabl(2,2,ialgo)=-999. r=-999. xmae=-999. rmse=-999. rmsen=-999. cgag=-999. calg=-999. rmaxg=-999. rmaxa=-999. gagemn=-999. algomn=-999. diff=-999. endif write(iout,'(f8.0)') ctabl(1,1,ialgo) write(iout,'(f8.0)') ctabl(2,1,ialgo) write(iout,'(f8.0)') ctabl(1,2,ialgo) write(iout,'(f8.0)') ctabl(2,2,ialgo) write(iout,'(f8.3)') r write(iout,'(f8.2)') xmae write(iout,'(f8.2)') rmse write(iout,'(f8.2)') rmsen write(iout,'(f8.2)') cgag write(iout,'(f8.2)') calg write(iout,'(f8.2)') rmaxg write(iout,'(f8.2)') rmaxa write(iout,'(f8.2)') gagemn write(iout,'(f8.2)') algomn write(iout,'(a8)') title(ialgo) c write(iout_gr,rec=(ndays-1)*10+7) r write(iout_gr,rec=(ndays-1)*10+8) xmae write(iout_gr,rec=(ndays-1)*10+9) rmse write(iout_gr,rec=(ndays-1)*10+10) diff close (iout) close (iout_gr) 500 continue return end c c------------------------------------------------------------------------------------- c subroutine vradar(xmsng,nx,ny,ix1,ix2,jy1,jy2,nalgo,algo,gag,gagn, x istag,ctabl,calgo,jdate,iexist,out_ascii,out_grads, x iout,iout_gr,a,b,ndays,title) c c------- Validates estimates using RADAR data as the reference standard c real*4 algo(nx,ny,nalgo), gag(nx,ny), xmask(nx,ny), x ctabl(2,2,nalgo), a(nx*ny), b(nx*ny), gagn(nx,ny) c integer*4 iexist(nalgo), istag(nalgo) logical there c character*72 out_ascii character*72 out_grads character*8 jdate character*6 calgo (nalgo) character*8 title(nalgo) call contin(xmsng,nx,ny,ix1,ix2,jy1,jy2,nalgo,algo,gag,gagn,istag, x ctabl) do 500 ialgo=1,nalgo write(*,*) 'in "vradar", "iout" = ',iout write(iout,*) calgo(ialgo) if(calgo(ialgo).eq.'radar') go to 500 out_ascii(19:24)=calgo(ialgo) out_ascii(26:31)=calgo(ialgo) out_ascii(42:49)=jdate out_grads(31:36)=calgo(ialgo) out_grads(37:42)='.grads' open (iout,file=out_ascii) open (iout_gr,file=out_grads,access='direct',recl=1) write(iout,*) calgo(ialgo) if(iexist(ialgo).eq.0) then do istat=1,10 write(iout_gr,rec=(ndays-1)*10+istat) xmsng enddo do istat=2,25 write(iout,'(f8.0)') xmsng enddo go to 500 endif c c write(*,*) write(*,*) '********************************************' write(*,*) write(*,*) write(*,*) ' RADAR' write(*,*) write(*,*) ' >=1 <1' write(*,*) write(*,'(8x,''>=1'',2f8.0)') (ctabl(i,1,ialgo),i=1,2) write(*,*) calgo(ialgo) write(*,'(9x,''<1'',2f8.0)') (ctabl(i,2,ialgo),i=1,2) write(*,*) c call precor(algo,gag,gagn,istag,nx,ny,ix1,ix2,jy1,jy2,nalgo,ialgo, x xmsng,a,b,npts) call corr(a,b,npts,r,slope,yint,0) write(*,'(''Correlation :'',f8.3)') r c call stats(algo,gag,gagn,istag,nx,ny,ix1,ix2,jy1,jy2,nalgo,ialgo, x xmsng,xmae,rmse,cgag,calg, rmaxg,rmaxa,algomn,gagemn) write(*,'(''Mean Abs Error :'',f8.3)') xmae write(*,'(''RMSE :'',f8.3)') rmse write(*,'(''Normalized (gauge)RMSE :'',f8.3)') rmse/gagemn c call scores(nalgo,ialgo,ctabl,calgo,iout,iout_gr,ndays,iexist) write(*,'(''Conditional rain rate :'',2f8.1)') cgag,calg write(*,'(''Maximum rain value :'',2f8.1)') rmaxg,rmaxa write(*,'(''Mean rain :'',2f8.1)') gagemn,algomn c c.... Write statistics to file for Graphics c write(iout,'(f8.0)') ctabl(1,1,ialgo) write(iout,'(f8.0)') ctabl(2,1,ialgo) write(iout,'(f8.0)') ctabl(1,2,ialgo) write(iout,'(f8.0)') ctabl(2,2,ialgo) write(iout,'(f8.3)') r write(iout,'(f8.2)') xmae write(iout,'(f8.2)') rmse write(iout,'(f8.2)') rmse/gagemn write(iout,'(f8.2)') cgag write(iout,'(f8.2)') calg write(iout,'(f8.2)') rmaxg write(iout,'(f8.2)') rmaxa write(iout,'(f8.2)') gagemn write(iout,'(f8.2)') algomn write(iout,'(a8)') title(ialgo) c write(iout_gr,rec=(ndays-1)*10+7) r write(iout_gr,rec=(ndays-1)*10+8) xmae write(iout_gr,rec=(ndays-1)*10+9) rmse write(iout_gr,rec=(ndays-1)*10+10) algomn-gagemn close (iout) close (iout_gr) 500 continue return end