!*******************************************************************************
       subroutine bias_rmse_cv_mpi(cgrid,ista,iend,jsta,jend,nx,ny,igust,ivis,ipblh,idist, &
                               itotmodel,imodel,cmodel,rmusecv,rmuseb, &
                               lpadjust,mype,npe,fnamesuffix)


       
       use mpi

       use cressanl_common, only: nobsmax, xlocs, ylocs, hgt0s, &
                             hobs, rmuses, &
                             rtvts, cstations, &
                             obstypes, insubdoms, dups, bckgs, xberrs, &
                             ipointer,uvpointer,nmax,nflds, &
                             kps,kts,kqs,kus,kvs,kugrds,kvgrds,kws,kw2s,kwds, &
                             ktds,kgusts,kvis,kpblhs,kdists


       implicit none

!      note: consider mpi version
!
!      output
!      bias(1)  , rmse(1)  - psfc bias and rmse 
!      bias(2)  , rmse(2)  - 2m t bias and rmse of guess
!      bias(3)  , rmse(3)  - 2m q bias and rmse of guess
!      bias(4)  , rmse(4)  - 10m u bias and rmse of guess
!      bias(5)  , rmse(5)  - 10m v bias and rmse of guess
!      bias(6)  , rmse(6)  - 10m wind bias and rmse of guess
!      bias(7)  , rmse(7)  - 10m wind2 bias and rmse of guess
!      bias(8)  , rmse(8)  - 10m wdir bias and rmse of guess
!      bias(9)  , rmse(9)  - 2m td bias and rmse of guess
!      bias(10) , rmse(10) - 10m gust bias and rmse of guess
!      bias(11) , rmse(11) - vis bias and rmse of guess
!      bias(12) , rmse(12) - pblh bias and rmse of guess
!      bias(13) , rmse(13) - dist bias and rmse of guess
!
!      note: order for nflds is psfc,t,q,u,v,w,w2,wdir,td,gust,vis,pblh,dist
!------------------------------------------------------------------------------
! Declare passed variables
       integer(4),intent(in):: ista,iend,jsta,jend
       integer(4),intent(in):: nx,ny
       integer(4),intent(in):: igust,ivis,ipblh,idist
       integer(4),intent(in):: itotmodel
       integer(4),intent(in)::imodel
       integer(4),intent(in)::mype,npe
       real(4),intent(in)::rmusecv
       real(4),intent(in)::rmuseb
       logical,intent(in)::lpadjust
       character(60),intent(in)::cgrid
       character(60),intent(in)::cmodel
       character(*),intent(in)::fnamesuffix
       
! Declare local parameters
       integer(4),parameter::nflds0=13  !same as nflds
       integer(4),parameter::itotmodel0=5 !if itmodel > itmodel0 must increase itotmodel0
       real(4),parameter :: gravity=9.81
       real(4),parameter :: epsilon=1.e-03
       real(8),parameter :: r360_d=360._8
       real(4),parameter::spval=-999.

