program routinecv

!======================================================================
! compile as xlf90 -L/nwprod/lib -lw3_4
! Choose: a) number of cross-validation sets ncvs
!         b) total number of cycles ncycs. polulate
!            chh 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    !p,t,q,u,v,w,w2,wdir,td + gust,vis,pblh,dist

                                           !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 cycles run 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 sigf03,siganl
        real(4) biasin(nvars,3),absein(nvars,3),rmsein(nvars,3) !and sigcress
        character(80) cvar(nvars)

        character(2) chh(ncycs)         !store hours 00,01,02,03,04,05, ... 20,21,22,23
        character(8) cday(nt)
        character(80) cpath,cpathfull
        character(10) cycle
        character(6) cdaymonth
        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 chh,cday and cvar
                vbias2_hdvar,vrmse2_hdvar,vabse2_hdvar, &  !2dvar /per chh,cday and cvar
                vbias3_hdvar,vrmse3_hdvar,vabse3_hdvar     !cressman /per chh,cday and cvar

        real(4),dimension(ncycs,nvars):: & 
                vbias1_hvar,vrmse1_hvar,vabse1_hvar, &     !guess /per chh and cvar
                vbias2_hvar,vrmse2_hvar,vabse2_hvar, &     !2dvar /per chh and cvar
                vbias3_hvar,vrmse3_hvar,vabse3_hvar        !cressman /per chh and cvar

        real(4),dimension(nt,nvars):: & 
                vbias1_dvar,vrmse1_dvar,vabse1_dvar, &     !guess /per cday and cvar
                vbias2_dvar,vrmse2_dvar,vabse2_dvar, &     !2dvar /per cday and cvar
                vbias3_dvar,vrmse3_dvar,vabse3_dvar        !cressman /per cday 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) /'psfc_statistics'/
        data cvar(2) /'t2m_statistics'/
        data cvar(3) /'q2m_statistics'/
        data cvar(4) /'u10m_statistics'/
        data cvar(5) /'v10m_statistics'/
        data cvar(6) /'w10m_statistics'/
        data cvar(7) /'w210m_statistics'/
        data cvar(8) /'wdir_statistics'/
        data cvar(9) /'td2m_statistics'/
        data cvar(10) /'gust_statistics'/
        data cvar(11) /'vis_statistics'/
        data cvar(12) /'pblh_statistics'/
        data cvar(13) /'dist_statistics'/

        if (ncycs==24) ndcyc=1
        if (ncycs==8)  ndcyc=3
        do m=1,ncycs 
           write(chh(m),'(i2.2)') (m-1)*ndcyc
           print*,'m,chh(m)=',chh(m)
        enddo

        cpath='/ptmp/wx20pm/TEST/2010/'
