program routinecv !====================================================================== ! compile as xlf90 -L/nwprod/lib -lw3_4 ! Choose: a) number of cross-validation sets ncvs ! b) total number of Z times per day ncycs. polulate ! cyc accordingly ! c) total number of days nt. can just decide to ! choose a large, fixed value, eg. nt=100. !====================================================================== implicit none ! character(80),parameter::filename='rawstats.dat_for_rmuse_pos1' character(80) filename integer(4),parameter::nvars=9+4+8 !p,t,q,u,v,w,w2,wdir,td + gust,vis,pblh,dist + !wspd10m,td2m,mxtm,mitm,pmsl,howv,tcamt,lcbas !note: w1 and w2 account for the two !distinct manners used to compute !rmse for total wind include 'params_nt_ncycs_ncvs.incl' ! integer(4),parameter::nt=8 !total number of days ! integer(4),parameter::ncycs=24 !total number of Z times per day ! integer(4),parameter::ncvs=1 !total number of cross-validation sets integer(4) k,m,n,nn integer(4) nobsin(nvars,3) !the 3 is for sigges,siganl real(4) biasin(nvars,3),absein(nvars,3),rmsein(nvars,3) !and sigcress character(20) cvar(nvars) character(2) cyc(ncycs) !store hours 00,01,02,03,04,05, ... 20,21,22,23 character(8) pdy(nt) character(120) cpath,cpathfull character(10) cdate character(6) cmonth character(4) cyear character(2) ch2 integer(4) n0,n1,n2,n3,npts,ncount,ndcyc integer(4) nfound,nmiss logical fexist,lcheck,goodcycle integer(4) nobs(ncvs,ncycs,nt,nvars) integer(4) nobs1(ncvs,ncycs,nt,nvars) integer(4) nobs2(ncvs,ncycs,nt,nvars) integer(4) nobs3(ncvs,ncycs,nt,nvars) real(4),dimension(ncvs,ncycs,nt,nvars):: & vbias1,vrmse1,vabse1, & !guess vbias2,vrmse2,vabse2, & !2dvar vbias3,vrmse3,vabse3 !cressman real(4),dimension(ncycs,nt,nvars):: & vbias1_hdvar,vrmse1_hdvar,vabse1_hdvar, & !guess /per cyc,day and cvar vbias2_hdvar,vrmse2_hdvar,vabse2_hdvar, & !2dvar /per cyc,day and cvar vbias3_hdvar,vrmse3_hdvar,vabse3_hdvar !cressman /per cyc,day and cvar real(4),dimension(ncycs,nvars):: & vbias1_hvar,vrmse1_hvar,vabse1_hvar, & !guess /per cyc and cvar vbias2_hvar,vrmse2_hvar,vabse2_hvar, & !2dvar /per cyc and cvar vbias3_hvar,vrmse3_hvar,vabse3_hvar !cressman /per cyc and cvar real(4),dimension(nt,nvars):: & vbias1_dvar,vrmse1_dvar,vabse1_dvar, & !guess /per day and cvar vbias2_dvar,vrmse2_dvar,vabse2_dvar, & !2dvar /per day and cvar vbias3_dvar,vrmse3_dvar,vabse3_dvar !cressman /per day and cvar real(4),dimension(nvars):: & vbias1_var,vrmse1_var,vabse1_var, & !guess /per cvar vbias2_var,vrmse2_var,vabse2_var, & !2dvar /per cvar vbias3_var,vrmse3_var,vabse3_var !cressman /per cvar real(4),parameter::spval=-1111. real(4),parameter::r_0001=1.e-03 real(4) spval_hdvar(ncycs,nt,nvars) real(4) spval_hvar(ncycs,nvars) real(4) spval_dvar(nt,nvars) real(4) spval_var(nvars) real(4),dimension(nvars):: & vbias1_hdvaravg,vrmse1_hdvaravg,vabse1_hdvaravg, & vbias2_hdvaravg,vrmse2_hdvaravg,vabse2_hdvaravg, & vbias3_hdvaravg,vrmse3_hdvaravg,vabse3_hdvaravg real(4) s1,s2,s3 integer(4) iyear0,imonth0,iday0 integer(4) iyear,imonth,iday,idaywk,idayyyr integer(4) iw3jdn,jldayn namelist/statscycleinfo/cpath,filename,iyear0,imonth0,iday0 ! data cpath /'meso/noscrub/wx20pm/CVRAWDATA_2p5'/ data filename /'rawstats.dat_for_rmuse_pos1'/ data iyear0 /2010/ data imonth0 /6/ data iday0 /29/ data cvar(1) /'p'/ data cvar(2) /'t'/ data cvar(3) /'q'/ data cvar(4) /'u'/ data cvar(5) /'v'/ data cvar(6) /'w'/ data cvar(7) /'w2'/ data cvar(8) /'wd'/ data cvar(9) /'td'/ data cvar(10) /'gust'/ data cvar(11) /'vis'/ data cvar(12) /'pblh'/ data cvar(13) /'dist'/ data cvar(14) /'wspd10m'/ data cvar(15) /'td2m'/ data cvar(16) /'mxtm'/ data cvar(17) /'mitm'/ data cvar(18) /'pmsl'/ data cvar(19) /'howv'/ data cvar(20) /'tcamt'/ data cvar(21) /'lcbas'/ if (ncycs==24) ndcyc=1 if (ncycs==8) ndcyc=3 do m=1,ncycs write(cyc(m),'(i2.2)') (m-1)*ndcyc print*,'m,cyc(m)=',cyc(m) enddo cpath='/ptmpp1/Manuel.Pondeca/TEST/2010/' print*,'nt,ncycs,ncvs,=',nt,ncycs,ncvs print*,' ' inquire (file='statscycleinfo_input',exist=fexist) if (fexist) then open (55,file='statscycleinfo_input',form='formatted') read(55,statscycleinfo) close(55) endif print*,'iyear0,imonth0,iday0=',iyear0,imonth0,iday0 jldayn=iw3jdn(iyear0,imonth0,iday0) do n=1,nt if (n > 1) jldayn=jldayn+1 print*,'jldayn=',jldayn call w3fs26(jldayn,iyear,imonth,iday,idaywk,idayyyr) write(pdy(n)(1:4),'(i4.4)') iyear write(pdy(n)(5:6),'(i2.2)') imonth write(pdy(n)(7:8),'(i2.2)') iday print*,'n,pdy(n)=',n,trim(pdy(n)) enddo !***************************************************************************************** nobs1=0 ; nobs2=0 ; nobs3=0 vbias1=0. ; vbias2=0. ; vbias3=0. vrmse1=0. ; vrmse2=0. ; vrmse3=0. vabse1=0. ; vabse2=0. ; vabse3=0. nfound=0 nmiss=0 do n=1,nt do m=1,ncycs cdate=pdy(n)//cyc(m) cmonth=pdy(n)(1:6) cyear=pdy(n)(1:4) do k=1,ncvs write(ch2,'(i2.2)') k ! cpathfull=trim(cpath)//pdy(n)//'/'//cdate//'/cvset_'//ch2 ! cpathfull=trim(cpath)//pdy(n)//'/'//cdate cpathfull=trim(cpath)//cyear//'/'//cmonth//'/'//pdy(n)//'/'//cdate print*,'day,cdate,cvset,input file=',n,m,k,trim(cpathfull)//'/'//trim(filename) inquire (file=trim(cpathfull)//'/'//trim(filename),exist=fexist) if (fexist) then call read_rawstats(cpathfull,filename,nobsin, & biasin,absein,rmsein,nvars,3,goodcycle) if (.not.goodcycle) cycle nfound=nfound+1 do nn=1,nvars nobs1(k,m,n,nn)=nobsin(nn,1) nobs2(k,m,n,nn)=nobsin(nn,2) nobs3(k,m,n,nn)=nobsin(nn,3) vbias1(k,m,n,nn)=biasin(nn,1) vbias2(k,m,n,nn)=biasin(nn,2) vbias3(k,m,n,nn)=biasin(nn,3) vrmse1(k,m,n,nn)=rmsein(nn,1) vrmse2(k,m,n,nn)=rmsein(nn,2) vrmse3(k,m,n,nn)=rmsein(nn,3) vabse1(k,m,n,nn)=absein(nn,1) vabse2(k,m,n,nn)=absein(nn,2) vabse3(k,m,n,nn)=absein(nn,3) enddo else nmiss=nmiss+1 print*,'input file does not exist. nobs will be zero' endif enddo enddo enddo print*,'number of input files expected=',nt*ncycs*ncvs print*,'number of input files found=',nfound print*,'number of missing input files=',nmiss print*,'note that number of missing input files include files found with invalid data' ! for each cycle and cv-set, make sure that the ! number of obs used to verify a given variable is ! the same for all models (sigges,siganl,sigcress) lcheck=.true. do n=1,nt do m=1,ncycs do k=1,ncvs do nn=1,nvars n1=nobs1(k,m,n,nn)-nobs2(k,m,n,nn) n2=nobs1(k,m,n,nn)-nobs3(k,m,n,nn) n3=nobs2(k,m,n,nn)-nobs3(k,m,n,nn) if (n1/=0 .or. n2/=0 .or. n3/=0) then lcheck=.false. print*,' trouble in ob count: n,m,k,nn,n1,n2,n3=',& n,m,k,nn,n1,n2,n3 endif enddo enddo enddo enddo if (.not.lcheck) then print*,' aborting ... ' call abort endif nobs=nobs1 ! average over number of cross-validation sets vbias1_hdvar(:,:,:)=0. vbias2_hdvar(:,:,:)=0. vbias3_hdvar(:,:,:)=0. vrmse1_hdvar(:,:,:)=0. vrmse2_hdvar(:,:,:)=0. vrmse3_hdvar(:,:,:)=0. vabse1_hdvar(:,:,:)=0. vabse2_hdvar(:,:,:)=0. vabse3_hdvar(:,:,:)=0. spval_hdvar(:,:,:)=0. do nn=1,nvars do n=1,nt do m=1,ncycs npts=0 do k=1,ncvs npts=npts+nobs(k,m,n,nn) enddo if (npts==0) spval_hdvar(m,n,nn)=spval npts=max(1,npts) s1=0. s2=0. s3=0. do k=1,ncvs n0=nobs(k,m,n,nn) vbias1_hdvar(m,n,nn)=vbias1_hdvar(m,n,nn)+n0*vbias1(k,m,n,nn)/float(npts) vbias2_hdvar(m,n,nn)=vbias2_hdvar(m,n,nn)+n0*vbias2(k,m,n,nn)/float(npts) vbias3_hdvar(m,n,nn)=vbias3_hdvar(m,n,nn)+n0*vbias3(k,m,n,nn)/float(npts) vabse1_hdvar(m,n,nn)=vabse1_hdvar(m,n,nn)+n0*vabse1(k,m,n,nn)/float(npts) vabse2_hdvar(m,n,nn)=vabse2_hdvar(m,n,nn)+n0*vabse2(k,m,n,nn)/float(npts) vabse3_hdvar(m,n,nn)=vabse3_hdvar(m,n,nn)+n0*vabse3(k,m,n,nn)/float(npts) s1=s1+n0*vrmse1(k,m,n,nn)*vrmse1(k,m,n,nn)/float(npts) s2=s2+n0*vrmse2(k,m,n,nn)*vrmse2(k,m,n,nn)/float(npts) s3=s3+n0*vrmse3(k,m,n,nn)*vrmse3(k,m,n,nn)/float(npts) enddo if (s1 > 0.) vrmse1_hdvar(m,n,nn)=sqrt(s1) if (s2 > 0.) vrmse2_hdvar(m,n,nn)=sqrt(s2) if (s3 > 0.) vrmse3_hdvar(m,n,nn)=sqrt(s3) enddo enddo enddo ! average over number of cross-validation sets and number of days vbias1_hvar(:,:)=0. vbias2_hvar(:,:)=0. vbias3_hvar(:,:)=0. vrmse1_hvar(:,:)=0. vrmse2_hvar(:,:)=0. vrmse3_hvar(:,:)=0. vabse1_hvar(:,:)=0. vabse2_hvar(:,:)=0. vabse3_hvar(:,:)=0. spval_hvar(:,:)=0. do nn=1,nvars do m=1,ncycs npts=0 do n=1,nt do k=1,ncvs npts=npts+nobs(k,m,n,nn) enddo enddo if (npts==0) spval_hvar(m,nn)=spval npts=max(1,npts) s1=0. s2=0. s3=0. do n=1,nt do k=1,ncvs n0=nobs(k,m,n,nn) vbias1_hvar(m,nn)=vbias1_hvar(m,nn)+n0*vbias1(k,m,n,nn)/float(npts) vbias2_hvar(m,nn)=vbias2_hvar(m,nn)+n0*vbias2(k,m,n,nn)/float(npts) Vbias3_hvar(m,nn)=vbias3_hvar(m,nn)+n0*vbias3(k,m,n,nn)/float(npts) vabse1_hvar(m,nn)=vabse1_hvar(m,nn)+n0*vabse1(k,m,n,nn)/float(npts) vabse2_hvar(m,nn)=vabse2_hvar(m,nn)+n0*vabse2(k,m,n,nn)/float(npts) Vabse3_hvar(m,nn)=vabse3_hvar(m,nn)+n0*vabse3(k,m,n,nn)/float(npts) s1=s1+n0*vrmse1(k,m,n,nn)*vrmse1(k,m,n,nn)/float(npts) s2=s2+n0*vrmse2(k,m,n,nn)*vrmse2(k,m,n,nn)/float(npts) s3=s3+n0*vrmse3(k,m,n,nn)*vrmse3(k,m,n,nn)/float(npts) enddo enddo if (s1 > 0.) vrmse1_hvar(m,nn)=sqrt(s1) if (s2 > 0.) vrmse2_hvar(m,nn)=sqrt(s2) if (s3 > 0.) vrmse3_hvar(m,nn)=sqrt(s3) enddo enddo ! average over number of cross-validation sets and cycles vbias1_dvar(:,:)=0. vbias2_dvar(:,:)=0. vbias3_dvar(:,:)=0. vrmse1_dvar(:,:)=0. vrmse2_dvar(:,:)=0. vrmse3_dvar(:,:)=0. vabse1_dvar(:,:)=0. vabse2_dvar(:,:)=0. vabse3_dvar(:,:)=0. spval_dvar(:,:)=0. do nn=1,nvars do n=1,nt npts=0 do m=1,ncycs do k=1,ncvs npts=npts+nobs(k,m,n,nn) enddo enddo if (npts==0) spval_dvar(n,nn)=spval npts=max(1,npts) s1=0. s2=0. s3=0. do m=1,ncycs do k=1,ncvs n0=nobs(k,m,n,nn) vbias1_dvar(n,nn)=vbias1_dvar(n,nn)+n0*vbias1(k,m,n,nn)/float(npts) vbias2_dvar(n,nn)=vbias2_dvar(n,nn)+n0*vbias2(k,m,n,nn)/float(npts) vbias3_dvar(n,nn)=vbias3_dvar(n,nn)+n0*vbias3(k,m,n,nn)/float(npts) vabse1_dvar(n,nn)=vabse1_dvar(n,nn)+n0*vabse1(k,m,n,nn)/float(npts) vabse2_dvar(n,nn)=vabse2_dvar(n,nn)+n0*vabse2(k,m,n,nn)/float(npts) vabse3_dvar(n,nn)=vabse3_dvar(n,nn)+n0*vabse3(k,m,n,nn)/float(npts) s1=s1+n0*vrmse1(k,m,n,nn)*vrmse1(k,m,n,nn)/float(npts) s2=s2+n0*vrmse2(k,m,n,nn)*vrmse2(k,m,n,nn)/float(npts) s3=s3+n0*vrmse3(k,m,n,nn)*vrmse3(k,m,n,nn)/float(npts) enddo enddo if (s1 > 0.) vrmse1_dvar(n,nn)=sqrt(s1) if (s2 > 0.) vrmse2_dvar(n,nn)=sqrt(s2) if (s3 > 0.) vrmse3_dvar(n,nn)=sqrt(s3) enddo enddo ! average over number of cross-validation sets and number of days and times vbias1_var(:)=0. vbias2_var(:)=0. vbias3_var(:)=0. vrmse1_var(:)=0. vrmse2_var(:)=0. vrmse3_var(:)=0. vabse1_var(:)=0. vabse2_var(:)=0. vabse3_var(:)=0. spval_var(:)=0. do nn=1,nvars npts=0 do m=1,ncycs do n=1,nt do k=1,ncvs npts=npts+nobs(k,m,n,nn) enddo enddo enddo if (npts==0) spval_var(nn)=spval npts=max(1,npts) s1=0. s2=0. s3=0. do m=1,ncycs do n=1,nt do k=1,ncvs n0=nobs(k,m,n,nn) vbias1_var(nn)=vbias1_var(nn)+n0*vbias1(k,m,n,nn)/float(npts) vbias2_var(nn)=vbias2_var(nn)+n0*vbias2(k,m,n,nn)/float(npts) vbias3_var(nn)=vbias3_var(nn)+n0*vbias3(k,m,n,nn)/float(npts) vabse1_var(nn)=vabse1_var(nn)+n0*vabse1(k,m,n,nn)/float(npts) vabse2_var(nn)=vabse2_var(nn)+n0*vabse2(k,m,n,nn)/float(npts) vabse3_var(nn)=vabse3_var(nn)+n0*vabse3(k,m,n,nn)/float(npts) s1=s1+n0*vrmse1(k,m,n,nn)*vrmse1(k,m,n,nn)/float(npts) s2=s2+n0*vrmse2(k,m,n,nn)*vrmse2(k,m,n,nn)/float(npts) s3=s3+n0*vrmse3(k,m,n,nn)*vrmse3(k,m,n,nn)/float(npts) enddo enddo enddo if (s1 > 0.) vrmse1_var(nn)=sqrt(s1) if (s2 > 0.) vrmse2_var(nn)=sqrt(s2) if (s3 > 0.) vrmse3_var(nn)=sqrt(s3) enddo ! if no obs available, write spval to correspondent entry in array vbias1_hdvar = vbias1_hdvar + spval_hdvar vbias2_hdvar = vbias2_hdvar + spval_hdvar vbias3_hdvar = vbias3_hdvar + spval_hdvar vrmse1_hdvar = vrmse1_hdvar + spval_hdvar vrmse2_hdvar = vrmse2_hdvar + spval_hdvar vrmse3_hdvar = vrmse3_hdvar + spval_hdvar vabse1_hdvar = vabse1_hdvar + spval_hdvar vabse2_hdvar = vabse2_hdvar + spval_hdvar vabse3_hdvar = vabse3_hdvar + spval_hdvar vbias1_hvar = vbias1_hvar + spval_hvar vbias2_hvar = vbias2_hvar + spval_hvar vbias3_hvar = vbias3_hvar + spval_hvar vrmse1_hvar = vrmse1_hvar + spval_hvar vrmse2_hvar = vrmse2_hvar + spval_hvar vrmse3_hvar = vrmse3_hvar + spval_hvar vabse1_hvar = vabse1_hvar + spval_hvar vabse2_hvar = vabse2_hvar + spval_hvar vabse3_hvar = vabse3_hvar + spval_hvar vbias1_dvar = vbias1_dvar + spval_dvar vbias2_dvar = vbias2_dvar + spval_dvar vbias3_dvar = vbias3_dvar + spval_dvar vrmse1_dvar = vrmse1_dvar + spval_dvar vrmse2_dvar = vrmse2_dvar + spval_dvar vrmse3_dvar = vrmse3_dvar + spval_dvar vabse1_dvar = vabse1_dvar + spval_dvar vabse2_dvar = vabse2_dvar + spval_dvar vabse3_dvar = vabse3_dvar + spval_dvar vbias1_var = vbias1_var + spval_var vbias2_var = vbias2_var + spval_var vbias3_var = vbias3_var + spval_var vrmse1_var = vrmse1_var + spval_var vrmse2_var = vrmse2_var + spval_var vrmse3_var = vrmse3_var + spval_var vabse1_var = vabse1_var + spval_var vabse2_var = vabse2_var + spval_var vabse3_var = vabse3_var + spval_var !==>compute avergage of time series that has cycle times in the absciss vbias1_hdvaravg(:)=0. ; vrmse1_hdvaravg(:)=0. ; vabse1_hdvaravg(:)=0. vbias2_hdvaravg(:)=0. ; vrmse2_hdvaravg(:)=0. ; vabse2_hdvaravg(:)=0. vbias3_hdvaravg(:)=0. ; vrmse3_hdvaravg(:)=0. ; vabse3_hdvaravg(:)=0. do nn=1,nvars ncount=0 do n=1,nt do m=1,ncycs if (abs(spval_hdvar(m,n,nn)-spval) .gt. r_0001) then ncount=ncount+1 vbias1_hdvaravg(nn)=vbias1_hdvaravg(nn)+vbias1_hdvar(m,n,nn) vbias2_hdvaravg(nn)=vbias2_hdvaravg(nn)+vbias2_hdvar(m,n,nn) vbias3_hdvaravg(nn)=vbias3_hdvaravg(nn)+vbias3_hdvar(m,n,nn) vabse1_hdvaravg(nn)=vabse1_hdvaravg(nn)+vabse1_hdvar(m,n,nn) vabse2_hdvaravg(nn)=vabse2_hdvaravg(nn)+vabse2_hdvar(m,n,nn) vabse3_hdvaravg(nn)=vabse3_hdvaravg(nn)+vabse3_hdvar(m,n,nn) vrmse1_hdvaravg(nn)=vrmse1_hdvaravg(nn)+vrmse1_hdvar(m,n,nn) vrmse2_hdvaravg(nn)=vrmse2_hdvaravg(nn)+vrmse2_hdvar(m,n,nn) vrmse3_hdvaravg(nn)=vrmse3_hdvaravg(nn)+vrmse3_hdvar(m,n,nn) endif enddo enddo if (ncount > 0) then vbias1_hdvaravg(nn)=vbias1_hdvaravg(nn)/float(ncount) vbias2_hdvaravg(nn)=vbias2_hdvaravg(nn)/float(ncount) vbias3_hdvaravg(nn)=vbias3_hdvaravg(nn)/float(ncount) vabse1_hdvaravg(nn)=vabse1_hdvaravg(nn)/float(ncount) vabse2_hdvaravg(nn)=vabse2_hdvaravg(nn)/float(ncount) vabse3_hdvaravg(nn)=vabse3_hdvaravg(nn)/float(ncount) vrmse1_hdvaravg(nn)=vrmse1_hdvaravg(nn)/float(ncount) vrmse2_hdvaravg(nn)=vrmse2_hdvaravg(nn)/float(ncount) vrmse3_hdvaravg(nn)=vrmse3_hdvaravg(nn)/float(ncount) endif enddo do nn=1,nvars print*,'statistics for ',trim(cvar(nn)) print*,' ' do m=1,ncycs print*,'hour,abse1,abse2,abse3=',cyc(m),vabse1_hvar(m,nn),vabse2_hvar(m,nn),vabse3_hvar(m,nn) enddo print*,'avg,abse1,abse2,abse3=','avg',vabse1_var(nn),vabse2_var(nn),vabse3_var(nn) print*,' ' do m=1,ncycs print*,'hour,bias1,bias2,bias3=',cyc(m),vbias1_hvar(m,nn),vbias2_hvar(m,nn),vbias3_hvar(m,nn) enddo print*,'avg,bias1,bias2,bias3=','avg',vbias1_var(nn),vbias2_var(nn),vbias3_var(nn) print*,' ' do m=1,ncycs print*,'hour,rmse1,rmse2,rmse3=',cyc(m),vrmse1_hvar(m,nn),vrmse2_hvar(m,nn),vrmse3_hvar(m,nn) enddo print*,'avg,rmse1,rmse2,rmse3=','avg',vrmse1_var(nn),vrmse2_var(nn),vrmse3_var(nn) print*,'=====================================================' print*,'=====================================================' enddo 112 format(a4,f12.6,2x,f12.6,6x,f12.6) !==>write out to file for grads display !============================================================ !==> averages by cycle time for each variable !============================================================ open (11,file='vbias1_hvar.dat',form='unformatted') write(11) vbias1_hvar close(11) open (11,file='vbias2_hvar.dat',form='unformatted') write(11) vbias2_hvar close(11) open (11,file='vbias3_hvar.dat',form='unformatted') write(11) vbias3_hvar close(11) open (11,file='vrmse1_hvar.dat',form='unformatted') write(11) vrmse1_hvar close(11) open (11,file='vrmse2_hvar.dat',form='unformatted') write(11) vrmse2_hvar close(11) open (11,file='vrmse3_hvar.dat',form='unformatted') write(11) vrmse3_hvar close(11) open (11,file='vabse1_hvar.dat',form='unformatted') write(11) vabse1_hvar close(11) open (11,file='vabse2_hvar.dat',form='unformatted') write(11) vabse2_hvar close(11) open (11,file='vabse3_hvar.dat',form='unformatted') write(11) vabse3_hvar close(11) ! !============================================================ !==> time series with days in the absciss ordinate !============================================================ open (11,file='vbias1_dvar.dat',form='unformatted') write(11) vbias1_dvar close(11) open (11,file='vbias2_dvar.dat',form='unformatted') write(11) vbias2_dvar close(11) open (11,file='vbias3_dvar.dat',form='unformatted') write(11) vbias3_dvar close(11) open (11,file='vrmse1_dvar.dat',form='unformatted') write(11) vrmse1_dvar close(11) open (11,file='vrmse2_dvar.dat',form='unformatted') write(11) vrmse2_dvar close(11) open (11,file='vrmse3_dvar.dat',form='unformatted') write(11) vrmse3_dvar close(11) open (11,file='vabse1_dvar.dat',form='unformatted') write(11) vabse1_dvar close(11) open (11,file='vabse2_dvar.dat',form='unformatted') write(11) vabse2_dvar close(11) open (11,file='vabse3_dvar.dat',form='unformatted') write(11) vabse3_dvar close(11) ! !============================================================ !==> time series with cycle times in the absciss ordinate !============================================================ open (11,file='vbias1_hdvar.dat',form='unformatted') write(11) vbias1_hdvar close(11) open (11,file='vbias2_hdvar.dat',form='unformatted') write(11) vbias2_hdvar close(11) open (11,file='vbias3_hdvar.dat',form='unformatted') write(11) vbias3_hdvar close(11) open (11,file='vrmse1_hdvar.dat',form='unformatted') write(11) vrmse1_hdvar close(11) open (11,file='vrmse2_hdvar.dat',form='unformatted') write(11) vrmse2_hdvar close(11) open (11,file='vrmse3_hdvar.dat',form='unformatted') write(11) vrmse3_hdvar close(11) open (11,file='vabse1_hdvar.dat',form='unformatted') write(11) vabse1_hdvar close(11) open (11,file='vabse2_hdvar.dat',form='unformatted') write(11) vabse2_hdvar close(11) open (11,file='vabse3_hdvar.dat',form='unformatted') write(11) vabse3_hdvar close(11) call dumpavgs(nvars,cvar, & vbias1_hdvaravg,vrmse1_hdvaravg,vabse1_hdvaravg, & vbias2_hdvaravg,vrmse2_hdvaravg,vabse2_hdvaravg, & vbias3_hdvaravg,vrmse3_hdvaravg,vabse3_hdvaravg ,& iyear0,imonth0,iday0,filename) end !=============================================================== !=============================================================== subroutine dumpavgs(nvars,cvar, & vbias1_hdvaravg,vrmse1_hdvaravg,vabse1_hdvaravg, & vbias2_hdvaravg,vrmse2_hdvaravg,vabse2_hdvaravg, & vbias3_hdvaravg,vrmse3_hdvaravg,vabse3_hdvaravg ,& iyear0,imonth0,iday0,filename) implicit none character(3),parameter::cyc0='00Z' integer(4),intent(in):: nvars integer(4),intent(in):: iyear0,imonth0,iday0 real(4),dimension(nvars),intent(in):: & vbias1_hdvaravg,vrmse1_hdvaravg,vabse1_hdvaravg, & vbias2_hdvaravg,vrmse2_hdvaravg,vabse2_hdvaravg, & vbias3_hdvaravg,vrmse3_hdvaravg,vabse3_hdvaravg character(80),intent(in):: filename character(20),intent(in):: cvar(nvars) integer(4) nlen character(3) cmonth0 character(60) starttime character(4) cpos1neg9 integer(4) n real(4) scale if (imonth0==1) cmonth0='JAN' if (imonth0==2) cmonth0='FEB' if (imonth0==3) cmonth0='MAR' if (imonth0==4) cmonth0='APR' if (imonth0==5) cmonth0='MAY' if (imonth0==6) cmonth0='JUN' if (imonth0==7) cmonth0='JUL' if (imonth0==8) cmonth0='AUG' if (imonth0==9) cmonth0='SEP' if (imonth0==10) cmonth0='OCT' if (imonth0==11) cmonth0='NOV' if (imonth0==12) cmonth0='DEC' write(starttime,'(a3,1x,i2.2,1x,a3,1x,i4)') cyc0,iday0,cmonth0,iyear0 print*,'starttime=','"',trim(starttime),'"' nlen=len_trim(filename) cpos1neg9(1:4)=filename(nlen-3:nlen) open (12,file='allavgs_of_hdvar.txt',form='formatted') write(12,*) 'starttime=','"',trim(starttime),'"' write(12,*) 'filetype=','"',trim(cpos1neg9),'"' write(12,*) ' ' write(12,*)'p_bckgbiasavg="AVG =',vbias1_hdvaravg(1), 'hPa"' write(12,*)'t_bckgbiasavg="AVG =',vbias1_hdvaravg(2), 'K"' write(12,*)'q_bckgbiasavg="AVG =',vbias1_hdvaravg(3), 'g/kg"' write(12,*)'u_bckgbiasavg="AVG =',vbias1_hdvaravg(4), 'm/s"' write(12,*)'v_bckgbiasavg="AVG =',vbias1_hdvaravg(5), 'm/s"' write(12,*)'w_bckgbiasavg="AVG =',vbias1_hdvaravg(6), 'm/s"' write(12,*)'w2_bckgbiasavg="AVG =',vbias1_hdvaravg(7), 'm/s"' write(12,*)'wd_bckgbiasavg="AVG =',vbias1_hdvaravg(8), 'deg"' write(12,*)'td_bckgbiasavg="AVG =',vbias1_hdvaravg(9), 'K"' write(12,*)'p_bckgabseavg="AVG =',vabse1_hdvaravg(1), 'hPa"' write(12,*)'t_bckgabseavg="AVG =',vabse1_hdvaravg(2), 'K"' write(12,*)'q_bckgabseavg="AVG =',vabse1_hdvaravg(3), 'g/kg"' write(12,*)'u_bckgabseavg="AVG =',vabse1_hdvaravg(4), 'm/s"' write(12,*)'v_bckgabseavg="AVG =',vabse1_hdvaravg(5), 'm/s"' write(12,*)'w_bckgabseavg="AVG =',vabse1_hdvaravg(6), 'm/s"' write(12,*)'w2_bckgabseavg="AVG =',vabse1_hdvaravg(7), 'm/s"' write(12,*)'wd_bckgabseavg="AVG =',vabse1_hdvaravg(8), 'deg"' write(12,*)'td_bckgabseavg="AVG =',vabse1_hdvaravg(9), 'K"' write(12,*)'p_bckgrmseavg="AVG =',vrmse1_hdvaravg(1), 'hPa"' write(12,*)'t_bckgrmseavg="AVG =',vrmse1_hdvaravg(2), 'K"' write(12,*)'q_bckgrmseavg="AVG =',vrmse1_hdvaravg(3), 'g/kg"' write(12,*)'u_bckgrmseavg="AVG =',vrmse1_hdvaravg(4), 'm/s"' write(12,*)'v_bckgrmseavg="AVG =',vrmse1_hdvaravg(5), 'm/s"' write(12,*)'w_bckgrmseavg="AVG =',vrmse1_hdvaravg(6), 'm/s"' write(12,*)'w2_bckgrmseavg="AVG =',vrmse1_hdvaravg(7), 'm/s"' write(12,*)'wd_bckgrmseavg="AVG =',vrmse1_hdvaravg(8), 'deg"' write(12,*)'td_bckgrmseavg="AVG =',vrmse1_hdvaravg(9), 'K"' write(12,*) ' ' write(12,*)'p_anlbiasavg="AVG =',vbias2_hdvaravg(1), 'hPa"' write(12,*)'t_anlbiasavg="AVG =',vbias2_hdvaravg(2), 'K"' write(12,*)'q_anlbiasavg="AVG =',vbias2_hdvaravg(3), 'g/kg"' write(12,*)'u_anlbiasavg="AVG =',vbias2_hdvaravg(4), 'm/s"' write(12,*)'v_anlbiasavg="AVG =',vbias2_hdvaravg(5), 'm/s"' write(12,*)'w_anlbiasavg="AVG =',vbias2_hdvaravg(6), 'm/s"' write(12,*)'w2_anlbiasavg="AVG =',vbias2_hdvaravg(7), 'm/s"' write(12,*)'wd_anlbiasavg="AVG =',vbias2_hdvaravg(8), 'deg"' write(12,*)'td_anlbiasavg="AVG =',vbias2_hdvaravg(9), 'K"' write(12,*)'p_anlabseavg="AVG =',vabse2_hdvaravg(1), 'hPa"' write(12,*)'t_anlabseavg="AVG =',vabse2_hdvaravg(2), 'K"' write(12,*)'q_anlabseavg="AVG =',vabse2_hdvaravg(3), 'g/kg"' write(12,*)'u_anlabseavg="AVG =',vabse2_hdvaravg(4), 'm/s"' write(12,*)'v_anlabseavg="AVG =',vabse2_hdvaravg(5), 'm/s"' write(12,*)'w_anlabseavg="AVG =',vabse2_hdvaravg(6), 'm/s"' write(12,*)'w2_anlabseavg="AVG =',vabse2_hdvaravg(7), 'm/s"' write(12,*)'wd_anlabseavg="AVG =',vabse2_hdvaravg(8), 'deg"' write(12,*)'td_anlabseavg="AVG =',vabse2_hdvaravg(9), 'K"' write(12,*)'p_anlrmseavg="AVG =',vrmse2_hdvaravg(1), 'hPa"' write(12,*)'t_anlrmseavg="AVG =',vrmse2_hdvaravg(2), 'K"' write(12,*)'q_anlrmseavg="AVG =',vrmse2_hdvaravg(3), 'g/kg"' write(12,*)'u_anlrmseavg="AVG =',vrmse2_hdvaravg(4), 'm/s"' write(12,*)'v_anlrmseavg="AVG =',vrmse2_hdvaravg(5), 'm/s"' write(12,*)'w_anlrmseavg="AVG =',vrmse2_hdvaravg(6), 'm/s"' write(12,*)'w2_anlrmseavg="AVG =',vrmse2_hdvaravg(7), 'm/s"' write(12,*)'wd_anlrmseavg="AVG =',vrmse2_hdvaravg(8), 'deg"' write(12,*)'td_anlrmseavg="AVG =',vrmse2_hdvaravg(9), 'K"' write(12,*) ' ' write(12,*)'p_cressbiasavg="AVG =',vbias3_hdvaravg(1), 'hPa"' write(12,*)'t_cressbiasavg="AVG =',vbias3_hdvaravg(2), 'K"' write(12,*)'q_cressbiasavg="AVG =',vbias3_hdvaravg(3), 'g/kg"' write(12,*)'u_cressbiasavg="AVG =',vbias3_hdvaravg(4), 'm/s"' write(12,*)'v_cressbiasavg="AVG =',vbias3_hdvaravg(5), 'm/s"' write(12,*)'w_cressbiasavg="AVG =',vbias3_hdvaravg(6), 'm/s"' write(12,*)'w2_cressbiasavg="AVG =',vbias3_hdvaravg(7), 'm/s"' write(12,*)'wd_cressbiasavg="AVG =',vbias3_hdvaravg(8), 'deg"' write(12,*)'td_cressbiasavg="AVG =',vbias3_hdvaravg(9), 'K"' write(12,*)'p_cressabseavg="AVG =',vabse3_hdvaravg(1), 'hPa"' write(12,*)'t_cressabseavg="AVG =',vabse3_hdvaravg(2), 'K"' write(12,*)'q_cressabseavg="AVG =',vabse3_hdvaravg(3), 'g/kg"' write(12,*)'u_cressabseavg="AVG =',vabse3_hdvaravg(4), 'm/s"' write(12,*)'v_cressabseavg="AVG =',vabse3_hdvaravg(5), 'm/s"' write(12,*)'w_cressabseavg="AVG =',vabse3_hdvaravg(6), 'm/s"' write(12,*)'w2_cressabseavg="AVG =',vabse3_hdvaravg(7), 'm/s"' write(12,*)'wd_cressabseavg="AVG =',vabse3_hdvaravg(8), 'deg"' write(12,*)'td_cressabseavg="AVG =',vabse3_hdvaravg(9), 'K"' write(12,*)'p_cressrmseavg="AVG =',vrmse3_hdvaravg(1), 'hPa"' write(12,*)'t_cressrmseavg="AVG =',vrmse3_hdvaravg(2), 'K"' write(12,*)'q_cressrmseavg="AVG =',vrmse3_hdvaravg(3), 'g/kg"' write(12,*)'u_cressrmseavg="AVG =',vrmse3_hdvaravg(4), 'm/s"' write(12,*)'v_cressrmseavg="AVG =',vrmse3_hdvaravg(5), 'm/s"' write(12,*)'w_cressrmseavg="AVG =',vrmse3_hdvaravg(6), 'm/s"' write(12,*)'w2_cressrmseavg="AVG =',vrmse3_hdvaravg(7), 'm/s"' write(12,*)'wd_cressrmseavg="AVG =',vrmse3_hdvaravg(8), 'deg"' write(12,*)'td_cressrmseavg="AVG =',vrmse3_hdvaravg(9), 'K"' close(12) !==> averages do n=1,nvars dfactor=1. ; if (n==11) scale=1./1609.344 open (12,file=trim(cvar(n))'_avgs_of_hdvar.txt',form='formatted') write(12,'(a11,a15,a1)') 'starttime="',trim(starttime),'"' write(12,'(a9,a1,a4,a1)') 'filetype=','"',trim(cpos1neg9),'"' write(12,'(a22,f6.2,a1)') trim(cvar(n))'_bckgrmseavg="AVG =',vrmse1_hdvaravg(n)*scale, '"' write(12,'(a21,f6.2,a1)') trim(cvar(n))'_anlrmseavg="AVG =',vrmse2_hdvaravg(n)*scale, '"' write(12,'(a23,f6.2,a1)') trim(cvar(n))'_cressrmseavg="AVG =',vrmse3_hdvaravg(n)*scale, '"' write(12,'(a22,f6.2,a1)') trim(cvar(n))'_bckgbiasavg="AVG =',vbias1_hdvaravg(n)*scale, '"' write(12,'(a21,f6.2,a1)') trim(cvar(n))'_anlbiasavg="AVG =',vbias2_hdvaravg(n)*scale, '"' write(12,'(a23,f6.2,a1)') trim(cvar(n))'_cressbiasavg="AVG =',vbias3_hdvaravg(n)*scale, '"' write(12,'(a22,f6.2,a1)') trim(cvar(n))'_bckgabseavg="AVG =',vabse1_hdvaravg(n)*scale, '"' write(12,'(a21,f6.2,a1)') trim(cvar(n))'_anlabseavg="AVG =',vabse2_hdvaravg(n)*scale, '"' write(12,'(a23,f6.2,a1)') trim(cvar(n))'_cressabseavg="AVG =',vabse3_hdvaravg(n)*scale, '"' close(12) enddo return end !=============================================================== !=============================================================== subroutine read_rawstats(cpath,filename,nobsout, & biasout,abseout,rmseout,nflds,nmodels,lgood) implicit none ! Declare passed variables character(120),intent(in)::cpath character(80),intent(in)::filename integer(4),intent(in)::nflds integer(4),intent(in)::nmodels integer(4),intent(out)::nobsout(nflds,nmodels) real(4),intent(out)::biasout(nflds,nmodels) real(4),intent(out)::abseout(nflds,nmodels) real(4),intent(out)::rmseout(nflds,nmodels) logical,intent(out)::lgood ! Declare local parameters ! Declare local variables character(60) cmodel integer(4) nflds0,nobsmax0 integer(4),allocatable,dimension(:):: nobs real(4),allocatable,dimension(:):: bias,rmse,abse integer(4) imodel,n,nn logical lres1,lres2,lres3 !*********************************************************** !*********************************************************** allocate(nobs(nflds)) allocate(bias(nflds),rmse(nflds),abse(nflds)) print*,'in read_rawstats: nmodels=',nmodels lgood=.true. open (18,file=trim(cpath)//'/'//trim(filename),form='unformatted') do 100 imodel=1,nmodels read(18) cmodel read(18) nflds0,nobsmax0 print*,'in read_rawstats: imodel=',imodel print*,'in read_rawstats: cmodel=',trim(cmodel) print*,'in read_rawstats: nflds0,nobsmax0=',nflds0,nobsmax0 read(18) nobs,bias,abse,rmse lres1=any(bias(1:9) > 9998.) !any(bias(1:nflds) > 9998.) lres2=any(abse(1:9) > 9998.) !any(abse(1:nflds) > 9998.) lres3=any(rmse(1:9) > 9998.) !any(rmse(1:nflds) > 9998.) if (lres1 .or. lres2 .or. lres3) then print*,'in read_rawstats : start block' print*,'this run may have used an empty prepbufr file' print*,'set bias,abse,rmse to zero' bias=0. ; abse=0.; rmse=0. lgood=.false. print*,'in read_rawstats : end block' endif nobsout(:,imodel)=nobs(:) biasout(:,imodel)=bias(:) rmseout(:,imodel)=rmse(:) abseout(:,imodel)=abse(:) 100 continue close(18) deallocate(nobs) deallocate(bias,rmse,abse) return end !=============================================================== !=============================================================== subroutine read_rawstats_ccs_version(cpath,filename,nobsout, & biasout,abseout,rmseout,nflds,nmodels,lgood) implicit none ! Declare passed variables character(80),intent(in)::cpath,filename integer(4),intent(in)::nflds integer(4),intent(in)::nmodels integer(4),intent(out)::nobsout(nflds,nmodels) real(4),intent(out)::biasout(nflds,nmodels) real(4),intent(out)::abseout(nflds,nmodels) real(4),intent(out)::rmseout(nflds,nmodels) logical,intent(out)::lgood ! Declare local parameters integer(4),parameter::nobsmax=100000 ! Declare local variables character(60) cmodel integer(4) nflds0,nobsmax0 integer(4),allocatable,dimension(:):: nobs real(4),allocatable,dimension(:):: bias,rmse,abse real(4),allocatable::fno(:,:) real(4),allocatable::fn1(:,:) character(8),allocatable::cstation_all(:,:) character(8),allocatable::cprovider_all(:,:) character(8),allocatable::csubprovider_all(:,:) integer(4),allocatable:: itype_all(:,:) real(4),allocatable:: rlat_all(:,:) real(4),allocatable:: rlon_all(:,:) real(4),allocatable:: xx_all(:,:) real(4),allocatable:: yy_all(:,:) real(4),allocatable:: hgt0_all(:,:) real(4),allocatable:: hgt_all(:,:) real(4),allocatable:: slm_all(:,:) integer(4),allocatable::i1(:) real(4),allocatable::v0(:),v1(:) character(8),allocatable::c1(:) integer(4) imodel,n,nn logical lres1,lres2,lres3 !*********************************************************** !*********************************************************** allocate(nobs(nflds)) allocate(bias(nflds),rmse(nflds),abse(nflds)) allocate(fno(nobsmax,nflds)) allocate(fn1(nobsmax,nflds)) allocate(cstation_all(nobsmax,nflds)) allocate(cprovider_all(nobsmax,nflds)) allocate(csubprovider_all(nobsmax,nflds)) allocate(itype_all(nobsmax,nflds)) allocate(rlat_all(nobsmax,nflds)) allocate(rlon_all(nobsmax,nflds)) allocate(xx_all(nobsmax,nflds)) allocate(yy_all(nobsmax,nflds)) allocate(hgt0_all(nobsmax,nflds)) allocate(hgt_all(nobsmax,nflds)) allocate(slm_all(nobsmax,nflds)) print*,'in read_rawstats: nmodels=',nmodels lgood=.true. open (18,file=trim(cpath)//'/'//trim(filename),form='unformatted') do 100 imodel=1,nmodels read(18) cmodel read(18) nflds0,nobsmax0 print*,'in read_rawstats: imodel=',imodel print*,'in read_rawstats: cmodel=',trim(cmodel) print*,'in read_rawstats: nflds0,nobsmax0=',nflds0,nobsmax0 read(18) nobs,bias,abse,rmse lres1=any(bias(1:9) > 9998.) !any(bias(1:nflds) > 9998.) lres2=any(abse(1:9) > 9998.) !any(abse(1:nflds) > 9998.) lres3=any(rmse(1:9) > 9998.) !any(rmse(1:nflds) > 9998.) if (lres1 .or. lres2 .or. lres3) then print*,'in read_rawstats : start block' print*,'this run may have used an empty prepbufr file' print*,'set bias,abse,rmse to zero' bias=0. ; abse=0.; rmse=0. lgood=.false. print*,'in read_rawstats : end block' endif do nn=1,nflds n=nobs(nn) allocate(v0(n)) allocate(v1(n)) read(18) v0,v1 fno(1:n,nn)=v0(1:n) fn1(1:n,nn)=v1(1:n) deallocate(v0) deallocate(v1) enddo do nn=1,nflds n=nobs(nn) allocate(c1(n)) allocate(i1(n)) allocate(v1(n)) read(18) c1 ; cstation_all(1:n,nn)=c1(1:n) read(18) c1 ; cprovider_all(1:n,nn)=c1(1:n) read(18) c1 ; csubprovider_all(1:n,nn)=c1(1:n) read(18) i1 ; itype_all(1:n,nn)=i1(1:n) read(18) v1 ; rlat_all(1:n,nn)=v1(1:n) read(18) v1 ; rlon_all(1:n,nn)=v1(1:n) read(18) v1 ; xx_all(1:n,nn)=v1(1:n) read(18) v1 ; yy_all(1:n,nn)=v1(1:n) read(18) v1 ; hgt0_all(1:n,nn)=v1(1:n) read(18) v1 ; hgt_all(1:n,nn)=v1(1:n) read(18) v1 ; slm_all(1:n,nn)=v1(1:n) deallocate(c1) deallocate(i1) deallocate(v1) enddo !==> can call subroutine here to compute more specific ! stats. example, regional based, network based, etc. nobsout(:,imodel)=nobs(:) biasout(:,imodel)=bias(:) rmseout(:,imodel)=rmse(:) abseout(:,imodel)=abse(:) 100 continue deallocate(nobs) deallocate(bias,rmse,abse) deallocate(fno) deallocate(fn1) deallocate(cstation_all) deallocate(cprovider_all) deallocate(csubprovider_all) deallocate(itype_all) deallocate(rlat_all) deallocate(rlon_all) deallocate(xx_all) deallocate(yy_all) deallocate(hgt0_all) deallocate(hgt_all) deallocate(slm_all) return end !=============================================================== !===============================================================