! Declare local variables
       integer(4),allocatable,dimension(:):: i1_info
       integer(4),allocatable,dimension(:):: i2_info
       integer(4),allocatable,dimension(:):: j1_info
       integer(4),allocatable,dimension(:):: j2_info
       integer(4),allocatable,dimension(:):: iaux1,iaux2,iaux3,iaux4

       integer(4) i1,i2,j1,j2
       integer(4) m1,m2,n1,n2

       integer(4) mype1

       integer(4) i,j,k,kk,n,nn,iu,iv
       integer(4) ierror,islab,nobsmax0

       integer(4),allocatable,dimension(:,:):: ipointer0
       real(4),allocatable,dimension(:):: xlocs0, ylocs0,hobs0,rmuses0,rtvts0,hgt00s
       logical,allocatable,dimension(:):: insubdoms0

       integer(4),    save:: allnobs(nflds0,3,itotmodel0)
       real(4),       save:: allbias(nflds0,3,itotmodel0) !the 3 accounts for rmuse=rmusecv.,1., and (rmusecv .or. 1.).
       real(4),       save:: allrmse(nflds0,3,itotmodel0)
       real(4),       save:: allabse(nflds0,3,itotmodel0)
       character(60), save:: allcmodel(itotmodel0)

       integer(4) nobs0(nflds)
       real(4) bias0(nflds)
       real(4) rmse0(nflds)
       real(4) abse0(nflds)
       character(60) cmodel0

       integer(4) nobs(nflds,3)
       real(4)    bias(nflds,3)
       real(4)    rmse(nflds,3)
       real(4)    abse(nflds,3)

       real(8)    s1(nflds,3)
       real(8)    s2(nflds,3)
       real(8)    s3(nflds,3)
       real(8)    s4(nflds,3)

       integer(4) nobs99(nflds,3)
       real(8)    s99(nflds,3)

       integer(4) rtime(6),nlon,nlat,nsig

       real(4),allocatable,dimension(:,:)::glon
       real(4),allocatable,dimension(:,:)::fis1,psfc1,t1,q1,u1,v1,w1,wd1,td1
       real(4),allocatable,dimension(:,:)::gust1,vis1,pblh1,dist1
       real(4),allocatable,dimension(:,:)::tv1
       real(4),allocatable,dimension(:,:)::ue,ve
       real(4),allocatable,dimension(:,:)::fis0
       real(4),allocatable,dimension(:,:)::tv0

       integer(4),allocatable,dimension(:,:):: iauxfield
       real(4),allocatable,dimension(:,:)::auxfield,auxfield2 
       real(4),allocatable,dimension(:,:,:)::field0

       character(8) cstation,cprovider,csubprovider,obstype
       integer(4) itype

       integer(4) npts
       integer(4) kobs0(13)
       real(4) amin0,amin,amax0,amax

       real(4) rlat,rlon,xx,yy,ob,ob_model,rmuse
       real(4) ds,rtvts9,zsges,tges
       real(4) fint
       real(4) uob,uob_model,vob,vob_model
       real(8) diff,diff2,du2,dv2
       real(8) zeroone(3)
       logical insubdom
       character(60) fname
       character(10) clun
       character(10) chvar

       logical fexist
       logical lstats1,lstats2,lstats3
       logical lanlerr