!       cpath='/meso/noscrub/wx20pm/RETRIEVE_HPSS_HRES_FLDS/2010/tw/201006/'
!       cpath='/ptmp/wx20pm/newdump/'

        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(cday(n)(1:4),'(i4.4)') iyear
           write(cday(n)(5:6),'(i2.2)') imonth
           write(cday(n)(7:8),'(i2.2)') iday
           print*,'n,cday(n)=',n,trim(cday(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 
              cycle=cday(n)//chh(m)
              cdaymonth=cday(n)(1:6)
              cyear=cday(n)(1:4)
              do k=1,ncvs
                 write(ch2,'(i2.2)') k
!                cpathfull=trim(cpath)//cday(n)//'/'//cycle//'/cvset_'//ch2
!                cpathfull=trim(cpath)//cday(n)//'/'//cycle
                 cpathfull=trim(cpath)//cyear//'/'//cdaymonth//'/'//cday(n)//'/'//cycle
                 print*,'day,cycle,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 (sigf03,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=',chh(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=',chh(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=',chh(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, &
                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, &
                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

        integer(4) nlen
        character(3) cmonth0
        character(60) starttime
        character(4) cpos1neg9

        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)=trim(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 for surface pressure
        open (12,file='p_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,'(a20,f5.2,a1)')'p_bckgrmseavg="AVG =',vrmse1_hdvaravg(1), '"'
        write(12,'(a19,f5.2,a1)')'p_anlrmseavg="AVG =',vrmse2_hdvaravg(1), '"'
        write(12,'(a21,f5.2,a1)')'p_cressrmseavg="AVG =',vrmse3_hdvaravg(1), '"'

        write(12,'(a20,f5.2,a1)')'p_bckgbiasavg="AVG =',vbias1_hdvaravg(1), '"'
        write(12,'(a19,f5.2,a1)')'p_anlbiasavg="AVG =',vbias2_hdvaravg(1), '"'
        write(12,'(a21,f5.2,a1)')'p_cressbiasavg="AVG =',vbias3_hdvaravg(1), '"'

        write(12,'(a20,f5.2,a1)')'p_bckgabseavg="AVG =',vabse1_hdvaravg(1), '"'
        write(12,'(a19,f5.2,a1)')'p_anlabseavg="AVG =',vabse2_hdvaravg(1), '"'
        write(12,'(a21,f5.2,a1)')'p_cressabseavg="AVG =',vabse3_hdvaravg(1), '"'
        close(12)

!==> averages for 2m-T
        open (12,file='t_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,'(a20,f5.2,a1)')'t_bckgrmseavg="AVG =',vrmse1_hdvaravg(2), '"'
        write(12,'(a19,f5.2,a1)')'t_anlrmseavg="AVG =',vrmse2_hdvaravg(2), '"'
        write(12,'(a21,f5.2,a1)')'t_cressrmseavg="AVG =',vrmse3_hdvaravg(2), '"'

        write(12,'(a20,f5.2,a1)')'t_bckgbiasavg="AVG =',vbias1_hdvaravg(2), '"'
        write(12,'(a19,f5.2,a1)')'t_anlbiasavg="AVG =',vbias2_hdvaravg(2), '"'
        write(12,'(a21,f5.2,a1)')'t_cressbiasavg="AVG =',vbias3_hdvaravg(2), '"'

        write(12,'(a20,f5.2,a1)')'t_bckgabseavg="AVG =',vabse1_hdvaravg(2), '"'
        write(12,'(a19,f5.2,a1)')'t_anlabseavg="AVG =',vabse2_hdvaravg(2), '"'
        write(12,'(a21,f5.2,a1)')'t_cressabseavg="AVG =',vabse3_hdvaravg(2), '"'
        close(12)

!==> averages for 2m-q
        open (12,file='q_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,'(a20,f5.2,a1)')'q_bckgrmseavg="AVG =',vrmse1_hdvaravg(3), '"'
        write(12,'(a19,f5.2,a1)')'q_anlrmseavg="AVG =',vrmse2_hdvaravg(3), '"'
        write(12,'(a21,f5.2,a1)')'q_cressrmseavg="AVG =',vrmse3_hdvaravg(3), '"'

        write(12,'(a20,f5.2,a1)')'q_bckgbiasavg="AVG =',vbias1_hdvaravg(3), '"'
        write(12,'(a19,f5.2,a1)')'q_anlbiasavg="AVG =',vbias2_hdvaravg(3), '"'
        write(12,'(a21,f5.2,a1)')'q_cressbiasavg="AVG =',vbias3_hdvaravg(3), '"'

        write(12,'(a20,f5.2,a1)')'q_bckgabseavg="AVG =',vabse1_hdvaravg(3), '"'
        write(12,'(a19,f5.2,a1)')'q_anlabseavg="AVG =',vabse2_hdvaravg(3), '"'
        write(12,'(a21,f5.2,a1)')'q_cressabseavg="AVG =',vabse3_hdvaravg(3), '"'


        close(12)

!==> averages for 10m-u
        open (12,file='u_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,'(a20,f5.2,a1)')'u_bckgrmseavg="AVG =',vrmse1_hdvaravg(4), '"'
        write(12,'(a19,f5.2,a1)')'u_anlrmseavg="AVG =',vrmse2_hdvaravg(4), '"'
        write(12,'(a21,f5.2,a1)')'u_cressrmseavg="AVG =',vrmse3_hdvaravg(4), '"'

        write(12,'(a20,f5.2,a1)')'u_bckgbiasavg="AVG =',vbias1_hdvaravg(4), '"'
        write(12,'(a19,f5.2,a1)')'u_anlbiasavg="AVG =',vbias2_hdvaravg(4), '"'
        write(12,'(a21,f5.2,a1)')'u_cressbiasavg="AVG =',vbias3_hdvaravg(4), '"'

        write(12,'(a20,f5.2,a1)')'u_bckgabseavg="AVG =',vabse1_hdvaravg(4), '"'
        write(12,'(a19,f5.2,a1)')'u_anlabseavg="AVG =',vabse2_hdvaravg(4), '"'
        write(12,'(a21,f5.2,a1)')'u_cressabseavg="AVG =',vabse3_hdvaravg(4), '"'
        close(12)

!==> averages for 10m-v
        open (12,file='v_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,'(a20,f5.2,a1)')'v_bckgrmseavg="AVG =',vrmse1_hdvaravg(5), 'm/s"'
        write(12,'(a19,f5.2,a1)')'v_anlrmseavg="AVG =',vrmse2_hdvaravg(5), 'm/s"'
        write(12,'(a21,f5.2,a1)')'v_cressrmseavg="AVG =',vrmse3_hdvaravg(5), 'm/s"'

        write(12,'(a20,f5.2,a1)')'v_bckgbiasavg="AVG =',vbias1_hdvaravg(5), 'm/s"'
        write(12,'(a19,f5.2,a1)')'v_anlbiasavg="AVG =',vbias2_hdvaravg(5), 'm/s"'
        write(12,'(a21,f5.2,a1)')'v_cressbiasavg="AVG =',vbias3_hdvaravg(5), 'm/s"'

        write(12,'(a20,f5.2,a1)')'v_bckgabseavg="AVG =',vabse1_hdvaravg(5), 'm/s"'
        write(12,'(a19,f5.2,a1)')'v_anlabseavg="AVG =',vabse2_hdvaravg(5), 'm/s"'
        write(12,'(a21,f5.2,a1)')'v_cressabseavg="AVG =',vabse3_hdvaravg(5), 'm/s"'
        close(12)

!==> averages for 10m-w
        open (12,file='w_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,'(a20,f5.2,a1)')'w_bckgrmseavg="AVG =',vrmse1_hdvaravg(6), '"'
        write(12,'(a19,f5.2,a1)')'w_anlrmseavg="AVG =',vrmse2_hdvaravg(6), '"'
        write(12,'(a21,f5.2,a1)')'w_cressrmseavg="AVG =',vrmse3_hdvaravg(6), '"'

        write(12,'(a20,f5.2,a1)')'w_bckgbiasavg="AVG =',vbias1_hdvaravg(6), '"'
        write(12,'(a19,f5.2,a1)')'w_anlbiasavg="AVG =',vbias2_hdvaravg(6), '"'
        write(12,'(a21,f5.2,a1)')'w_cressbiasavg="AVG =',vbias3_hdvaravg(6), '"'

        write(12,'(a20,f5.2,a1)')'w_bckgabseavg="AVG =',vabse1_hdvaravg(6), '"'
        write(12,'(a19,f5.2,a1)')'w_anlabseavg="AVG =',vabse2_hdvaravg(6), '"'
        write(12,'(a21,f5.2,a1)')'w_cressabseavg="AVG =',vabse3_hdvaravg(6), '"'
        close(12)

!==> averages for 10m-w2
        open (12,file='w2_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,'(a21,f5.2,a1)')'w2_bckgrmseavg="AVG =',vrmse1_hdvaravg(7), '"'
        write(12,'(a20,f5.2,a1)')'w2_anlrmseavg="AVG =',vrmse2_hdvaravg(7), '"'
        write(12,'(a22,f5.2,a1)')'w2_cressrmseavg="AVG =',vrmse3_hdvaravg(7), '"'

        write(12,'(a21,f5.2,a1)')'w2_bckgbiasavg="AVG =',vbias1_hdvaravg(7), '"'
        write(12,'(a20,f5.2,a1)')'w2_anlbiasavg="AVG =',vbias2_hdvaravg(7), '"'
        write(12,'(a22,f5.2,a1)')'w2_cressbiasavg="AVG =',vbias3_hdvaravg(7), '"'

        write(12,'(a21,f5.2,a1)')'w2_bckgabseavg="AVG =',vabse1_hdvaravg(7), '"'
        write(12,'(a20,f5.2,a1)')'w2_anlabseavg="AVG =',vabse2_hdvaravg(7), '"'
        write(12,'(a22,f5.2,a1)')'w2_cressabseavg="AVG =',vabse3_hdvaravg(7), '"'
        close(12)

!==> averages for 10m-wd
        open (12,file='wd_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,'(a21,f5.2,a1)')'wd_bckgrmseavg="AVG =',vrmse1_hdvaravg(8), '"'
        write(12,'(a20,f5.2,a1)')'wd_anlrmseavg="AVG =',vrmse2_hdvaravg(8), '"'
        write(12,'(a22,f5.2,a1)')'wd_cressrmseavg="AVG =',vrmse3_hdvaravg(8), '"'

        write(12,'(a21,f5.2,a1)')'wd_bckgbiasavg="AVG =',vbias1_hdvaravg(8), '"'
        write(12,'(a20,f5.2,a1)')'wd_anlbiasavg="AVG =',vbias2_hdvaravg(8), '"'
        write(12,'(a22,f5.2,a1)')'wd_cressbiasavg="AVG =',vbias3_hdvaravg(8), '"'

        write(12,'(a21,f5.2,a1)')'wd_bckgabseavg="AVG =',vabse1_hdvaravg(8), '"'
        write(12,'(a20,f5.2,a1)')'wd_anlabseavg="AVG =',vabse2_hdvaravg(8), '"'
        write(12,'(a22,f5.2,a1)')'wd_cressabseavg="AVG =',vabse3_hdvaravg(8), '"'
        close(12)

!==> averages for 2m-td
        open (12,file='td_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,'(a21,f5.2,a1)')'td_bckgrmseavg="AVG =',vrmse1_hdvaravg(9), '"'
        write(12,'(a20,f5.2,a1)')'td_anlrmseavg="AVG =',vrmse2_hdvaravg(9), '"'
        write(12,'(a22,f5.2,a1)')'td_cressrmseavg="AVG =',vrmse3_hdvaravg(9), '"'

        write(12,'(a21,f5.2,a1)')'td_bckgbiasavg="AVG =',vbias1_hdvaravg(9), '"'
        write(12,'(a20,f5.2,a1)')'td_anlbiasavg="AVG =',vbias2_hdvaravg(9), '"'
        write(12,'(a22,f5.2,a1)')'td_cressbiasavg="AVG =',vbias3_hdvaravg(9), '"'

        write(12,'(a21,f5.2,a1)')'td_bckgabseavg="AVG =',vabse1_hdvaravg(9), '"'
        write(12,'(a20,f5.2,a1)')'td_anlabseavg="AVG =',vabse2_hdvaravg(9), '"'
        write(12,'(a22,f5.2,a1)')'td_cressabseavg="AVG =',vabse3_hdvaravg(9), '"'
        close(12)

!==> averages for 10m-gust
        open (12,file='gust_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,'(a23,f5.2,a1)')'gust_bckgrmseavg="AVG =',vrmse1_hdvaravg(10), '"'
        write(12,'(a22,f5.2,a1)')'gust_anlrmseavg="AVG =',vrmse2_hdvaravg(10), '"'
        write(12,'(a24,f5.2,a1)')'gust_cressrmseavg="AVG =',vrmse3_hdvaravg(10), '"'

        write(12,'(a23,f5.2,a1)')'gust_bckgbiasavg="AVG =',vbias1_hdvaravg(10), '"'
        write(12,'(a22,f5.2,a1)')'gust_anlbiasavg="AVG =',vbias2_hdvaravg(10), '"'
        write(12,'(a24,f5.2,a1)')'gust_cressbiasavg="AVG =',vbias3_hdvaravg(10), '"'

        write(12,'(a23,f5.2,a1)')'gust_bckgabseavg="AVG =',vabse1_hdvaravg(10), '"'
        write(12,'(a22,f5.2,a1)')'gust_anlabseavg="AVG =',vabse2_hdvaravg(10), '"'
        write(12,'(a24,f5.2,a1)')'gust_cressabseavg="AVG =',vabse3_hdvaravg(10), '"'
        close(12)

!==> averages for 10m-vis
        open (12,file='vis_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,f7.1,a1)')'vis_bckgrmseavg="AVG =',vrmse1_hdvaravg(11), '"'
        write(12,'(a21,f7.1,a1)')'vis_anlrmseavg="AVG =',vrmse2_hdvaravg(11), '"'
        write(12,'(a23,f7.1,a1)')'vis_cressrmseavg="AVG =',vrmse3_hdvaravg(11), '"'

        write(12,'(a22,f7.1,a1)')'vis_bckgbiasavg="AVG =',vbias1_hdvaravg(11), '"'
        write(12,'(a21,f7.1,a1)')'vis_anlbiasavg="AVG =',vbias2_hdvaravg(11), '"'
        write(12,'(a23,f7.1,a1)')'vis_cressbiasavg="AVG =',vbias3_hdvaravg(11), '"'

        write(12,'(a22,f7.1,a1)')'vis_bckgabseavg="AVG =',vabse1_hdvaravg(11), '"'
        write(12,'(a21,f7.1,a1)')'vis_anlabseavg="AVG =',vabse2_hdvaravg(11), '"'
        write(12,'(a23,f7.1,a1)')'vis_cressabseavg="AVG =',vabse3_hdvaravg(11), '"'
        close(12)

        return
        end
!===============================================================
!===============================================================
      subroutine read_rawstats(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
!===============================================================
!===============================================================