!*******************************************************************************

       if (itotmodel > itotmodel0) then 
          if (mype.eq.0) then 
               print*,'in bias_rmse_cv:'
               print*,'itotmodel,itotmodel0=',itotmodel,itotmodel0
               print*,' itotmodel must be <= itotmodel0. ... returning'
          endif
          return
       endif

       if (trim(cgrid)/='conus'  .and. trim(cgrid)/='alaska' .and. & 
           trim(cgrid)/='prico'  .and. trim(cgrid)/='hawaii' .and. & 
           trim(cgrid)/='cohres' .and. trim(cgrid)/='akhres' .and. &
           trim(cgrid)/='hrrr'   .and. trim(cgrid)/='guam'   .and. & 
           trim(cgrid)/='juneau' .and. trim(cgrid)/='cohresext' .and. &
           trim(cgrid)/='dtc' )  then

              if (mype.eq.0) print*,& 
                   'in bias_rmse_cv: unknown grid ', cgrid, 'returning ...'
         return
       endif

      if(mype.eq.0) then
       print*,'in bias_rmse_cv:: cgrid,nx,ny=',trim(cgrid),nx,ny
       print*,'in bias_rmse_cv:: cmodel,rmusecv,rmuseb=',trim(cmodel),rmusecv,rmuseb
       print*,'in bias_rmse_cv:: fnamesuffix=',trim(fnamesuffix)
      endif

      allocate(i1_info(npe))
      allocate(i2_info(npe))
      allocate(j1_info(npe))
      allocate(j2_info(npe))
      allocate(iaux1(npe))
      allocate(iaux2(npe))
      allocate(iaux3(npe))
      allocate(iaux4(npe))

      i1_info(:)=0      ;   iaux1(:)=0
      i2_info(:)=0      ;   iaux2(:)=0
      j1_info(:)=0      ;   iaux3(:)=0
      j2_info(:)=0      ;   iaux4(:)=0

      i1=max(1,ista-1)  ;   iaux1(mype+1)=i1
      i2=min(nx,iend+1) ;   iaux2(mype+1)=i2
      j1=max(1,jsta-1)  ;   iaux3(mype+1)=j1
      j2=min(ny,jend+1) ;   iaux4(mype+1)=j2


      call mpi_allreduce(iaux1,i1_info,npe,mpi_integer, &
                          mpi_sum,mpi_comm_world,ierror)

      call mpi_allreduce(iaux2,i2_info,npe,mpi_integer, &
                          mpi_sum,mpi_comm_world,ierror)

      call mpi_allreduce(iaux3,j1_info,npe,mpi_integer, &
                          mpi_sum,mpi_comm_world,ierror)

      call mpi_allreduce(iaux4,j2_info,npe,mpi_integer, &
                          mpi_sum,mpi_comm_world,ierror)


       allocate ( iauxfield (nx,ny) )
       allocate ( auxfield  (nx,ny) )
       allocate ( auxfield2 (nx,ny) )

       allocate ( glon  (i1:i2,j1:j2) )
       allocate ( psfc1 (i1:i2,j1:j2) )
       allocate ( fis1  (i1:i2,j1:j2) )
       allocate ( t1    (i1:i2,j1:j2) )
       allocate ( tv1   (i1:i2,j1:j2) )
       allocate ( q1    (i1:i2,j1:j2) )
       allocate ( u1    (i1:i2,j1:j2) ) ; allocate ( ue(i1:i2,j1:j2) )
       allocate ( v1    (i1:i2,j1:j2) ) ; allocate ( ve(i1:i2,j1:j2) )
       allocate ( w1    (i1:i2,j1:j2) )
       allocate ( wd1   (i1:i2,j1:j2) )
       allocate ( td1   (i1:i2,j1:j2) )
       allocate ( gust1 (i1:i2,j1:j2) )
       allocate ( vis1  (i1:i2,j1:j2) )
       allocate ( pblh1 (i1:i2,j1:j2) )
       allocate ( dist1 (i1:i2,j1:j2) )

       gust1=0.
       vis1=0.
       pblh1=0.
       dist1=0.

       lanlerr=trim(cmodel)=='errfield.dat'.or.trim(cmodel)=='errfield.dat_precval'

       open (52,file=trim(cmodel),form='unformatted')

       if (.not.lanlerr) then
          read(52) rtime,nlon,nlat,nsig
          read(52) auxfield,auxfield2    !glat,dx
          read(52) auxfield,auxfield2    !glon,dy

          if (mype.eq.0) print*,'in bias_rmse_cv nlon,nlat,nsig=',nlon,nlat,nsig

          glon (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)

          read(52) auxfield  ; psfc1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)
          read(52) auxfield  ; fis1  (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)
          read(52) auxfield  ; t1    (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)
          read(52) auxfield  ; q1    (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)
          read(52) auxfield  ; u1    (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)
          read(52) auxfield  ; v1    (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)


          if (igust>0 .or. ivis>0 .or. ipblh>0) then
             if (trim(cmodel)=='siganl' .or. trim(cmodel)=='sigfupdate02' .or. trim(cmodel)=='sigfupdate03') then
                do n=1,14               !must jump 14 real records to get to gust1
                   read(52) auxfield    !(see subroutine convert_binary_2d)
                enddo                   
              else
                do n=1,12               !must jump 12 real+integer records to get to gust1. this is for sigges
                   if (n==4 .or. n==5) then 
                     read(52) iauxfield 
                    else
                     read(52) auxfield
                   endif
                enddo                   
             endif

             if (igust > 0 ) read(52) auxfield  ;  gust1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)
             if (ivis  > 0 ) read(52) auxfield  ;  vis1  (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)
             if (ipblh > 0 ) read(52) auxfield  ;  pblh1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)
          endif

          do j=j1,j2
          do i=i1,i2
             tv1(i,j)=t1(i,j)+(1.+0.608*max(0.,q1(i,j)))
          enddo
          enddo

          call wspdwdir(cgrid,glon,u1,v1,ue,ve,w1,wd1,i1,i2,j1,j2,mype,npe)

          inquire(file='td.dat_'//trim(cmodel),exist=fexist) 
          if (fexist) then 
             open(57,file='td.dat_'//trim(cmodel),form='unformatted')
             read (57) auxfield  ; td1(i1:i2,j1:j2)=auxfield(i1:i2,j1:j2)
             close(57)
           else
             call get_dewpt(psfc1,q1,t1,td1,i1,i2,j1,j2)
          endif
       endif

       if (lanlerr) then
          read(52) auxfield  ; psfc1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)  !perr
          read(52) auxfield  ; t1    (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)  !terr
          read(52) auxfield  ; td1   (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)  !tderr
          read(52) auxfield  ; u1    (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)  !uerr
          read(52) auxfield  ; v1    (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)  !verr
          read(52) auxfield  ; q1    (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)  !qerr
          read(52) auxfield  ; wd1   (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)  !wdirerr2
          read(52) auxfield  ; w1    (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)  !wspderr
          read(52) auxfield  ; ue    (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)  !uerr2
          read(52) auxfield  ; ve    (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2)  !verr2
          read(52) auxfield                                                  !wdirerr

          if (igust > 0 ) read(52) auxfield  ; gust1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !gusterr
          if (ivis  > 0 ) read(52) auxfield  ; vis1  (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !viserr
          if (ipblh > 0 ) read(52) auxfield  ; pblh1 (i1:i2,j1:j2) = auxfield (i1:i2,j1:j2) !pblherr

          tv1(i1:i2,j1:j2)=t1(i1:i2,j1:j2)
       endif

       close(52)


       if (mype==0 ) print*,'*************in bias_rmse_cv********************'
       do nn=1,14

          if ( nn==11 .and. igust <= 0 ) cycle
          if ( nn==12 .and. ivis  <= 0 ) cycle
          if ( nn==13 .and. ipblh <= 0 ) cycle
          if ( nn==14 .and. idist <= 0 ) cycle

          if ( nn == 1  )  then ; amin0=minval(psfc1) ; amax0=maxval(psfc1) ; chvar='psfc1' ; endif
          if ( nn == 2  )  then ; amin0=minval(fis1)  ; amax0=maxval(fis1)  ; chvar='fis1'  ; endif
          if ( nn == 3  )  then ; amin0=minval(t1)    ; amax0=maxval(t1)    ; chvar='t1'    ; endif
          if ( nn == 4  )  then ; amin0=minval(tv1)   ; amax0=maxval(tv1)   ; chvar='tv1'   ; endif
          if ( nn == 5  )  then ; amin0=minval(td1)   ; amax0=maxval(td1)   ; chvar='td1'   ; endif
          if ( nn == 6  )  then ; amin0=minval(ue)    ; amax0=maxval(ue)    ; chvar='ue'    ; endif
          if ( nn == 7  )  then ; amin0=minval(ve)    ; amax0=maxval(ve)    ; chvar='ve'    ; endif
          if ( nn == 8  )  then ; amin0=minval(w1)    ; amax0=maxval(w1)    ; chvar='w1'    ; endif
          if ( nn == 9  )  then ; amin0=minval(wd1)   ; amax0=maxval(wd1)   ; chvar='wd1'   ; endif
          if ( nn == 10 )  then ; amin0=minval(q1)    ; amax0=maxval(q1)    ; chvar='q1'    ; endif
          if ( nn == 11 )  then ; amin0=minval(gust1) ; amax0=maxval(gust1) ; chvar='gust1' ; endif
          if ( nn == 12 )  then ; amin0=minval(vis1)  ; amax0=maxval(vis1)  ; chvar='vis1'  ; endif
          if ( nn == 13 )  then ; amin0=minval(pblh1) ; amax0=maxval(pblh1) ; chvar='pblh1' ; endif
          if ( nn == 14 )  then ; amin0=minval(dist1) ; amax0=maxval(dist1) ; chvar='dist1' ; endif

          call mpi_allreduce(amin0,amin,1,mpi_real,mpi_min,mpi_comm_world,ierror)
          call mpi_allreduce(amax0,amax,1,mpi_real,mpi_max,mpi_comm_world,ierror)

          if (mype==0) print*,trim(chvar),'min,max=',amin,amax
       enddo
       if (mype==0) print*,'************************************************'


       mype1=mype+1

       if (imodel==1) then
         allnobs(:,:,:)=9999
         allbias(:,:,:)=9999.000
         allrmse(:,:,:)=9999.000
         allabse(:,:,:)=9999.000
       endif

       nobs(:,:)=0
       s1(:,:)=0._8
       s2(:,:)=0._8
       s3(:,:)=0._8
       s4(:,:)=0._8

       fis1=fis1/gravity   !terrain field

       do 300 islab=1,npe
          m1=i1_info(islab)
          m2=i2_info(islab)
          n1=j1_info(islab)
          n2=j2_info(islab)
  
          allocate(field0 (m1:m2,n1:n2,nflds))
          allocate(fis0  (m1:m2,n1:n2))
          allocate(tv0   (m1:m2,n1:n2))
 
          if (mype1==islab) then
             do j=n1,n2
             do i=m1,m2
                field0 (i,j,1)  = psfc1(i,j)
                field0 (i,j,2)  = t1   (i,j)
                field0 (i,j,3)  = q1   (i,j)
                field0 (i,j,4)  = ue   (i,j)
                field0 (i,j,5)  = ve   (i,j)
                field0 (i,j,6)  = w1   (i,j)  !use same
                field0 (i,j,7)  = w1   (i,j)  !use same
                field0 (i,j,8)  = wd1  (i,j)
                field0 (i,j,9)  = td1  (i,j)
                field0 (i,j,10) = gust1(i,j)
                field0 (i,j,11) = vis1 (i,j)
                field0 (i,j,12) = pblh1(i,j)
                field0 (i,j,13) = dist1(i,j)
 
                fis0  (i,j)    = fis1(i,j)
                tv0   (i,j)    = tv1 (i,j)
             enddo
             enddo

             nobsmax0=nobsmax

             kobs0(1)=kps
             kobs0(2)=kts
             kobs0(3)=kqs
             kobs0(4)=kus
             kobs0(5)=kvs
             kobs0(6)=kws
             kobs0(7)=kw2s
             kobs0(8)=kwds
             kobs0(9)=ktds
             kobs0(10)=kgusts
             kobs0(11)=kvis
             kobs0(12)=kpblhs
             kobs0(13)=kdists
          endif

          npts=(n2-n1+1)*(m2-m1+1)*nflds
          call mpi_bcast (field0, npts,  mpi_real,  islab-1, mpi_comm_world, ierror)

          npts=(n2-n1+1)*(m2-m1+1)
          call mpi_bcast (fis0, npts,  mpi_real,  islab-1, mpi_comm_world, ierror)

          npts=(n2-n1+1)*(m2-m1+1)
          call mpi_bcast (tv0, npts,  mpi_real,  islab-1, mpi_comm_world, ierror)

          call mpi_bcast (nobsmax0, 1,  mpi_integer,  islab-1, mpi_comm_world, ierror)

          call mpi_bcast (kobs0, 13,  mpi_integer,  islab-1, mpi_comm_world, ierror)

          allocate(xlocs0(nobsmax0))
          allocate(ylocs0(nobsmax0))
          allocate(hobs0(nobsmax0))
          allocate(rmuses0(nobsmax0))
          allocate(rtvts0(nobsmax0))
          allocate(hgt00s(nobsmax0))
          allocate(insubdoms0(nobsmax0))
          allocate(ipointer0(nmax,nflds))

          if (mype1==islab) then
             xlocs0(:)=xlocs(:)
             ylocs0(:)=ylocs(:)
             hobs0(:)=hobs(:)
             rmuses0(:)=rmuses(:)
             rtvts0(:)=rtvts(:)
             hgt00s(:)=hgt0s(:)
             insubdoms0(:)=insubdoms(:)
             ipointer0(:,:)=ipointer(:,:)
          endif
          call mpi_bcast (xlocs0,  nobsmax0, mpi_real,  islab-1, mpi_comm_world, ierror)
          call mpi_bcast (ylocs0,  nobsmax0, mpi_real,  islab-1, mpi_comm_world, ierror)
          call mpi_bcast (hobs0,   nobsmax0, mpi_real,  islab-1, mpi_comm_world, ierror)
          call mpi_bcast (rmuses0, nobsmax0, mpi_real,  islab-1, mpi_comm_world, ierror)
          call mpi_bcast (rtvts0,  nobsmax0, mpi_real,  islab-1, mpi_comm_world, ierror)
          call mpi_bcast (hgt00s,   nobsmax0, mpi_real,  islab-1, mpi_comm_world, ierror)
          call mpi_bcast (insubdoms0, nobsmax0, mpi_logical,  islab-1, mpi_comm_world, ierror)

          npts=nmax*nflds
          call mpi_bcast (ipointer0, npts,  mpi_integer,  islab-1, mpi_comm_world, ierror)

          do 200 nn=1,nflds

             do 100 kk=1,kobs0(nn)
                if (mype==mod(kk-1,npe)) then
                   n=ipointer0(kk,nn)
                   if (.not.insubdoms0(n)) cycle                

                   xx=xlocs0(n) 
                   yy=ylocs0(n)
                   ob=hobs0(n)
                   rmuse=rmuses0(n)

                   call bilinear_2d0v2(field0(m1:m2,n1:n2,nn),m1,m2,n1,n2,fint,yy,xx)

                   if (nn==1) then  
                      if (lpadjust .and. .not.lanlerr) then 
                           call bilinear_2d0v2(fis0,m1,m2,n1,n2,zsges,yy,xx)

                           rtvts9=1. !assume this for now
                           if (rtvts9==0.) then
                               call bilinear_2d0v2(tv0,m1,m2,n1,n2,tges,yy,xx)
                             else
                               call bilinear_2d0v2(field0(m1:m2,n1:n2,2),m1,m2,n1,n2,tges,yy,xx)
                           endif

                           call adjust_pressure(zsges,tges,hgt00s(n),spval,spval,fint)
                      endif

                      fint=fint/100. !convert to hPa
                      ob=ob/100.     !convert to hPa
                   endif

                   if (nn==2) then
                       if(rtvts0(n)==0. .and. .not.lanlerr)  &
                          call bilinear_2d0v2(tv0,m1,m2,n1,n2,fint,yy,xx)
                   endif

                   if (nn==3) then
                       fint=fint*1000.     !use g/kg
                       ob=ob*1000.         !use g/kg
                   endif

                   if (lanlerr) then 
                         diff=dble(fint)
                      else
                         diff=dble(fint-ob)
   
                         if (n==8 .and. abs(diff) > 0._8) then 
                            diff2=min(abs(diff),abs(diff+r360_d),abs(diff-r360_d))
                            diff=diff2*diff/abs(diff)
                         endif
                   endif

                   lstats1=abs(rmuse-rmusecv)<=epsilon
                   lstats2=abs(rmuse-rmuseb)<=epsilon
                   lstats3=lstats1.or.lstats2

                   zeroone(:)=0._8

                   if (lstats1) then ; nobs(nn,1)=nobs(nn,1)+1 ; zeroone(1)=1._8 ; endif
                   if (lstats2) then ; nobs(nn,2)=nobs(nn,2)+1 ; zeroone(2)=1._8 ; endif
                   if (lstats3) then ; nobs(nn,3)=nobs(nn,3)+1 ; zeroone(3)=1._8 ; endif

                   do k=1,3
                      s1(nn,k)=s1(nn,k)+zeroone(k)*abs(diff)
                      s2(nn,k)=s2(nn,k)+zeroone(k)*diff
                      s3(nn,k)=s3(nn,k)+zeroone(k)*diff*diff
                   enddo

                   if (.not.lanlerr) then 
                      if (nn==7) then
                         call bilinear_2d0v2(field0(m1:m2,n1:n2,4),m1,m2,n1,n2,uob_model,yy,xx)
                         call bilinear_2d0v2(field0(m1:m2,n1:n2,5),m1,m2,n1,n2,vob_model,yy,xx) 
   
                         iu=ipointer0(kk,4)  ; uob=hobs0(iu)
                         iv=ipointer0(kk,5)  ; vob=hobs0(iv)

                         du2=dble(uob_model-uob)*dble(uob_model-uob) !(delta u squared)
                         dv2=dble(vob_model-vob)*dble(vob_model-vob) !(delta v squared)

                         do k=1,3
                            s4(nn,k)=s4(nn,k)+zeroone(k)*(du2+dv2)
                         enddo
                      endif
                   endif
                endif  !mype condition
100          continue
             call mpi_barrier(mpi_comm_world,ierror)

200       continue

          deallocate (field0)
          deallocate (fis0)
          deallocate (tv0)
          deallocate (xlocs0)
          deallocate (ylocs0)
          deallocate (hobs0)
          deallocate(rmuses0)
          deallocate(rtvts0)
          deallocate(hgt00s)
          deallocate(insubdoms0)
          deallocate (ipointer0)
300    continue

       npts=nflds*3
       nobs99=nobs ; call mpi_allreduce(nobs99, nobs, npts, mpi_integer, mpi_sum, mpi_comm_world, ierror)
       s99=s1      ; call mpi_allreduce(s99,    s1,   npts, mpi_real8,   mpi_sum, mpi_comm_world, ierror)
       s99=s2      ; call mpi_allreduce(s99,    s2,   npts, mpi_real8,   mpi_sum, mpi_comm_world, ierror)
       s99=s3      ; call mpi_allreduce(s99,    s3,   npts, mpi_real8,   mpi_sum, mpi_comm_world, ierror)
       s99=s4      ; call mpi_allreduce(s99,    s4,   npts, mpi_real8,   mpi_sum, mpi_comm_world, ierror)

       if (.not.lanlerr) s3(7,:)=s4(7,:)

       bias(:,:)=9999.000
       rmse(:,:)=9999.000
       abse(:,:)=9999.000

       do nn=1,nflds
          do k=1,3
             n=nobs(nn,k)
             if (n > 0) then
                abse(nn,k)=sngl(s1(nn,k)/dble(n))
                bias(nn,k)=sngl(s2(nn,k)/dble(n))
                rmse(nn,k)=sngl(sqrt(s3(nn,k)/dble(n)))
             endif
          enddo
       enddo

       abse(8,:)=20.              !force this for now. rethink bilinear interpolation for wdir
       bias(8,:)=10.              !force this for now. rethink bilinear interpolation for wdir
       rmse(8,:)=20.              !force this for now. rethink bilinear interpolation for wdir

       allcmodel(imodel)=cmodel
       allnobs(:,:,imodel)=nobs(:,:)
       allabse(:,:,imodel)=abse(:,:)
       allbias(:,:,imodel)=bias(:,:)
       allrmse(:,:,imodel)=rmse(:,:)

       if (imodel==itotmodel) then

          do 500 k=1,3 !accounts for rmuse=rmusecv.,1., and (rmusecv. & 1.).
             if (k==1) then 
                if (rmusecv==-2.)  fname='for_rmuse_neg2'
                if (rmusecv==-3.)  fname='for_rmuse_neg3'
             endif
             if (k==2)  fname='for_rmuse_pos1'
             if (k==3)  fname='for_rmuse_neg3pos1'

             fname=trim(fname)//trim(fnamesuffix)

             if (mype==0) then
                open (17,file='stats_'//trim(fname),form='formatted')
                open (18,file='rawstats.dat_'//trim(fname),form='unformatted')

                print*,'in bias_rmse_cv: k,fname=',k,trim(fname)

                do i=1,itotmodel
                   cmodel0=allcmodel(i)
                   nobsmax0=maxval(allnobs(:,k,i))  !no need really. but leave to conform with content of older files
                   nobs0(:)=allnobs(:,k,i)
                   bias0(:)=allbias(:,k,i)
                   abse0(:)=allabse(:,k,i)
                   rmse0(:)=allrmse(:,k,i)

                   print*,'           '
                   print*,'verification statistics for ',trim(cmodel0)
                   write(17,*) trim(cmodel0)
                   do n=1,nflds
                      print 123,'n,nobs,abse,bias,rmse=',n,nobs0(n),abse0(n),bias0(n),rmse0(n)
                      write(17,123) 'n,nobs,abse,bias,rmse=',n,nobs0(n),abse0(n),bias0(n),rmse0(n)
                   enddo

                   write(18) cmodel0
                   write(18) nflds,nobsmax0
                   write(18) nobs0,bias0,abse0,rmse0
                enddo
             endif
             call mpi_barrier(mpi_comm_world,ierror)
500       continue
          close(17)
          close(18)
       endif
123    format(a23,i3,i7,3(f14.6))

       deallocate(i1_info)
       deallocate(i2_info)
       deallocate(j1_info)
       deallocate(j2_info)
       deallocate(iaux1)
       deallocate(iaux2)
       deallocate(iaux3)
       deallocate(iaux4)

       deallocate(glon)
       deallocate(psfc1)
       deallocate(fis1)
       deallocate(t1)
       deallocate(tv1)
       deallocate(q1)
       deallocate(u1)  ; deallocate(ue)
       deallocate(v1)  ; deallocate(ve)
       deallocate(w1)
       deallocate(wd1)
       deallocate(td1)
       deallocate(gust1)
       deallocate(vis1)
       deallocate(pblh1)
       deallocate(dist1)
       deallocate(iauxfield)
       deallocate(auxfield)
       deallocate(auxfield2)
end subroutine bias_rmse_cv_mpi
!------------------------------------------------------
    subroutine adjust_pressure(zsges,tges0,hgt0,tob,spval,fint)

    implicit none

!Declare passed variables
    real(4),intent(in):: zsges
    real(4),intent(in):: tges0
    real(4),intent(in):: hgt0
    real(4),intent(in):: tob
    real(4),intent(in):: spval
    real(4),intent(inout):: fint

!Declare local variables
    real(4) rdp,rdelz,pges,tges
    real(4) half,half_tlapse,g_over_rd

    half=0.5
    half_tlapse=0.00325
    g_over_rd=9.81/287.04

    pges=fint

    rdp=0.

    rdelz=hgt0-zsges
    tges=tges0

    if (abs(tob-spval).gt.1.e-03) then
       tges = half*(tges0+tob)
      else
       if(rdelz < 0.)then
         tges=tges0-half_tlapse*rdelz
       endif
     endif

     rdp = g_over_rd*rdelz/tges

     !Adjust hydrostatically
     fint=1000.*exp(log(pges/1000.) - rdp) !/100. !hPa

end subroutine adjust_pressure
!------------------------------------------------------
!------------------------------------------------------
