 program anl_error
  
!********************************************************************
! abstract: use barnes or cressman analysis to obtain gridded field *
!           of analysis uncertainty                                 *
!                                                                   *
! program history log:                                              *
!   2005-10-08  pondeca                                             *
!                                                                   *
!********************************************************************
  use mpi
  use kinds, only: i_kind


  implicit none

! Declare local variables
  integer(i_kind) mype,npe,ierror
  logical projway
  logical anlqlty
  logical lobjanl
  logical climatology
  logical cvbasedrecalibration
  character(60) cvbasedcmodel0
  integer(4) miter0

  namelist/anlerrmethod/projway,anlqlty,lobjanl,climatology,cvbasedrecalibration, & 
                        cvbasedcmodel0,miter0

! MPI setup
  call mpi_init(ierror)
  call mpi_comm_size(mpi_comm_world,npe,ierror)
  call mpi_comm_rank(mpi_comm_world,mype,ierror)

  call run_biascor(mype,npe)
   
  call mpi_barrier(mpi_comm_world,ierror)

  projway=.true.
  climatology=.false.
  anlqlty=.false.
  lobjanl=.false.
  cvbasedrecalibration=.false.
  cvbasedcmodel0='siganl'
  miter0=2  !default number of outer loops in gsi 
   
  open (55,file='lanczosparm.anl',form='formatted')
  read(55,anlerrmethod)
  close(55)

  projway=.true.  !must force this after reading namelist, since
                  !projway=.false. is no longer a valid option /18Oct2012

  if (mype==0) print*,'projway=',projway
  if (mype==0) print*,'anlqlty=',anlqlty
  if (mype==0) print*,'lobjanl=',lobjanl
  if (mype==0) print*,'climatology=',climatology
  if (mype==0) print*,'cvbasedrecalibration=',cvbasedrecalibration
  if (mype==0) print*,'cvbasedcmodel0=',cvbasedcmodel0
  if (mype==0) print*,'miter0=',miter0

  if (projway) call lanczos_anlerr(anlqlty,lobjanl,climatology,cvbasedrecalibration, & 
                                   cvbasedcmodel0,miter0,mype,npe)
! if (.not.projway) call simplified_anlerr(mype,npe)

  call mpi_finalize(ierror)

  end
!=================================================================================
!=================================================================================
 subroutine lanczos_anlerr(anlqlty,lobjanl,climatology,cvbasedrecalibration, & 
                           cvbasedcmodel0,miter0,mype,npe)

!********************************************************************
! abstract: use barnes or cressman analysis to obtain gridded field *
!           of analysis uncertainty                                 *
! Must run with a minimum of 5 tasks                                *
! Change this later so that less tasks also possible                *
!                                                                   *
! program history log:                                              *
!   2006-11-20  pondeca                                             *
!                                                                   *
!********************************************************************
  use mpi
  use kinds, only: r_single,r_kind,r_double,i_kind
  use constants, only: zero,one,tiny_r_kind,tiny_single,half

  use cressanl_common, only: nobsmax, xlocs, ylocs, hgt0s, hgts, &
                             hobs, rmuses, oberrs, &
                             dtimes, cstations, & 
                             obstypes, jpointer, & 
                             nflds,kflds, & 
                             kps,kts,kqs,kus,kvs,kugrds,kvgrds,kws,kw2s,kwds, & 
                             ktds,kgusts,kvis,kpblhs,kdists, & 
                             create_cressanl_common, & 
                             load_cressanl_common, & 
                             destroy_cressanl_common                             

  use errs_common, only: psierr,chierr,uerr,verr,uerr2,verr2, & 
                         wspderr,wdirerr,wdirerr2,terr,tderr, & 
                         qerr,perr,gusterr,viserr,pblherr,disterr
  use errs_common, only: create_errs_common,destroy_errs_common

  implicit none

! Declare passed variables
  integer(i_kind),intent(in):: mype,npe
  logical,intent(in):: anlqlty
  logical,intent(in):: lobjanl
  logical,intent(in):: climatology
  logical,intent(in):: cvbasedrecalibration
  character(60),intent(in):: cvbasedcmodel0
  integer(4),intent(in):: miter0

! Declare local variables
  character(60) cgrid
  integer(i_kind) ierror,ierror2
  integer(i_kind) ista,iend,jsta,jend
  integer(i_kind) nlon,nlat,nsig,jpch,npred,jtype,npredp,jiter,nv,nvused
  integer(i_kind) nlon1,nlat1,nsig1,jpch1,npred1,jtype1,npredp1,jiter1,nv1
  integer(i_kind) nrf3_sf,nrf3_vp,nrf3_t,nrf3_q,nrf3_oz,nrf3_cw
  integer(i_kind) nrf2_ps,nrf2_sst,nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist
  integer(i_kind) nrf1_rad,nrf1_prec
  integer(i_kind) i,j,k,kall,lun1,lun2,m,n,ifield,nsym,kk
  integer(i_kind) num_slabs,num_pad,islab,islab_prev
  integer(i_kind) isub,jsub,ijsub
  integer(i_kind) nx,ny,nk
  integer(i_kind) npointer(12) !total number of variables
  integer(i_kind),allocatable,dimension(:)::iscnt,isdisp,ircnt,irdisp
  integer(i_kind),allocatable,dimension(:)::ijpts_in_sub
  integer(i_kind),allocatable,dimension(:)::ista_info,iend_info,jsta_info,jend_info
  integer(i_kind),allocatable,dimension(:)::iaux0,iaux1,iaux2,iaux3,iaux4
  real(r_double) innerprod,dotp0,dotp
  real(r_kind),allocatable,dimension(:,:)::h_slab,g_slab
  real(r_kind),allocatable,dimension(:)::h_rad,h_prec,g_rad,g_prec
  real(r_kind),allocatable,dimension(:)::hfield,gfield
  real(r_kind),allocatable,dimension(:,:,:,:)::hall_loc,gall_loc
  real(r_kind),allocatable,dimension(:,:,:,:)::hall_loc2,gall_loc2
  real(r_double),allocatable,dimension(:)::cnorm

  character(2) clun1
  character(3) clun2
  character(3) clun3
  character(20) varname1,varname2
  logical twodvar_regional
  logical fexist1
  logical,parameter:: leigenvector_write=.false.

  real(r_kind),allocatable,dimension(:,:,:)::psi
  real(r_kind),allocatable,dimension(:,:,:)::chi
  real(r_kind),allocatable,dimension(:,:,:)::t
  real(r_kind),allocatable,dimension(:,:,:)::q
  real(r_kind),allocatable,dimension(:,:,:)::oz
  real(r_kind),allocatable,dimension(:,:,:)::cwmr
  real(r_kind),allocatable,dimension(:,:)::sfcp
  real(r_kind),allocatable,dimension(:,:)::sfct
  real(r_kind),allocatable,dimension(:,:)::gust,vis,pblh,dist
  real(r_kind),allocatable,dimension(:,:,:,:)::u,u2
  real(r_kind),allocatable,dimension(:,:,:,:)::v,v2
  real(r_kind),allocatable,dimension(:,:,:,:)::wdir,wdir2
  real(r_kind),allocatable,dimension(:,:,:)::usum,usum2
  real(r_kind),allocatable,dimension(:,:,:)::vsum,vsum2
  real(r_kind),allocatable,dimension(:,:,:)::wdirsum,wdirsum2
  real(r_kind),allocatable,dimension(:,:)::uaux1
  real(r_kind),allocatable,dimension(:,:)::vaux1
  real(r_kind),allocatable,dimension(:,:)::uaux2
  real(r_kind),allocatable,dimension(:,:)::vaux2
  real(r_kind),allocatable,dimension(:,:)::coeffx
  real(r_kind),allocatable,dimension(:,:)::coeffy

  real(r_kind),allocatable,dimension(:,:)::grad,hrad
  real(r_kind),allocatable,dimension(:,:)::gprec,hprec
  real(r_kind),allocatable,dimension(:,:)::grad2,hrad2
  real(r_kind),allocatable,dimension(:,:)::gprec2,hprec2
  real(r_double),allocatable,dimension(:,:)::hgdotaux,hgdot
  real(r_kind),allocatable,dimension(:,:)::tmatrix
  real(r_kind),allocatable,dimension(:,:)::slab1,slab2
  real(r_kind),allocatable,dimension(:)::alpha,gamma
  real(r_kind),allocatable,dimension(:)::beta,delta,theta
  real(r_kind),allocatable,dimension(:)::atest,eval
  real(r_kind),allocatable,dimension(:,:)::rtest
  real(r_kind),allocatable,dimension(:,:)::tempb8,tempc8
  real(r_single),allocatable,dimension(:,:)::tempa,tempb
  real(r_single),allocatable,dimension(:,:)::glon
  real(r_single),allocatable,dimension(:,:):: climerr 
                 
  real(r_single) ds0
  real(r_kind),parameter::ubar=10.
  real(r_kind),parameter::vbar=10.
  real(r_kind),parameter::r25=25._r_kind
  real(r_kind),parameter::r210=210._r_kind
  real(r_kind),parameter::r265=265._r_kind
  real(r_kind),parameter::r38_5=38.5_r_kind
  real(r_kind),parameter::r262_5=262.5_r_kind
  real(r_kind),parameter::r225=225._r_kind
  real(r_kind) dg2rad,xn,elonv,angle2,sinx2,cosx2
  real(r_kind) alpha0,gamma0
  integer(i_kind),parameter::nvmin=4

  real(r_single) pbiascor,tbiascor,qbiascor,ubiascor,vbiascor,tdbiascor, &
          gustbiascor,visbiascor,pblhbiascor
  logical lpbiascor,ltbiascor,lqbiascor,lubiascor,lvbiascor,ltdbiascor ,& 
          lgustbiascor,lvisbiascor,lpblhbiascor
  logical lbiascor
  logical mkrjlists
  logical hwrfblend

  integer(i_kind) icolor, ikey, mpi_comm_new, mypenew, npenew
  integer(i_kind) mypecutoff

  real(r_single),parameter::ribuffer_km=250.
  real(r_single),parameter::rjbuffer_km=250.

  character(60) tdfname
  character(60) filename

  integer(i_kind) ista2,iend2,jsta2,jend2   !memory dimensions
  real(8) amin8,amax8,amin82,amax82

  namelist/gridname/cgrid,lbiascor,lpbiascor,pbiascor,ltbiascor,tbiascor, &
                    lqbiascor,qbiascor,lubiascor,ubiascor,lvbiascor,vbiascor, &
                    ltdbiascor,tdbiascor,lgustbiascor,gustbiascor,lvisbiascor, &
                    visbiascor,lpblhbiascor,pblhbiascor,mkrjlists,hwrfblend
!*************************************************************************
!*************************************************************************
  data cgrid/'conus'/

  open (55,file='gridname_input',form='formatted')
  read(55,gridname)
  close(55)

  if (mype==0) print*,'in lanczos_anlerr:, cgrid=',cgrid

  if (trim(cgrid) /= 'conus' .and. &
      trim(cgrid) /= 'alaska' .and. &
      trim(cgrid) /= 'hawaii' .and. &
      trim(cgrid) /= 'guam' .and. &
      trim(cgrid) /= 'prico' .and. &
      trim(cgrid) /= 'cohres' .and. &
      trim(cgrid) /= 'akhres' .and. &
      trim(cgrid) /= 'hrrr' .and. &
      trim(cgrid) /= 'dtc' .and. &
      trim(cgrid) /= 'juneau' .and. &
      trim(cgrid) /= 'cohresext') then
      print*,'in anl_error: unknown grid ',trim(cgrid),'...aborting'
      call mpi_abort(mpi_comm_world,ierror2,ierror)
  endif

  call domain_dims(cgrid,nx,ny,ds0)
  if (mype==0) print*,'in lanczos_anlerr:, ds0=',ds0,nx,ny

  lun1=19
  lun2=29

  if (mype.eq.0) then
    open (lun1,file='grady.dat',form='unformatted')
    read (lun1) nlon,nlat,nsig,jpch,npred,jtype,npredp,jiter,nv,alpha0,gamma0, & 
                nrf3_sf,nrf3_vp,nrf3_t,nrf3_q,nrf3_oz,nrf3_cw, &
                nrf2_ps,nrf2_sst,nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist

    print*,' -----------in lanczos_anlerr-----------------'

    print*,'nlon,nlat,nsig,jpch,npred,jtype,npredp,jiter,nv=', &
            nlon,nlat,nsig,jpch,npred,jtype,npredp,jiter,nv

    print*,'nrf3_sf,nrf3_vp,nrf3_t,nrf3_q,nrf3_oz,nrf3_cw=', & 
            nrf3_sf,nrf3_vp,nrf3_t,nrf3_q,nrf3_oz,nrf3_cw

    print*,'nrf2_ps,nrf2_sst,nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist=',&
            nrf2_ps,nrf2_sst,nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist


    nrf1_rad=-1 
    if ( jpch*npred  > 0 ) nrf1_rad=10  !any positive integer

    nrf1_prec=-1
    if (jtype*npredp > 0 ) nrf1_prec=11 !any positive integer
    close(lun1) 

    !may have used less iterations than original 'nv=niter(1)'. address that possibility:

    open (lun1,file='used_iterations.dat',form='formatted')
    read (lun1,*) nvused
    nvused=nvused-1
    close(lun1) 
    print*,'nv,nvused,nvmin=',nv,nvused,nvmin
    nv=min(nvused,30)                               !added 19Nov2007
  endif


  call mp_flush(1)
  call mpi_bcast ( nlon      , 1 , mpi_integer , 0 , mpi_comm_world,ierror)
  call mpi_bcast ( nlat      , 1 , mpi_integer , 0 , mpi_comm_world,ierror)
  call mpi_bcast ( nsig      , 1 , mpi_integer , 0 , mpi_comm_world,ierror)
  call mpi_bcast ( jpch      , 1 , mpi_integer , 0 , mpi_comm_world,ierror)
  call mpi_bcast ( npred     , 1 , mpi_integer , 0 , mpi_comm_world,ierror)
  call mpi_bcast ( jtype     , 1 , mpi_integer , 0 , mpi_comm_world,ierror)
  call mpi_bcast ( npredp    , 1 , mpi_integer , 0 , mpi_comm_world,ierror)
  call mpi_bcast ( jiter     , 1 , mpi_integer , 0 , mpi_comm_world,ierror)
  call mpi_bcast ( nv        , 1 , mpi_integer , 0 , mpi_comm_world,ierror)

  call mpi_bcast ( nrf3_sf   , 1 , mpi_integer , 0 , mpi_comm_world,ierror )
  call mpi_bcast ( nrf3_vp   , 1 , mpi_integer , 0 , mpi_comm_world,ierror )
  call mpi_bcast ( nrf3_t    , 1 , mpi_integer , 0 , mpi_comm_world,ierror )
  call mpi_bcast ( nrf3_q    , 1 , mpi_integer , 0 , mpi_comm_world,ierror )
  call mpi_bcast ( nrf3_oz   , 1 , mpi_integer , 0 , mpi_comm_world,ierror )
  call mpi_bcast ( nrf3_cw   , 1 , mpi_integer , 0 , mpi_comm_world,ierror )
  call mpi_bcast ( nrf2_ps   , 1 , mpi_integer , 0 , mpi_comm_world,ierror )
  call mpi_bcast ( nrf2_sst  , 1 , mpi_integer , 0 , mpi_comm_world,ierror )
  call mpi_bcast ( nrf2_gust , 1 , mpi_integer , 0 , mpi_comm_world,ierror )
  call mpi_bcast ( nrf2_vis  , 1 , mpi_integer , 0 , mpi_comm_world,ierror )
  call mpi_bcast ( nrf2_pblh , 1 , mpi_integer , 0 , mpi_comm_world,ierror )
  call mpi_bcast ( nrf2_dist , 1 , mpi_integer , 0 , mpi_comm_world,ierror )
  call mpi_bcast ( nrf1_rad  , 1 , mpi_integer , 0 , mpi_comm_world,ierror )
  call mpi_bcast ( nrf1_prec , 1 , mpi_integer , 0 , mpi_comm_world,ierror )
!***********************************************************
!***********************************************************
!==>Divide up full horizontal domain into subdomains

  call horiz_domain_partition(nlon,nlat,mype,npe, & 
                                   ista,iend,jsta,jend)
   
  print*,'ista,iend,jsta,jend=',ista,iend,jsta,jend
!***********************************************************
!***********************************************************
  if (anlqlty .or. cvbasedrecalibration)  then
     call get_streamlined_stats_cv_mpi(cgrid,ista,iend,jsta,jend,ribuffer_km,mype,npe)

     tdfname='sigges' ;        call td_flds(nlon,nlat,tdfname,mype,npe)
     tdfname='siganl' ;        call td_flds(nlon,nlat,tdfname,mype,npe)

     if (cvbasedrecalibration) call td_flds(nlon,nlat,cvbasedcmodel0,mype,npe)

     call create_cressanl_common(ista,iend,jsta,jend,nlon,nlat,npe)
     call load_cressanl_common(mype)

     call mpi_barrier(mpi_comm_world,ierror)
  endif
!***********************************************************
!***********************************************************
  if (climatology .or. (nv < nvmin)) goto 3500
 
!***********************************************************
!***********************************************************
!==>Let every processor know about the number of horizontal 
! points and the starting and ending i and j values on each 
! subdomain

  allocate(ijpts_in_sub(npe))
  allocate(ista_info(npe))
  allocate(iend_info(npe))
  allocate(jsta_info(npe))
  allocate(jend_info(npe))
  allocate(iaux0(npe))
  allocate(iaux1(npe))
  allocate(iaux2(npe))
  allocate(iaux3(npe))
  allocate(iaux4(npe))


  ijpts_in_sub(:)=0   ;   iaux0(:)=0
  ista_info(:)=0      ;   iaux1(:)=0
  iend_info(:)=0      ;   iaux2(:)=0
  jsta_info(:)=0      ;   iaux3(:)=0
  jend_info(:)=0      ;   iaux4(:)=0
 
  iaux0(mype+1)=(iend-ista+1)*(jend-jsta+1)
  iaux1(mype+1)=ista
  iaux2(mype+1)=iend
  iaux3(mype+1)=jsta
  iaux4(mype+1)=jend

  call mpi_allreduce(iaux0,ijpts_in_sub,npe,mpi_integer, &
                      mpi_sum,mpi_comm_world,ierror)

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

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

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

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

  if (mype.eq.0) then
   do i=1,npe
      print*,'i,ijpts_in_sub=',i,ijpts_in_sub(i)
      print*,'i,ista_info,iend_info=',i,ista_info(i),iend_info(i)
      print*,'i,jsta_info,jend_info=',i,jsta_info(i),jend_info(i)
   enddo
  endif

!***********************************************************
!==>Total number of horizontal slabs in one gradient file
!excluding gradient with respect to rad or prec

  num_slabs=0
  if (nrf3_sf   > 0 )  num_slabs = num_slabs + nsig
  if (nrf3_vp   > 0 )  num_slabs = num_slabs + nsig
  if (nrf3_t    > 0 )  num_slabs = num_slabs + nsig
  if (nrf3_q    > 0 )  num_slabs = num_slabs + nsig
  if (nrf3_oz   > 0 )  num_slabs = num_slabs + nsig
  if (nrf3_cw   > 0 )  num_slabs = num_slabs + nsig

  if (nrf2_ps   > 0 )  num_slabs = num_slabs + 1
  if (nrf2_sst  > 0 )  num_slabs = num_slabs + 1
  if (nrf2_gust > 0 )  num_slabs = num_slabs + 1
  if (nrf2_vis  > 0 )  num_slabs = num_slabs + 1
  if (nrf2_pblh > 0 )  num_slabs = num_slabs + 1
  if (nrf2_dist > 0 )  num_slabs = num_slabs + 1

  if (mype .eq. 0) &  
  print*,'# of horizontal slabs in gradient file=',num_slabs
!***********************************************************
!==>prepare for alltoallv

  if (mod(num_slabs,npe) == 0) then 
     num_pad=num_slabs
   else
     num_pad=(num_slabs/npe +1)*npe
  endif

  if (mype .eq. 0) print*,'num_pad=',num_pad

  allocate(hall_loc(ista:iend,jsta:jend,num_pad,0:nv))
  allocate(gall_loc(ista:iend,jsta:jend,num_pad,0:nv))

  if (nrf1_rad > 0)  then 
     allocate(hrad(jpch*npred,0:nv))
     allocate(grad(jpch*npred,0:nv))
  endif

  if (nrf1_prec > 0) then
     allocate(hprec(jtype*npredp,0:nv))
     allocate(gprec(jtype*npredp,0:nv))
  endif

  allocate(iscnt(npe))
  allocate(isdisp(npe))
  allocate(ircnt(npe))
  allocate(irdisp(npe))

  do i=1,npe
   iscnt(i)=ijpts_in_sub(i)
  enddo

  isdisp(1)=0
  do i=1,npe
   if (i /=1 ) isdisp(i)=isdisp(i-1)+iscnt(i-1)
  enddo

  do i=1,npe
   ircnt(i)=ijpts_in_sub(mype+1)
   irdisp(i)=(i-1)*ijpts_in_sub(mype+1)
  enddo

!***********************************************************
!***********************************************************
!==>for each grad field, have all processors read in the 
!   data slabs. evoke alltoallv to send fields to subdomains 
!   everytime a group of npe slabs have been read in or when 
!   end of slab count has been reached

  allocate(h_slab(nlon,nlat))
  allocate(hfield(nlon*nlat))
  allocate(h_rad(jpch*npred))
  allocate(h_prec(jtype*npredp))

  allocate(g_slab(nlon,nlat))
  allocate(gfield(nlon*nlat))
  allocate(g_rad(jpch*npred))
  allocate(g_prec(jtype*npredp))

  allocate(alpha(0:nv))
  allocate(gamma(-1:nv-1))

  write (clun1(1:2),'(i2.2)') jiter

  do 400 n=0,nv

!   print*,'gradient field number=',n

    write (clun2(1:3),'(i3.3)') n

    open (lun1,file='grady.dat_'//clun1//'_'//clun2,form='unformatted')
    open (lun2,file='gradx.dat_'//clun1//'_'//clun2,form='unformatted')

    rewind (lun1)
    rewind (lun2)

    read (lun1) nlon1,nlat1,nsig1,jpch1,npred1, & 
                jtype1,npredp1,jiter1,nv1,alpha(n),gamma(n-1)

    read (lun2) nlon1,nlat1,nsig1,jpch1,npred1, & 
                jtype1,npredp1,jiter1,nv1,alpha(n),gamma(n-1)

!   print*,'nlon1,nlat1,nsig1,jpch1,npred1,jtype1,npredp1,jiter1,nv1=',&
!          nlon1,nlat1,nsig1,jpch1,npred1,jtype1,npredp1,jiter1,nv1

    islab_prev=1
    do 200 islab=1,num_slabs

       if(mype==mod(islab-1,npe)) then
         read(lun1) h_slab
         call vectorform(h_slab,hfield,nlon,nlat, & 
                    ista_info,iend_info,jsta_info,jend_info,npe)

         read(lun2) g_slab
         call vectorform(g_slab,gfield,nlon,nlat, & 
                    ista_info,iend_info,jsta_info,jend_info,npe)

!        print*,'min,max,g_slab,gfield=', & 
!          minval(g_slab),minval(gfield),maxval(g_slab),maxval(gfield)
!        print*,'min,max,h_slab,hfield=', & 
!          minval(h_slab),minval(hfield),maxval(h_slab),maxval(hfield)
        else
         read(lun1)
         read(lun2)
       end if

       if(mod(islab,npe) == 0 .or. islab == num_slabs) then
         call mp_flush(1)
         call mpi_alltoallv(hfield,iscnt,isdisp,mpi_real8, &
              hall_loc(ista,jsta,islab_prev,n),ircnt,irdisp, & 
                               mpi_real8,mpi_comm_world,ierror)

         call mpi_alltoallv(gfield,iscnt,isdisp,mpi_real8, &
              gall_loc(ista,jsta,islab_prev,n),ircnt,irdisp, & 
                               mpi_real8,mpi_comm_world,ierror)

       islab_prev=islab+1
       end if
200 continue
    read(lun1) h_rad
    read(lun1) h_prec

    read(lun2) g_rad
    read(lun2) g_prec

    if (nrf1_rad > 0) then 
       hrad(:,n)=h_rad(:)
       grad(:,n)=g_rad(:)
    endif

    if (nrf1_prec > 0) then 
       hprec(:,n)=h_prec(:)
       gprec(:,n)=g_prec(:)
    endif

    close (lun1)
    close (lun2)
400 continue

!***********************************************************
!***********************************************************
!==>compute renormalization coefficients cnorm

   isub=(iend-ista+1)
   jsub=(jend-jsta+1)
   ijsub=isub*jsub

   allocate(slab1(1:isub,1:jsub))
   allocate(slab2(1:isub,1:jsub))

   allocate(hgdotaux(0:nv,0:nv))
   allocate(hgdot(0:nv,0:nv))
   
   hgdotaux(0:nv,0:nv)=0._r_double
   do n=0,nv
      do m=0,nv
         do islab=1,num_slabs
            do j=1,jsub
            do i=1,isub
               slab1(i,j)=hall_loc(i+ista-1,j+jsta-1,islab,n)
               slab2(i,j)=gall_loc(i+ista-1,j+jsta-1,islab,m)
            enddo
            enddo
            hgdotaux(n,m)=hgdotaux(n,m)+innerprod(slab1,slab2,ijsub)
!           print*,'n,m,islab,min,max,slab1,slab2=', & 
!           n,m,islab,minval(slab1),minval(slab2),maxval(slab1),maxval(slab2)
         enddo
      enddo
   enddo


   hgdot(0:nv,0:nv)=0._r_double
   call mpi_allreduce(hgdotaux,hgdot,(nv+1)*(nv+1),mpi_real8, & 
                      mpi_sum,mpi_comm_world,ierror)

  if (nrf1_rad > 0) then
     do n=0,nv
        do m=0,nv
           hgdot(n,m)=hgdot(n,m)+innerprod(hrad(1,n),grad(1,m),jpch*npred)
        enddo
     enddo
  endif

  if (nrf1_prec > 0) then
     do n=0,nv
        do m=0,nv
           hgdot(n,m)=hgdot(n,m)+innerprod(hprec(1,n),gprec(1,m),jtype*npredp)
           if (mype.eq.0) print*,'n,m,hgdot=',n,m,hgdot(n,m)
        enddo
     enddo
  endif

  allocate(cnorm(0:nv))
  do n=0,nv
     if (hgdot(n,n) .le. tiny_r_kind*1._r_double) then
       print*,'n,hgdot(n,n)=',n,hgdot(n,n)
       print*,'trouble with normalization--must abort'
       call abort
      else
       cnorm(n)=1._r_double/sqrt(hgdot(n,n))
      endif
  enddo
  if (mype .eq. 0) print*,'min,max,cnorm=',minval(cnorm),maxval(cnorm)

  deallocate(ijpts_in_sub)
  deallocate(ista_info)
  deallocate(iend_info)
  deallocate(jsta_info)
  deallocate(jend_info)
  deallocate(iaux0)
  deallocate(iaux1)
  deallocate(iaux2)
  deallocate(iaux3)
  deallocate(iaux4)
  deallocate(iscnt)
  deallocate(isdisp)
  deallocate(ircnt)
  deallocate(irdisp)
  deallocate(h_slab)
  deallocate(g_slab)
  deallocate(hfield)
  deallocate(gfield)
  deallocate(h_rad)
  deallocate(g_rad)
  deallocate(h_prec)
  deallocate(g_prec)
  deallocate(hgdot)
  deallocate(hgdotaux)

!******************************************************************************
!      step-II : build the T matrix                                           *
!******************************************************************************
  allocate(beta(0:nv))
  allocate(delta(0:nv))
  allocate(theta(0:nv))

  beta(0:nv)=zero
  delta(0:nv)=zero
  theta(0:nv)=zero

  delta(0)=1./alpha(0)
  do n=0,nv 
   if (n.lt.nv) then
     delta(n+1)=one/alpha(n+1) + gamma(n)/alpha(n)
     beta(n+1)=-cnorm(n)/cnorm(n+1)/alpha(n)
   endif
     if (n.gt.0) theta(n)=-gamma(n-1)/alpha(n-1)*cnorm(n)/cnorm(n-1)
  enddo

  if (mype.eq.0) then
    print*,'delta(0)=',delta(0)
    do n=1,nv
      print*,'n,delta(n),beta(n),theta(n)=',n,delta(n),beta(n),theta(n)
    enddo
   endif

  allocate(tmatrix(1:nv,1:nv))

  tmatrix(:,:)=zero

  do n=1,nv
     tmatrix(n,n)=delta(n-1)
     if (n.lt.nv) tmatrix(n,n+1)=theta(n)
     if (n.gt.1)  tmatrix(n,n-1)=beta(n-1)
  enddo

  do n=1,nv
  do m=n+1,nv
     tmatrix(n,m)=tmatrix(m,n)
  enddo
  enddo

!******************************************************************************
!      step-III : compute eigenvectors and eigenvalues of matrix
!******************************************************************************

  nsym=(nv*nv-nv)/2+nv
  allocate(atest(nsym))
    
  j=0
  do n=1,nv
     do i=1,n
        j=j+1
        atest(j)=tmatrix(i,n) 
     enddo 
  enddo
  if (mype.eq.0) print*,'nsym,j=',nsym,j

  allocate(rtest(1:nv,1:nv))

  rtest(:,:)=zero
  call eigen_3(atest,rtest,nv,0)

  allocate(eval(nv))
  
  j=0
  do n=1,nv
     do i=1,n
        j=j+1
        if (i.eq.n) eval(n)=atest(j)
     enddo 
     if (mype==0) print*,'n,eval=',n,eval(n)
  enddo

! check orthogonality of eigenvectors
  do n=1,nv
  do m=1,nv
     dotp=zero
     do k=1,nv
        dotp=dotp+rtest(k,m)*rtest(k,n)
     enddo
!     if (mype==0) print*,'orthogonality of T eigvectors, m,n,dotp=',m,n,dotp
  enddo
  enddo

!******************************************************************************
!      step-V : compute the eigenvectors of the Hessian. Store in gall_loc
!               grad, and gprec
!******************************************************************************
   allocate(hall_loc2(ista:iend,jsta:jend,num_slabs,1:nv))
   allocate(gall_loc2(ista:iend,jsta:jend,num_slabs,1:nv))

   hall_loc2(ista:iend,jsta:jend,1:num_slabs,1:nv)=zero
   gall_loc2(ista:iend,jsta:jend,1:num_slabs,1:nv)=zero

   do n=1,nv
      do islab=1,num_slabs
      do j=jsta,jend
      do i=ista,iend
         do k=1,nv
            hall_loc2(i,j,islab,n)=hall_loc2(i,j,islab,n)+ & 
                   hall_loc(i,j,islab,k-1)*cnorm(k-1)*rtest(k,n)
            gall_loc2(i,j,islab,n)=gall_loc2(i,j,islab,n)+ & 
                   gall_loc(i,j,islab,k-1)*cnorm(k-1)*rtest(k,n)
         enddo
      enddo
      enddo
      enddo
   enddo

   if (nrf1_rad > 0) then 
      allocate(hrad2(jpch*npred,1:nv))
      allocate(grad2(jpch*npred,1:nv))

      hrad2(:,:)=zero
      grad2(:,:)=zero

      do n=1,nv
         do j=1,jpch*npred
            do k=1,nv
               hrad2(j,n)=hrad2(j,n)+hrad(j,k-1)*cnorm(k-1)*rtest(k,n)
               grad2(j,n)=grad2(j,n)+grad(j,k-1)*cnorm(k-1)*rtest(k,n)
            enddo
         enddo
      enddo
   endif

   if (nrf1_prec > 0) then 
      allocate(hprec2(jtype*npredp,1:nv))
      allocate(gprec2(jtype*npredp,1:nv))

      hprec2(:,:)=zero
      gprec2(:,:)=zero

      do n=1,nv
         do j=1,jtype*npredp
            do k=1,nv
               hprec2(j,n)=hprec2(j,n)+hprec(j,k-1)*cnorm(k-1)*rtest(k,n)
               gprec2(j,n)=gprec2(j,n)+gprec(j,k-1)*cnorm(k-1)*rtest(k,n)
            enddo
         enddo
      enddo
   endif
!******************************************************************************
!      step-VI : construct vector of error variance. Store in 
!                hall_loc(:,:,:,0), hrad(:,0), and  hprec(:,0)
!******************************************************************************
   do n=1,nv
     dotp0=0._r_double
     do islab=1,num_slabs
        do j=1,jsub
        do i=1,isub
           slab1(i,j)=hall_loc2(i+ista-1,j+jsta-1,islab,n)
           slab2(i,j)=gall_loc2(i+ista-1,j+jsta-1,islab,n)
        enddo
        enddo
        dotp0=dotp0+innerprod(slab1,slab2,ijsub)
     enddo

     dotp=0._r_double
     call mpi_allreduce(dotp0,dotp,1,mpi_real8, & 
                      mpi_sum,mpi_comm_world,ierror)

     if (nrf1_rad > 0) then 
        dotp=dotp+innerprod(hrad2(1,n),grad2(1,n),jpch*npred)
     endif

     if (nrf1_prec > 0) then 
        dotp=dotp+innerprod(hprec2(1,n),gprec2(1,n),jtype*npredp)
     endif

     if (dotp .gt. tiny_r_kind*1._r_double) then
       hall_loc2(ista:iend,jsta:jend,1:num_slabs,n)= & 
                hall_loc2(ista:iend,jsta:jend,1:num_slabs,n)/sqrt(dotp)

       gall_loc2(ista:iend,jsta:jend,1:num_slabs,n)= & 
                gall_loc2(ista:iend,jsta:jend,1:num_slabs,n)/sqrt(dotp)

       if (nrf1_rad > 0) then 
           hrad2(:,n)=hrad2(:,n)/sqrt(dotp)
           grad2(:,n)=grad2(:,n)/sqrt(dotp)
       endif

       if (nrf1_prec > 0) then 
           hprec2(:,n)=hprec2(:,n)/sqrt(dotp)
           gprec2(:,n)=gprec2(:,n)/sqrt(dotp)
        endif

      else
       print*,'for Hessian, getting null eigevector #',n
     endif
   enddo

!test L2-orthogonality of eigenvectors of the hessian
   do n=1,nv
      do m=1,nv
        dotp0=0._r_double
        do islab=1,num_slabs
           do j=1,jsub
           do i=1,isub
              slab1(i,j)=hall_loc2(i+ista-1,j+jsta-1,islab,n)
              slab2(i,j)=gall_loc2(i+ista-1,j+jsta-1,islab,m)
           enddo
           enddo
           dotp0=dotp0+innerprod(slab1,slab2,ijsub)
        enddo

        dotp=0._r_double
        call mpi_allreduce(dotp0,dotp,1,mpi_real8, & 
                      mpi_sum,mpi_comm_world,ierror)
!        if (mype.eq.0)  print*,'n,m,dotprod=',n,m,dotp
      enddo
   enddo

!==>  Perform variable transformation to satisfy specific needs (RTMA, NMM, GFS). For 
!  example, RTMA needs u and v variances. Hence, convert eigenvectors of hessian from
!  (psi,chi) to (u,v). 

   allocate(psi(ista:iend,jsta:jend,nsig))

   allocate(chi(ista:iend,jsta:jend,nsig))

   allocate(t(ista:iend,jsta:jend,nsig))

   allocate(q(ista:iend,jsta:jend,nsig))

   allocate(oz(ista:iend,jsta:jend,nsig))

   allocate(cwmr(ista:iend,jsta:jend,nsig))

   allocate(sfcp(ista:iend,jsta:jend))

   allocate(sfct(ista:iend,jsta:jend))

   allocate(u(ista:iend,jsta:jend,nsig,nv))

   allocate(v(ista:iend,jsta:jend,nsig,nv))

   allocate(u2(ista:iend,jsta:jend,nsig,nv))

   allocate(v2(ista:iend,jsta:jend,nsig,nv))

   allocate(wdir(ista:iend,jsta:jend,nsig,nv))

   allocate(wdir2(ista:iend,jsta:jend,nsig,nv))

   if (nrf2_gust > 0) allocate(gust(ista:iend,jsta:jend))
   if (nrf2_vis >  0) allocate(vis(ista:iend,jsta:jend))
   if (nrf2_pblh > 0) allocate(pblh(ista:iend,jsta:jend))
   if (nrf2_dist > 0) allocate(dist(ista:iend,jsta:jend))

   ista2=max(1,ista-1)
   iend2=min(nx,iend+1)
   jsta2=max(1,jsta-1)
   jend2=min(ny,jend+1)

   allocate(tempa(nlon,nlat))

   allocate ( coeffx (jsta2:jend2,ista2:iend2) )   !transposed
   allocate ( coeffy (jsta2:jend2,ista2:iend2) )   !transposed

   allocate ( glon   (ista2:iend2,jsta2:jend2) )   !non-transposed

   open (33,file='rtma_latlon_mpfactor.dat',form='unformatted', &
             convert='Little_Endian')
   read(33)
   read(33) tempa  !glat
   read(33) tempa  !glon

   if (mype==0) print*,'in lanczos_anlerr: glon,min,max=',minval(tempa),maxval(tempa)

   do j=jsta2,jend2
   do i=ista2,iend2
      glon(i,j)=tempa(i,j)
      if (glon(i,j) < 0._r_single)  &
          glon(i,j)=glon(i,j)+360._r_single! convert to eastern longitude
   enddo
   enddo

   rewind(33)
   read(33) tempa    !mapfact

   do j=jsta2,jend2
   do i=ista2,iend2
      coeffx(j,i)=half/(ds0*tempa(i,j)*one)
      coeffy(j,i)=coeffx(j,i)
   enddo
   enddo
   close(33)

   amin8=minval(one/coeffx) ; amax8=maxval(one/coeffx)
   call mpi_allreduce(amin8,amin82,1,mpi_real8,mpi_min,mpi_comm_world,ierror)
   call mpi_allreduce(amax8,amax82,1,mpi_real8,mpi_max,mpi_comm_world,ierror)

   if (mype==0) then
       print*,'in lanczos_anlerr: 1/coeffx,min,max=',amin82,amax82
       print*,'in lanczos_anlerr: 1/coeffy,min,max=',amin82,amax82
   endif

   dg2rad=atan(one)/45._r_kind

   if (trim(cgrid)=='alaska')    xn=one
   if (trim(cgrid)=='akhres')    xn=one
   if (trim(cgrid)=='juneau')    xn=one
   if (trim(cgrid)=='conus')     xn=sin(r25*dg2rad)
   if (trim(cgrid)=='cohres')    xn=sin(r25*dg2rad)
   if (trim(cgrid)=='hrrr')      xn=sin(r38_5*dg2rad)
   if (trim(cgrid)=='dtc')       xn=sin(47.49*dg2rad)
   if (trim(cgrid)=='cohresext') xn=sin(r25*dg2rad)

   if (trim(cgrid)=='alaska')    elonv=r210
   if (trim(cgrid)=='akhres')    elonv=r210
   if (trim(cgrid)=='juneau')    elonv=r225
   if (trim(cgrid)=='conus')     elonv=r265
   if (trim(cgrid)=='cohres')    elonv=r265
   if (trim(cgrid)=='hrrr')      elonv=r262_5
   if (trim(cgrid)=='dtc')       elonv=256.0
   if (trim(cgrid)=='cohresext') elonv=r265

   if (leigenvector_write) then 
      open (14,file='hessian_eigenvectors.dat',form='unformatted')
      write(14) nv
      write(14) eval

      allocate(tempb8(nlon,nlat))
      allocate(tempc8(nlon,nlat))
   endif

   do n=1,nv
      kall=0
      do ifield=1,6
         if (ifield==5 .and. nrf3_oz <= 0) cycle 
         if (ifield==6 .and. nrf3_cw <= 0) cycle
         do k=1,nsig
            kall=kall+1
            do j=jsta,jend
               do i=ista,iend
                  if (ifield.eq.1) psi(i,j,k)=hall_loc2(i,j,kall,n)
                  if (ifield.eq.2) chi(i,j,k)=hall_loc2(i,j,kall,n)
                  if (ifield.eq.3) t(i,j,k)=hall_loc2(i,j,kall,n)
                  if (ifield.eq.4) q(i,j,k)=hall_loc2(i,j,kall,n)
                  if (ifield.eq.5) oz(i,j,k)=hall_loc2(i,j,kall,n)
                  if (ifield.eq.6) cwmr(i,j,k)=hall_loc2(i,j,kall,n)
               enddo
             enddo
           enddo
      enddo

      kall=kall+1
      do j=jsta,jend
         do i=ista,iend
            sfcp(i,j)=hall_loc2(i,j,kall,n)
         enddo
      enddo

      if (nrf2_sst > 0 ) then
        kall=kall+1
        sfct(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n)
      endif

      if (nrf2_gust > 0 ) then 
         kall=kall+1
         gust(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n)
      endif

      if (nrf2_vis > 0 ) then 
         kall=kall+1
         vis(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n)
      endif

      if (nrf2_pblh > 0 ) then 
         kall=kall+1
         pblh(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n)
      endif

      if (nrf2_dist > 0 ) then 
         kall=kall+1
         dist(ista:iend,jsta:jend)=hall_loc2(ista:iend,jsta:jend,kall,n)
      endif


      ! note: rad and prec parts are hrad2(:,n) and hprec2(:,n)

      allocate(uaux1(ista:iend,jsta:jend))
      allocate(vaux1(ista:iend,jsta:jend))
      allocate(uaux2(ista:iend,jsta:jend))
      allocate(vaux2(ista:iend,jsta:jend))

      do k=1,nsig         ! (psi,chi) ==> (u,v) conversion

         call psichi2uv_eigenvec(psi(:,:,k),chi(:,:,k), & 
              coeffx,coeffy,glon,uaux1,vaux1,uaux2,vaux2, &
              dg2rad,xn,elonv,k,n,ista,iend,jsta,jend, &
              ista2,iend2,jsta2,jend2,nlon,nlat,cgrid,mype,npe)

         do j=jsta,jend
         do i=ista,iend
            u(i,j,k,n)=uaux1(i,j)
            v(i,j,k,n)=vaux1(i,j)
            wdir(i,j,k,n)=(u(i,j,k,n)/vbar-ubar/vbar**2*v(i,j,k,n))/(one+ubar**2/vbar**2)/dg2rad
            u2(i,j,k,n)=uaux2(i,j)
            v2(i,j,k,n)=vaux2(i,j)
            wdir2(i,j,k,n)=(u2(i,j,k,n)/vbar-ubar/vbar**2*v2(i,j,k,n))/(one+ubar**2/vbar**2)/dg2rad
         enddo
         enddo
      enddo

      deallocate(uaux1)
      deallocate(vaux1)
      deallocate(uaux2)
      deallocate(vaux2)

      if (leigenvector_write) then 
         do ifield=1,8+6+4
            if ( ifield == 5  .and. nrf3_oz   <= 0 ) cycle
            if ( ifield == 6  .and. nrf3_cw   <= 0 ) cycle 
            if ( ifield == 15 .and. nrf2_gust <= 0 ) cycle 
            if ( ifield == 16 .and. nrf2_vis  <= 0 ) cycle 
            if ( ifield == 17 .and. nrf2_pblh <= 0 ) cycle 
            if ( ifield == 18 .and. nrf2_dist <= 0 ) cycle 

            nk=nsig
            if (ifield==7 .or. ifield==8 .or. ifield > 14) nk=1

            do k=1,nk
               tempb8=zero
               if (ifield == 1)  tempb8(ista:iend,jsta:jend) = psi(ista:iend,jsta:jend,k)
               if (ifield == 2)  tempb8(ista:iend,jsta:jend) = chi(ista:iend,jsta:jend,k)
               if (ifield == 3)  tempb8(ista:iend,jsta:jend) = t(ista:iend,jsta:jend,k)
               if (ifield == 4)  tempb8(ista:iend,jsta:jend) = q(ista:iend,jsta:jend,k)

               if (ifield == 5)  tempb8(ista:iend,jsta:jend) = oz(ista:iend,jsta:jend,k)
               if (ifield == 6)  tempb8(ista:iend,jsta:jend) = cwmr(ista:iend,jsta:jend,k)

               if (ifield == 7)  tempb8(ista:iend,jsta:jend) = sfcp(ista:iend,jsta:jend)
               if (ifield == 8)  tempb8(ista:iend,jsta:jend) = sfct(ista:iend,jsta:jend)

               if (ifield == 9)  tempb8(ista:iend,jsta:jend) = u(ista:iend,jsta:jend,k,n)
               if (ifield == 10) tempb8(ista:iend,jsta:jend) = v(ista:iend,jsta:jend,k,n)
               if (ifield == 11) tempb8(ista:iend,jsta:jend) = wdir(ista:iend,jsta:jend,k,n)
               if (ifield == 12) tempb8(ista:iend,jsta:jend) = u2(ista:iend,jsta:jend,k,n)
               if (ifield == 13) tempb8(ista:iend,jsta:jend) = v2(ista:iend,jsta:jend,k,n)
               if (ifield == 14) tempb8(ista:iend,jsta:jend) = wdir2(ista:iend,jsta:jend,k,n)
    
               if (ifield == 15) tempb8(ista:iend,jsta:jend) = gust(ista:iend,jsta:jend)
               if (ifield == 16) tempb8(ista:iend,jsta:jend) = vis(ista:iend,jsta:jend)
               if (ifield == 17) tempb8(ista:iend,jsta:jend) = pblh(ista:iend,jsta:jend)
               if (ifield == 18) tempb8(ista:iend,jsta:jend) = dist(ista:iend,jsta:jend)

               tempc8=zero
               call mpi_reduce(tempb8,tempc8,nlon*nlat,mpi_real8,&
                    mpi_sum,0,mpi_comm_world,ierror)
                    if (mype==0) write(14) tempc8 
            enddo
         enddo
      endif
   enddo

   if (leigenvector_write) then 
      deallocate(tempb8)
      deallocate(tempc8)
      close(14)
   endif
 
!==>Perform final matrix multiplication to get variances
   hall_loc(ista:iend,jsta:jend,1:num_pad,0)=zero
   do islab=1,num_slabs
      do j=jsta,jend
      do i=ista,iend
         do k=1,nv
           if (eval(k).gt.tiny_r_kind) then
              hall_loc(i,j,islab,0)=hall_loc(i,j,islab,0)+ & 
              hall_loc2(i,j,islab,k)*(-one/eval(k)+one/eval(k)**2-one/eval(k)**3)*hall_loc2(i,j,islab,k)*float(k*k)
!             hall_loc2(i,j,islab,k)*(one/eval(k)-one)*hall_loc2(i,j,islab,k)   !exact
           endif
         enddo
      enddo
      enddo
   enddo

   if (nrf1_rad > 0) then
      hrad(:,0)=zero
      do j=1,jpch*npred
         do k=1,nv
            hrad(j,0)=hrad(j,0)+hrad2(j,k)*(-one/eval(k)+one/eval(k)**2-one/eval(k)**3)*hrad2(j,k)*float(k*k)
!           hrad(j,0)=hrad(j,0)+hrad2(j,k)*(one/eval(k)-one)*hrad2(j,k)            !exact
         enddo
      enddo
   endif

   if (nrf1_prec > 0) then
      hprec(:,0)=zero
      do j=1,jtype*npredp
         do k=1,nv
            hprec(j,0)=hprec(j,0)+hprec2(j,k)*(-one/eval(k)+one/eval(k)**2-one/eval(k)**3)*hprec2(j,k)*float(k*k)
!           hprec(j,0)=hprec(j,0)+hprec2(j,k)*(one/eval(k)-one)*hprec2(j,k)        !exact
         enddo
      enddo
   endif

   allocate(usum(ista:iend,jsta:jend,1:nsig)) 
   allocate(vsum(ista:iend,jsta:jend,1:nsig)) 
   allocate(wdirsum(ista:iend,jsta:jend,1:nsig)) 
   allocate(usum2(ista:iend,jsta:jend,1:nsig)) 
   allocate(vsum2(ista:iend,jsta:jend,1:nsig)) 
   allocate(wdirsum2(ista:iend,jsta:jend,1:nsig)) 

   usum=zero
   vsum=zero
   wdirsum=zero
   usum2=zero
   vsum2=zero
   wdirsum2=zero

   do kk=1,nsig
      do j=jsta,jend
      do i=ista,iend
         do k=1,nv
           if (eval(k).gt.tiny_r_kind) then
              usum(i,j,kk)=usum(i,j,kk)+ & 
              u(i,j,kk,k)*(-one/eval(k)+one/eval(k)**2-one/eval(k)**3)*u(i,j,kk,k)*float(k*k)
!             u(i,j,kk,k)*(one/eval(k)-one)*u(i,j,kk,k)  !exact 

              vsum(i,j,kk)=vsum(i,j,kk)+ & 
              v(i,j,kk,k)*(-one/eval(k)+one/eval(k)**2-one/eval(k)**3)*v(i,j,kk,k)*float(k*k) 
!             v(i,j,kk,k)*(one/eval(k)-one)*v(i,j,kk,k)  !exact

              wdirsum(i,j,kk)=wdirsum(i,j,kk)+ & 
              wdir(i,j,kk,k)*(-one/eval(k)+one/eval(k)**2-one/eval(k)**3)*wdir(i,j,kk,k)*float(k*k) 
!             wdir(i,j,kk,k)*(one/eval(k)-one)*wdir(i,j,kk,k)  !exact

              usum2(i,j,kk)=usum2(i,j,kk)+ & 
              u2(i,j,kk,k)*(-one/eval(k)+one/eval(k)**2-one/eval(k)**3)*u2(i,j,kk,k)*float(k*k) 
!             u2(i,j,kk,k)*(one/eval(k)-one)*u2(i,j,kk,k)  !exact 

              vsum2(i,j,kk)=vsum2(i,j,kk)+ & 
              v2(i,j,kk,k)*(-one/eval(k)+one/eval(k)**2-one/eval(k)**3)*v2(i,j,kk,k)*float(k*k) 
!             v2(i,j,kk,k)*(one/eval(k)-one)*v2(i,j,kk,k)  !exact 

              wdirsum2(i,j,kk)=wdirsum2(i,j,kk)+ & 
              wdir2(i,j,kk,k)*(-one/eval(k)+one/eval(k)**2-one/eval(k)**3)*wdir2(i,j,kk,k)*float(k*k) 
!             wdir2(i,j,kk,k)*(one/eval(k)-one)*wdir2(i,j,kk,k)  !exact
           endif
         enddo
      enddo
      enddo
   enddo
!******************************************************************************
!      step-VIII : write out solution (valid only for twodvar_regional)
!******************************************************************************
   call create_errs_common(ista,iend,jsta,jend, & 
        nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist)  !allocates 
                                                 !psierr,chierr,uerr,verr,uerr2,verr2,
                                                 !wspderr,wdirerr,wdirerr2,terr,tderr, 
                                                 !qerr,perr,gusterr,viserr,pblherr,disterr

   do j=jsta,jend
   do i=ista,iend
      uerr(i,j)=usum(i,j,1)
      verr(i,j)=vsum(i,j,1)
      wdirerr(i,j)=wdirsum(i,j,1)

      uerr2(i,j)=usum2(i,j,1)
      verr2(i,j)=vsum2(i,j,1)
      wdirerr2(i,j)=wdirsum2(i,j,1)
   enddo
   enddo

   uerr=sqrt(max(tiny_single,-uerr))
   verr=sqrt(max(tiny_single,-verr))
   wdirerr=sqrt(max(tiny_single,-wdirerr))

   uerr2=sqrt(max(tiny_single,-uerr2))
   verr2=sqrt(max(tiny_single,-verr2))
   wdirerr2=sqrt(max(tiny_single,-wdirerr2))

   npointer(:)=-1
   n=0
   if (nrf3_sf   > 0 ) then ; n=n+1 ; npointer(1)  = n ; endif 
   if (nrf3_vp   > 0 ) then ; n=n+1 ; npointer(2)  = n ; endif
   if (nrf3_t    > 0 ) then ; n=n+1 ; npointer(3)  = n ; endif
   if (nrf3_q    > 0 ) then ; n=n+1 ; npointer(4)  = n ; endif
   if (nrf3_oz   > 0 ) then ; n=n+1 ; npointer(5)  = n ; endif
   if (nrf3_cw   > 0 ) then ; n=n+1 ; npointer(6)  = n ; endif
   if (nrf2_ps   > 0 ) then ; n=n+1 ; npointer(7)  = n ; endif
   if (nrf2_sst  > 0 ) then ; n=n+1 ; npointer(8)  = n ; endif
   if (nrf2_gust > 0 ) then ; n=n+1 ; npointer(9)  = n ; endif
   if (nrf2_vis  > 0 ) then ; n=n+1 ; npointer(10) = n ; endif
   if (nrf2_pblh > 0 ) then ; n=n+1 ; npointer(11) = n ; endif
   if (nrf2_dist > 0 ) then ; n=n+1 ; npointer(12) = n ; endif


   if (mype==0) open (15,file='errvar.dat',form='unformatted') 

   allocate(tempb(nlon,nlat))

   do n=1,num_slabs

      tempa=0._r_single
      tempb=0._r_single

      do j=jsta,jend
      do i=ista,iend
         tempb(i,j)=hall_loc(i,j,n,0)
      enddo
      enddo

      call mpi_allreduce(tempb,tempa,nlon*nlat,mpi_real4,&
           mpi_sum,mpi_comm_world,ierror)

      if (mype==0) write(15) tempa

      tempa=max(tiny_single,-tempa)
      tempa=sqrt(tempa)

      if (n==1) psierr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend)
      if (n==2) chierr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend)
      if (n==3) terr   (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend)
      if (n==4) qerr   (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend)

      if (n==npointer(7)) perr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend)

      if (nrf2_gust > 0 .and. n==npointer(9) ) gusterr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend)
      if (nrf2_vis  > 0 .and. n==npointer(10)) viserr  (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend)
      if (nrf2_pblh > 0 .and. n==npointer(11)) pblherr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend)
      if (nrf2_dist > 0 .and. n==npointer(12)) disterr (ista:iend,jsta:jend) = tempa(ista:iend,jsta:jend)

   enddo
   if (mype==0) close(15)

   call error_conversion(cgrid,cvbasedrecalibration,cvbasedcmodel0, &
              nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist,ista,iend,jsta,jend,nlon,nlat,mype,npe)

!  call error_conversion_ops(cgrid,cvbasedrecalibration,cvbasedcmodel0, &
!             nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist,ista,iend,jsta,jend,nlon,nlat,mype,npe)

   filename='errfield.dat'
   call writeout_errfields(filename,ista,iend,jsta,jend,nx,ny, &
                           nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist, & 
                           mype,npe)

   deallocate(tempa)
   deallocate(tempb)
   deallocate(glon)

   call destroy_errs_common(nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist)

                            !deallocates 
                            !psierr,chierr,uerr,verr,uerr2,verr2,
                            !wspderr,wdirerr,wdirerr2,terr,tderr, 
                            !qerr,perr,gusterr,viserr,pblherr,disterr

   deallocate(hall_loc)
   deallocate(gall_loc)
   deallocate(hall_loc2)
   deallocate(gall_loc2)

   if (nrf1_rad  > 0) then 
      deallocate(hrad)
      deallocate(hrad2)
      deallocate(grad)
      deallocate(grad2)
   endif

   if (nrf1_prec > 0) then 
       deallocate(hprec)
       deallocate(hprec2)
       deallocate(gprec)
       deallocate(gprec2)
   endif


   deallocate(alpha)
   deallocate(gamma)
   deallocate(beta)
   deallocate(delta)
   deallocate(theta)
   deallocate(cnorm)
   deallocate(tmatrix)
   deallocate(atest)
   deallocate(rtest)
   deallocate(eval)

   deallocate(slab1)
   deallocate(slab2)

   deallocate(psi)
   deallocate(chi)
   deallocate(t)
   deallocate(q)
   deallocate(oz)
   deallocate(cwmr)
   deallocate(sfcp)
   deallocate(sfct)
   deallocate(u)
   deallocate(v)
   deallocate(u2)
   deallocate(v2)
   deallocate(coeffx)
   deallocate(coeffy)
   deallocate(usum)
   deallocate(vsum)
   deallocate(usum2)
   deallocate(vsum2)
   deallocate(wdir)
   deallocate(wdir2)
   deallocate(wdirsum)
   deallocate(wdirsum2)
   if (nrf2_gust > 0) deallocate(gust)
   if (nrf2_vis  > 0) deallocate(vis)
   if (nrf2_pblh > 0) deallocate(pblh)
   if (nrf2_dist > 0) deallocate(dist)

!******************************************************************************
!==> force use of climatological error when explicitly so specified or when
!    number of gradient vectors is less than a chosen nvmin value
!******************************************************************************
3500 continue
   if (nv .lt. nvmin .or. climatology) then !use climatological error

       allocate(climerr(nlon,nlat))

       if (mype==0) then

          print*,' either number of iterations used in minimization is less than ', nvmin, &
                 'or climatology is set to .true.'

          print*,'skip lanczos procedure and use climatological error instead'

          open (16,file='errfield.dat_clim',form='unformatted')
          open (17,file='errfield.dat',form='unformatted')

          do n=1,11
             read(16) climerr 
             write(17) climerr 
          enddo
          if (nrf2_gust > 0 )  then ; read(16) climerr ; write(17) climerr ; endif
          if (nrf2_vis  > 0 )  then ; read(16) climerr ; write(17) climerr ; endif
          if (nrf2_pblh > 0 )  then ; read(16) climerr ; write(17) climerr ; endif
          if (nrf2_dist > 0 )  then ; read(16) climerr ; write(17) climerr ; endif

          close(16)
          close(17)
       endif

       deallocate(climerr)
       call mpi_barrier(mpi_comm_world,ierror)
   endif

!******************************************************************************
!      step-IX : create grib2 files for ges, anl, and anlerr
!******************************************************************************
   call cnv_to_grib2_mpi(cgrid,ista,iend,jsta,jend,nlon,nlat, & 
        lbiascor,nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist,mype,npe)

   call mpi_barrier(mpi_comm_world,ierror)


!******************************************************************************
!      step-X : generate observation text files + create dynamic reject lists
!******************************************************************************
   mypecutoff=miter0

   if (mype <=  mypecutoff) then 
      icolor=1
      ikey=mype 
     else
      icolor=2
      ikey=(mype-(mypecutoff+1))
   endif

   call mpi_comm_split(mpi_comm_world,icolor,ikey,mpi_comm_new,ierror)

   call mpi_comm_size(mpi_comm_new,npenew,ierror)
   call mpi_comm_rank(mpi_comm_new,mypenew,ierror)

   print*,'mpi_comm_new,npenew,mypenew=',mpi_comm_new,npenew,mypenew

   call mpi_barrier(mpi_comm_world,ierror)


   if (icolor==1) then
      do n=1,miter0+1
         if (mypenew==mod(n-1,npenew)) call get_ob_lists(n,miter0,cgrid, & 
                                       nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist)
      enddo
   endif

   if (npe <= mypecutoff+1) then
      if (icolor==1) call create_rjlist(cgrid,mypenew,npenew,mpi_comm_new)
    else
      if (icolor==2) call create_rjlist(cgrid,mypenew,npenew,mpi_comm_new)
   endif

   call mpi_barrier(mpi_comm_world,ierror)


!******************************************************************************
!      step-XI : perform optional cressman analysis and/or cross-validation
!******************************************************************************
   if (anlqlty) call anl_quality(mype,npe,ista,iend,jsta,jend,rjbuffer_km,ds0, &
                     nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist,lobjanl)

  if (anlqlty .or. cvbasedrecalibration) call destroy_cressanl_common

   
end subroutine lanczos_anlerr
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
 subroutine psichi2uv_eigenvec(psi,chi,coeffx,coeffy,glon, & 
                              uaux1,vaux1,uaux2,vaux2, & 
                              dg2rad,xn,elonv,k,n,ista,iend,jsta,jend, & 
                              ista2,iend2,jsta2,jend2,nlon,nlat,cgrid,mype,npe)
                                 
!****************************************************************************
! abstract: convert (psi,chi) eigenvectors into (u,v) eigenvectors          *
!                                                                           *
! program history log:                                                      *
!   2005-10-08  pondeca                                                     *
!                                                                           *
!****************************************************************************
  use mpi
  use kinds, only: i_kind,r_single,r_kind
  use constants, only: zero,one

  implicit none

! Declare passed variables
  character(60),intent(in):: cgrid
  integer(i_kind),intent(in):: mype,npe
  integer(i_kind),intent(in):: nlon,nlat                  !global dimensions
  integer(i_kind),intent(in):: ista,iend,jsta,jend        !tile dimensions
  integer(i_kind),intent(in):: ista2,iend2,jsta2,jend2    !memory dimensions
  integer(i_kind),intent(in):: k                          !vertical level #
  integer(i_kind),intent(in):: n                          !eigenvector number

  real(r_kind),intent(in)::dg2rad,xn,elonv

  real(r_single),intent(in):: glon(ista2:iend2,jsta2:jend2)

  real(r_kind),intent(in):: psi(ista:iend,jsta:jend) 
  real(r_kind),intent(in):: chi(ista:iend,jsta:jend)

  real(r_kind),intent(in):: coeffx(jsta2:jend2,ista2:iend2)   !transposed
  real(r_kind),intent(in):: coeffy(jsta2:jend2,ista2:iend2)   !transposed

  real(r_kind),intent(out):: uaux1(ista:iend,jsta:jend) 
  real(r_kind),intent(out):: vaux1(ista:iend,jsta:jend)
  real(r_kind),intent(out):: uaux2(ista:iend,jsta:jend) 
  real(r_kind),intent(out):: vaux2(ista:iend,jsta:jend)

! Declare local variables
  integer(i_kind) i,j,ierror
  real(r_kind) angle2,sinx2,cosx2
  real(r_kind) amin8,amax8,amin82,amax82
  real(r_kind),allocatable,dimension(:,:):: dtemp1,dtemp2
  real(r_kind),allocatable,dimension(:,:):: daux1,daux2,daux3,daux4
! 
!*****************************************************************************

  if (n==1) then
     open (10,file='psi_slabs.dat',form='unformatted')
     open (20,file='chi_slabs.dat',form='unformatted')

     open (11,file='u_slabs.dat',form='unformatted')
     open (21,file='v_slabs.dat',form='unformatted')
   else
     open (10,file='psi_slabs.dat',form='unformatted',position='append')
     open (20,file='chi_slabs.dat',form='unformatted',position='append')

     open (11,file='u_slabs.dat',form='unformatted',position='append')
     open (21,file='v_slabs.dat',form='unformatted',position='append')
  endif

  allocate(dtemp1(nlon,nlat))
  allocate(dtemp2(nlon,nlat))

  allocate(daux1(jsta2:jend2,ista2:iend2))     !transposed
  allocate(daux2(jsta2:jend2,ista2:iend2))     !transposed
  allocate(daux3(jsta2:jend2,ista2:iend2))     !transposed
  allocate(daux4(jsta2:jend2,ista2:iend2))     !transposed

  dtemp1(:,:)=zero
  do j=jsta,jend
  do i=ista,iend
     dtemp1(i,j)=psi(i,j)
  enddo
  enddo
  call mpi_allreduce(dtemp1,dtemp2,nlon*nlat,mpi_real8, &
                     mpi_sum,mpi_comm_world,ierror)

  do j=jsta2,jend2
  do i=ista2,iend2
     daux1(j,i)=dtemp2(i,j)    !holds transposed psi in memory subdomains
  enddo
  enddo
            !optional
            !------------------------------------------------------------------------
            if (mype==0) print*,'n,k,psi,min,max=',n,k,minval(dtemp2),maxval(dtemp2)
            if (mype==0) write(10) sngl(dtemp2)
            !------------------------------------------------------------------------
  

  dtemp1(:,:)=zero
  do j=jsta,jend
  do i=ista,iend
       dtemp1(i,j)=chi(i,j)
  enddo
  enddo

  call mpi_allreduce(dtemp1,dtemp2,nlon*nlat,mpi_real8, &
                       mpi_sum,mpi_comm_world,ierror)

  do j=jsta2,jend2
  do i=ista2,iend2
     daux2(j,i)=dtemp2(i,j)    !holds transposed chi in memory subdomains
  enddo
  enddo
            !optional
            !------------------------------------------------------------------------
            if (mype==0) print*,'n,k,chi,min,max=',n,k,minval(dtemp2),maxval(dtemp2)
            if (mype==0) write(20) sngl(dtemp2)
            !------------------------------------------------------------------------


  call psichi2uv_reg_V2(daux1,daux2,daux3,daux4,coeffx,coeffy, &
                        1,     nlat,  1,     nlon,  1, 1, &
                        jsta2, jend2, ista2, iend2, 1, 1, &
                        jsta , jend , ista , iend , 1, 1)


  deallocate(daux1)
  deallocate(daux2)

  allocate(daux1(ista2:iend2,jsta2:jend2))     !non-transposed
  allocate(daux2(ista2:iend2,jsta2:jend2))     !non-transposed

  daux1=transpose(daux3)
  daux2=transpose(daux4)

            !optional
            !------------------------------------------------------------------------
            dtemp1(:,:)=zero
            do j=jsta,jend
            do i=ista,iend
               dtemp1(i,j)=daux1(i,j)
            enddo
            enddo
            call mpi_allreduce(dtemp1,dtemp2,nlon*nlat,mpi_real8, &
                      mpi_sum,mpi_comm_world,ierror)

            if (mype==0) print*,'n,k,u,min,max=',n,k,minval(dtemp2),maxval(dtemp2)
            if (mype==0) write(11) sngl(dtemp2)


            dtemp1(:,:)=zero
            do j=jsta,jend
            do i=ista,iend
               dtemp1(i,j)=daux2(i,j)
            enddo
            enddo
            call mpi_allreduce(dtemp1,dtemp2,nlon*nlat,mpi_real8, &
                     mpi_sum,mpi_comm_world,ierror)
            if (mype==0) print*,'n,k,v,min,max=',n,k,minval(dtemp2),maxval(dtemp2)
            if (mype==0) write(21) sngl(dtemp2)
            !------------------------------------------------------------------------

  !get earth relative u,v components of eigenvectors and store in daux3,daux4

  deallocate(daux3)
  deallocate(daux4)

  allocate(daux3(ista2:iend2,jsta2:jend2))     !non-transposed
  allocate(daux4(ista2:iend2,jsta2:jend2))     !non-transposed

  daux3=daux1         
  daux4=daux2         

  if (trim(cgrid)=='conus' .or. trim(cgrid)=='alaska' .or. & 
      trim(cgrid)=='cohres' .or. trim(cgrid)=='akhres' .or. trim(cgrid)=='hrrr' .or. & 
      trim(cgrid)=='juneau' .or. trim(cgrid)=='cohresext') then
     do j=jsta2,jend2
     do i=ista2,iend2
        angle2=xn*(glon(i,j)*one-elonv)*dg2rad
        sinx2=sin(angle2)
        cosx2 = cos(angle2)
        daux3(i,j)=cosx2*daux1(i,j)-sinx2*daux2(i,j)
        daux4(i,j)=sinx2*daux1(i,j)+cosx2*daux2(i,j)
     enddo 
     enddo 
  endif

            !optional
            !------------------------------------------------------------------------
            amin8=minval(daux3) ; amax8=maxval(daux3)
            call mpi_allreduce(amin8,amin82,1,mpi_real8,mpi_min,mpi_comm_world,ierror)
            call mpi_allreduce(amax8,amax82,1,mpi_real8,mpi_max,mpi_comm_world,ierror)
            if (mype==0) print*,'n,k,uearth,min,max=',n,k,amin82,amax82

            amin8=minval(daux4) ; amax8=maxval(daux4)
            call mpi_allreduce(amin8,amin82,1,mpi_real8,mpi_min,mpi_comm_world,ierror)
            call mpi_allreduce(amax8,amax82,1,mpi_real8,mpi_max,mpi_comm_world,ierror)
            if (mype==0) print*,'n,k,vearth,min,max=',n,k,amin82,amax82
            !------------------------------------------------------------------------

  do j=jsta,jend
  do i=ista,iend
     uaux1(i,j)=daux1(i,j)
     vaux1(i,j)=daux2(i,j)
     uaux2(i,j)=daux3(i,j)
     vaux2(i,j)=daux4(i,j)
  enddo
  enddo

  close(10)
  close(20)
  close(11)
  close(21)

  deallocate(daux1)
  deallocate(daux2)
  deallocate(daux3)
  deallocate(daux4)

  deallocate(dtemp1)
  deallocate(dtemp2)
end subroutine psichi2uv_eigenvec 
!******************************************************************************
!******************************************************************************

       subroutine horiz_domain_partition(nx,ny,mype,npe, & 
                                   ista,iend,jsta,jend)
       use kinds, only:i_kind
       implicit none

       integer(i_kind) nx,ny,mype,npe,ista,iend,jsta,jend

       jsta=1
       jend=ny

       call para_range2(1,nx,npe,mype,ista,iend)

       return
       end subroutine horiz_domain_partition
!******************************************************************************
       subroutine para_range2(n1,n2,nprocs,irank,ista,iend)

       use kinds, only:i_kind
       implicit none
                                                                                                     
       integer(i_kind) n1,n2,nprocs,irank,ista,iend
       integer(i_kind) iwork1,iwork2
                                                                                                     
       iwork1 = (n2-n1+1) / nprocs
       iwork2 = mod(n2 - n1 + 1 , nprocs)
       ista = irank * iwork1 + n1 + min(irank, iwork2)
       iend = ista + iwork1 -1
       if (iwork2 > irank) iend = iend+1
       return
       end subroutine para_range2
!******************************************************************************
       subroutine para_range3(n1,n2,nparts,ista,iend)

       use kinds, only:i_kind
       implicit none

       integer(i_kind),intent(in)  :: n1,n2,nparts
       integer(i_kind),intent(out) :: ista(nparts),iend(nparts)
       integer(i_kind) iwork1,iwork2,n

       iwork1 = (n2-n1+1) / nparts
       iwork2 = mod(n2 - n1 + 1 , nparts)

       do n=0,nparts-1
          ista(n+1) = n * iwork1 + n1 + min(n, iwork2)
          iend(n+1) = ista(n+1) + iwork1 -1
          if (iwork2 > n) iend(n+1) = iend(n+1)+1
       enddo
       return
       end subroutine para_range3
!******************************************************************************
       subroutine vectorform(g,v,nx,ny,istart,iend,jstart,jend,nprocs)

       use kinds, only:r_kind,i_kind
       implicit none

!      Declare passed variables
       integer(i_kind) nx,ny,nprocs
       integer(i_kind) istart(nprocs)
       integer(i_kind) iend(nprocs)
       integer(i_kind) jstart(nprocs)
       integer(i_kind) jend(nprocs)
       real(r_kind) g(nx,ny), v(nx*ny) 

!      Declare local variables
       integer(i_kind) i,j,ij,n

       ij=0
       do n=1,nprocs
          do j=jstart(n),jend(n)
             do i=istart(n),iend(n)
                ij=ij+1
                v(ij)=g(i,j)
             enddo
           enddo
       enddo
!      print*,'in vectorform, ij=',ij

       return
       end subroutine vectorform
!******************************************************************************
real(r_double) function innerprod(dx,dy,n)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    innerprod   calculates dot product for data on subdomain
!   prgmmr: derber           org: np23                date: 2004-05-13
!
! abstract: calculates dot product for data on subdomain.  Note loops over
!           streamfunction, velocity potential, temperature, etc. Also, only 
!           over interior of subdomain.
!
! program history log:
!   2004-05-13  derber, document
!   2004-07-28  treadon - add only on use declarations; add intent in/out
!
!   input argument list:
!     dx       - input vector 1
!     dy       - input vector 2
!
!   output argument list
!     dplev    - dot product
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  
  use kinds, only: r_kind,r_double,i_kind
  use constants, only: zero
  implicit none

! Declare passed variables
  integer(i_kind) n
  real(r_kind),dimension(n),intent(in)::dx,dy

! Declare local variables
  integer(i_kind) i
  
  innerprod=0._r_double
  do i=1,n
     innerprod=innerprod+dx(i)*dy(i)
  end do
  
  return
end function innerprod
!******************************************************************************
       subroutine ax_equal_lambda_bx(a,b,evec,eval,n)
       use kinds, only:r_kind,r_double,i_kind
       use constants, only:tiny_r_kind
       implicit none

       integer(i_kind) n,i,j
       real(r_kind) a(n,n),b(n,n),evec(n,n),eval(n)
       real(r_kind) beta(n),aux(3*n)  !CHECK DIMENSION OF AUX
       complex(r_kind) alpha(n)
       complex(r_kind),allocatable,dimension(:,:)::z

       allocate(z(n,n))

       call dgegv(1,a,n,b,n,alpha,beta,z,n,n,aux,3*n)

       do j=1,n

          if (abs(beta(j)) .gt. tiny_r_kind) then
             eval(j)=dreal(alpha(j))/beta(j)  !CHECK DREAL
           else
             print*,'trouble in dgegv. eval=infinity for j=',j
             eval(j)=-9999._r_kind
           endif

           do i=1,n
              evec(i,j)=dreal(z(i,j)) !CHECK DREAL
           enddo
       enddo 

       deallocate(z)

       return
       end subroutine ax_equal_lambda_bx
!******************************************************************************
!******************************************************************************
       subroutine error_conversion(cgrid,cvbasedrecalibration,cvbasedcmodel0, &
                        igust,ivis,ipblh,idist,ista,iend,jsta,jend,nx,ny,mype,npe)

       use mpi
       use constants, only: tiny_single

       use errs_common, only: psierr,chierr,uerr,verr,uerr2,verr2, &
                              wspderr,wdirerr,wdirerr2,terr,tderr, &
                              qerr,perr,gusterr,viserr,pblherr,disterr

       implicit none

!declare passed variables

       character(60),intent(in):: cgrid
       integer(4),intent(in):: ista,iend,jsta,jend
       integer(4),intent(in):: nx,ny
       integer(4),intent(in):: mype,npe
       integer(4),intent(in):: igust,ivis,ipblh,idist
       logical,intent(in):: cvbasedrecalibration
       character(60),intent(in):: cvbasedcmodel0

!declare local parameters
       integer(4),parameter::nbckgfls=9
       integer(4),parameter::nerrfls=15
       real(4),parameter:: wdirberr0=15. !(dg)
       real(4),parameter::a=243.5
       real(4),parameter::alpha=440.8
       real(4),parameter::b=19.48
       real(4),parameter::eps=0.62197
       real(4),parameter::perrm=300.
       real(4),parameter::qerrm=6.e-3

!declare local variables
       integer(4) i,j,im,ip,jm,jp
       integer(4) i1,i2,j1,j2
       integer(4) n,ierror
       integer(4) rtime(6),nlon,nlat,nsig

       real(4) amin,amax,aave
       real(4) amin2,amax2,aave2
       real(4) amin3,amax3,aave3
       real(4) tempamin,tempamax,fieldmin
       real(4) wspderrmin,wspderrmax
       real(4) uberrmin,vberrmin
       real(4) tderrmin,tderrmax

       character(60) bckgfname(nbckgfls)
       real(4),allocatable,dimension(:,:)::psiberr
       real(4),allocatable,dimension(:,:)::chiberr
       real(4),allocatable,dimension(:,:)::tberr
       real(4),allocatable,dimension(:,:)::qberr
       real(4),allocatable,dimension(:,:)::pberr
       real(4),allocatable,dimension(:,:)::uberr
       real(4),allocatable,dimension(:,:)::vberr
       real(4),allocatable,dimension(:,:)::wdirberr
       real(4),allocatable,dimension(:,:)::gustberr
       real(4),allocatable,dimension(:,:)::visberr
       real(4),allocatable,dimension(:,:)::pblhberr
       real(4),allocatable,dimension(:,:)::distberr

       real(4),allocatable,dimension(:,:)::field
       real(4),allocatable,dimension(:,:)::tempa

       real(4) qs0,qanl0,panl0
       real(4) qv,dqv,e,de,dlne,terrmax

       real(4) vis0

       logical fexist

       real(4) psicverr,chicverr,tcverr,pcverr,qcverr,ucverr,vcverr, & 
               wdircverr,u2cverr,v2cverr,wdir2cverr, & 
               wspdcverr,tdcverr, & 
               gustcverr,viscverr,pblhcverr,distcverr

       character(7) cfldname(nerrfls) !error field names
       real(4) hoberr(nerrfls)        !roughly half the observation error
       real(4) avgcverr(nerrfls)      !sqrt of the cross-validation rmse
       real(4) errupper(nerrfls)      !the maximum allowable error
       real(4) hoberrwspd,errupperwspd
       real(4) hoberrtd,erruppertd
       real(4) errmax

       character(60) filename

       namelist/cverrorupdate/psicverr,chicverr,tcverr,pcverr, & 
                qcverr,ucverr,vcverr,wdircverr,u2cverr,v2cverr,wdir2cverr, & 
                wspdcverr,tdcverr, & 
                gustcverr,viscverr,pblhcverr,distcverr

       data cfldname(1)  /'psi'/    ;  data psicverr   /-999./
       data cfldname(2)  /'chi'/    ;  data chicverr   /-999./
       data cfldname(3)  /'t'/      ;  data tcverr     /2.48/
       data cfldname(4)  /'p'/      ;  data pcverr     /2.12/   !ps in hPa
       data cfldname(5)  /'q'/      ;  data qcverr     /0.0011/ !this is for true sphum (in kg/kg)
       data cfldname(6)  /'u'/      ;  data ucverr     /2.26/
       data cfldname(7)  /'v'/      ;  data vcverr     /3.25/
       data cfldname(8)  /'wdir'/   ;  data wdircverr  /20./
       data cfldname(9)  /'u2'/     ;  data u2cverr    /2.26/
       data cfldname(10) /'v2'/     ;  data v2cverr    /3.25/
       data cfldname(11) /'wdir2'/  ;  data wdir2cverr /20./
       data cfldname(12) /'gust'/   ;  data gustcverr  /3.25/
       data cfldname(13) /'vis'/    ;  data viscverr   /3000./
       data cfldname(14) /'pblh'/   ;  data pblhcverr  /-999./
       data cfldname(15) /'dist'/   ;  data distcverr  /-999./
                                       data wspdcverr  /2.60/
                                       data tdcverr    /5.0/

       data hoberr(1)  /-999./  ; data errupper(1)  /-999./
       data hoberr(2)  /-999./  ; data errupper(2)  /-999./
       data hoberr(3)  /0.5/    ; data errupper(3)  /4.0/
       data hoberr(4)  /0.5/    ; data errupper(4)  /3./      !ps in hPa
       data hoberr(5)  /0.0005/ ; data errupper(5)  /0.004/   !this is for true sphum in (kg/kg)
       data hoberr(6)  /0.5/    ; data errupper(6)  /3.0/
       data hoberr(7)  /0.5/    ; data errupper(7)  /3.0/
       data hoberr(8)  /5./     ; data errupper(8)  /30./
       data hoberr(9)  /0.5/    ; data errupper(9)  /3.0/
       data hoberr(10) /0.5/    ; data errupper(10) /3.0/
       data hoberr(11) /0.5/    ; data errupper(11) /30./
       data hoberr(12) /1.0/    ; data errupper(12) /5.0/
       data hoberr(13) /100./   ; data errupper(13) /4500./
       data hoberr(14) /-999./  ; data errupper(14) /-999./
       data hoberr(15) /-999./  ; data errupper(15) /-999./
       data hoberrwspd /0.5/    ; data errupperwspd /3.0/
       data hoberrtd   /1.0/    ; data erruppertd   /6.0/


       data bckgfname(1)/'bckgvar.dat_ps'/
       data bckgfname(2)/'bckgvar.dat_t'/
       data bckgfname(3)/'bckgvar.dat_q'/
       data bckgfname(4)/'bckgvar.dat_u'/
       data bckgfname(5)/'bckgvar.dat_v'/
       data bckgfname(6)/'bckgvar.dat_gust'/
       data bckgfname(7)/'bckgvar.dat_vis'/
       data bckgfname(8)/'bckgvar.dat_pblh'/
       data bckgfname(9)/'bckgvar.dat_dist'/
!==========================================================================================
!==> read in the background error standard deviations. they either come from derive_xbar
!==========================================================================================
       i1=ista
       i2=iend
       j1=jsta
       j2=jend

       allocate ( psiberr  (i1:i2,j1:j2) )
       allocate ( chiberr  (i1:i2,j1:j2) )
       allocate ( tberr    (i1:i2,j1:j2) )
       allocate ( qberr    (i1:i2,j1:j2) )
       allocate ( pberr    (i1:i2,j1:j2) )
       allocate ( uberr    (i1:i2,j1:j2) )
       allocate ( vberr    (i1:i2,j1:j2) )
       allocate ( wdirberr (i1:i2,j1:j2) )

       if (igust > 0 ) allocate ( gustberr (i1:i2,j1:j2) )
       if (ivis  > 0 ) allocate ( visberr  (i1:i2,j1:j2) )
       if (ipblh > 0 ) allocate ( pblhberr (i1:i2,j1:j2) )
       if (idist > 0 ) allocate ( distberr (i1:i2,j1:j2) )

       allocate(field(ny,nx))        !transposed

!==> psi
       open (94,file='bckgvar.dat_psi',form='unformatted')        
       read(94) field
       do j=j1,j2
       do i=i1,i2
          psiberr(i,j)=field(j,i)
       enddo
       enddo
       close(94)

       amin2=minval(field)
       amax2=maxval(field)
       aave2=sum(field)/float(nx*ny)

!==> chi
       open (94,file='bckgvar.dat_chi',form='unformatted')        
       read(94) field
       do j=j1,j2
       do i=i1,i2
          chiberr(i,j)=field(j,i)
       enddo
       enddo
       close(94)

       amin3=minval(field)
       amax3=maxval(field)
       aave3=sum(field)/float(nx*ny)

       deallocate(field)           !deallocate transposed field
       allocate(field(nx,ny))      !allocate non-transposed field

!============================================================================
!==>dump bckg errors for u,v,t,q,gust,vis,pblh,vis on disc. also dump
!   background satutation specific humidity qs
!============================================================================
!

       call derive_xbvar_mpi(cgrid,nx,ny,igust,ivis,ipblh,idist,mype,npe)

!
!============================================================================
!==>read in the background satutation specific humidity qs and compute the 
!   domain averaged value qs0. it is used to convert the bckg error for    
!   specific humidity into the background error for pseudo-relative humidty, 
!   which is the true analysis variable
!============================================================================
       open (94,file='bckg_qsat.dat',form='unformatted')
       read(94) field
       close(94)

       qs0=sum(field(:,:))/float(nx*ny)   !0.01
       if (mype==0) print*,'in error_conversion: qs0=',qs0

!============================================================================
!==>read in the bckg errors for u,v,t,q,gust,vis,pblh, and vis. also, 
!   assume constant value of backg error for wind direction
!============================================================================

       if (mype==0)  print*,'============= in error_conversion ================'
       if (mype==0) print*,'psiberr,min,max,avg='   , amin2,amax2,aave2
       if (mype==0) print*,'chiberr,min,max,avg='   , amin3,amax3,aave3

       wdirberr=wdirberr0 !assume const background error

       do n=1,nbckgfls
          if (n == 6 .and. igust <= 0 ) cycle
          if (n == 7 .and. ivis  <= 0 ) cycle
          if (n == 8 .and. ipblh <= 0 ) cycle
          if (n == 9 .and. idist <= 0 ) cycle

          if (mype==0) print*,'in error_conversion: n,bckgfname=', & 
                                                    n,trim(bckgfname(n))

          open (94,file=trim(bckgfname(n)),form='unformatted')
          read(94) field
          close(94)

          if (n==1) field=field/1000. !get in cbars
          if (n==3) field=field/qs0

          if (n==1) pberr    (i1:i2,j1:j2) = field (i1:i2,j1:j2)
          if (n==2) tberr    (i1:i2,j1:j2) = field (i1:i2,j1:j2)
          if (n==3) qberr    (i1:i2,j1:j2) = field (i1:i2,j1:j2)
          if (n==4) uberr    (i1:i2,j1:j2) = field (i1:i2,j1:j2) 
          if (n==5) vberr    (i1:i2,j1:j2) = field (i1:i2,j1:j2)
          if (n==6) gustberr (i1:i2,j1:j2) = field (i1:i2,j1:j2)
          if (n==7) visberr  (i1:i2,j1:j2) = field (i1:i2,j1:j2)
          if (n==8) pblhberr (i1:i2,j1:j2) = field (i1:i2,j1:j2)
          if (n==9) distberr (i1:i2,j1:j2) = field (i1:i2,j1:j2)

          amin=minval(field)
          amax=maxval(field)
          aave=sum(field)/float(nx*ny)

          if (mype==0)  then
             if (n==1) print*,'pberr,min,max,avg='   , amin,amax,aave
             if (n==2) print*,'tberr,min,max,avg='   , amin,amax,aave 
             if (n==3) print*,'qberr,min,max,avg='   , amin,amax,aave
             if (n==4) print*,'uberr,min,max,avg='   , amin,amax,aave
             if (n==5) then 
                 print*,'vberr,min,max,avg='   , amin,amax,aave
                 print*,'wdirberr,min,max,avg=',wdirberr0,wdirberr0,wdirberr0
             endif
             if (n==6) print*,'gustberr,min,max,avg=' , amin,amax,aave
             if (n==7) print*,'visberr,min,max,avg='  , amin,amax,aave
             if (n==8) print*,'pblhberr,min,max,avg=' , amin,amax,aave
             if (n==9) print*,'distberr,min,max,avg=' , amin,amax,aave
          endif

       enddo
       deallocate(field)

       if (mype==0) print*,'==================================================='

!==========================================================================================
!==> read in x the domain-averaged cross-validation standard deviation error info through
!    namelist if available. otherwise, use 'climatological' info from data statements above.
!    it's used to rescale the analysis uncertainty.
!==========================================================================================
       inquire(file='cverrorupdate_input',exist=fexist)
       if (fexist) then
           open (55,file='cverrorupdate_input',form='formatted')
           read(55,cverrorupdate)
           close(55)
       endif

       avgcverr(1)  = psicverr
       avgcverr(2)  = chicverr
       avgcverr(3)  = tcverr
       avgcverr(4)  = pcverr
       avgcverr(5)  = qcverr
       avgcverr(6)  = ucverr
       avgcverr(7)  = vcverr
       avgcverr(8)  = wdircverr
       avgcverr(9)  = u2cverr
       avgcverr(10) = v2cverr
       avgcverr(11) = wdir2cverr
       avgcverr(12) = gustcverr
       avgcverr(13) = viscverr
       avgcverr(14) = pblhcverr
       avgcverr(15) = distcverr

       !convert from mb to cbar
       hoberr(4)=hoberr(4)/10.
       errupper(4)=errupper(4)/10.
       avgcverr(4)=avgcverr(4)/10.

       !convert from sphm error to pseudo-rh error
       hoberr(5)=hoberr(5)/qs0
       errupper(5)=errupper(5)/qs0
       avgcverr(5)=avgcverr(5)/qs0

       !convert from vis error to log(10)vis
       vis0=1500.
       if (hoberr(13)   > 0.) hoberr(13)=hoberr(13)/(log(10.)*vis0)
       if (errupper(13) > 0.) errupper(13)=errupper(13)/(log(10.)*vis0)
       if (avgcverr(13) > 0.) avgcverr(13)=avgcverr(13)/(log(10.)*vis0)

       do n=1,nerrfls
          if (n == 12 .and. igust <= 0 ) cycle
          if (n == 13 .and. ivis  <= 0 ) cycle
          if (n == 14 .and. ipblh <= 0 ) cycle
          if (n == 15 .and. idist <= 0 ) cycle

          if (mype==0) & 
              print*,' in error_conversion: hoberr,avgcverr,errupper for ',cfldname(n),':', &
              hoberr(n),avgcverr(n),errupper(n)
       enddo

!==========================================================================================
!==> compute the analysis uncertainty
!==========================================================================================
       allocate(tempa(i1:i2,j1:j2))
       allocate(field(i1:i2,j1:j2))   !allocate as a subdomain field

       do n=1,nerrfls
          if (n == 12 .and. igust <= 0 ) cycle
          if (n == 13 .and. ivis  <= 0 ) cycle
          if (n == 14 .and. ipblh <= 0 ) cycle
          if (n == 15 .and. idist <= 0 ) cycle

          if (n==1)  then  ; field=psiberr   ; tempa=psierr   ; endif
          if (n==2)  then  ; field=chiberr   ; tempa=chierr   ; endif
          if (n==3)  then  ; field=tberr     ; tempa=terr     ; endif
          if (n==4)  then  ; field=pberr     ; tempa=perr     ; endif
          if (n==5)  then  ; field=qberr     ; tempa=qerr     ; endif
          if (n==6)  then  ; field=uberr     ; tempa=uerr     ; endif
          if (n==7)  then  ; field=vberr     ; tempa=verr     ; endif
          if (n==8)  then  ; field=wdirberr  ; tempa=wdirerr  ; endif
          if (n==9)  then  ; field=uberr     ; tempa=uerr2    ; endif !assume uberr2=uberr
          if (n==10) then  ; field=vberr     ; tempa=verr2    ; endif !assume vberr2=vberr
          if (n==11) then  ; field=wdirberr  ; tempa=wdirerr2 ; endif !assume wdirberr2=wdirberr
          if (n==12) then  ; field=gustberr  ; tempa=gusterr  ; endif
          if (n==13) then  ; field=visberr   ; tempa=viserr   ; endif
          if (n==14) then  ; field=pblhberr  ; tempa=pblherr  ; endif
          if (n==15) then  ; field=distberr  ; tempa=disterr  ; endif

          if (n==13) then                   !visibility exception. assume anlerr is
             tempa=sqrt(sqrt(sqrt(tempa)))  !this function of the 'std deviation' reduction
             amax=maxval(tempa)
             call mpi_allreduce(amax,tempamax,1,mpi_real,mpi_max,mpi_comm_world,ierror)
             field=tempamax
          endif

          amin=minval(tempa) ; amax=maxval(tempa)

          call mpi_allreduce(amin,tempamin,1,mpi_real,mpi_min,mpi_comm_world,ierror)
          call mpi_allreduce(amax,tempamax,1,mpi_real,mpi_max,mpi_comm_world,ierror)

          if (mype==0) print*,'in error_conversion: n,tempa,min,max=', & 
                               n,tempamin,tempamax

          tempa=field-tempa                                   !preferred way
!         tempa=field-sqrt(tempa)                             !alternate way
!         tempa=sqrt(max(0.000001,(field*field-tempa*tempa))) !most correct way

          amin=minval(tempa)
          call mpi_allreduce(amin,tempamin,1,mpi_real,mpi_min,mpi_comm_world,ierror)

          if (hoberr(n) > 0.) then 
             tempa=tempa-tempamin+hoberr(n) 
           else
             tempa(:,:)=max(0.,tempa(:,:))
          endif
          
          amin=minval(field)
          call mpi_allreduce(amin,fieldmin,1,mpi_real,mpi_min,mpi_comm_world,ierror)

          if (avgcverr(n) > 0.) then
             errmax=max(min(avgcverr(n),errupper(n)),fieldmin)
            else
             errmax=1.
          endif
          if (mype==0) print*,'in error_conversion: n, errmax=',n, errmax

          amax=maxval(tempa)
          call mpi_allreduce(amax,tempamax,1,mpi_real,mpi_max,mpi_comm_world,ierror)

          tempa=tempa/tempamax*errmax

          if (n==1)  psierr   = tempa
          if (n==2)  chierr   = tempa
          if (n==3)  terr     = tempa
          if (n==4)  perr     = tempa
          if (n==5)  qerr     = tempa
          if (n==6)  uerr     = tempa
          if (n==7)  verr     = tempa
          if (n==8)  wdirerr  = tempa
          if (n==9)  uerr2    = tempa
          if (n==10) verr2    = tempa
          if (n==11) wdirerr2 = tempa
          if (n==12) gusterr  = tempa 
          if (n==13) viserr   = tempa
          if (n==14) pblherr  = tempa
          if (n==15) disterr  = tempa
       enddo

       amin=minval(uberr)
       call mpi_allreduce(amin,uberrmin,1,mpi_real,mpi_min,mpi_comm_world,ierror)

       amin=minval(vberr)
       call mpi_allreduce(amin,vberrmin,1,mpi_real,mpi_min,mpi_comm_world,ierror)


       wspderr=sqrt(uerr*uerr+verr*verr)            !a very wild approach

       amin=minval(wspderr)
       call mpi_allreduce(amin,wspderrmin,1,mpi_real,mpi_min,mpi_comm_world,ierror)

       wspderr=wspderr-wspderrmin+hoberrwspd
       errmax=max(min(wspdcverr,errupperwspd),0.5*(uberrmin+vberrmin)) !kinda reasonable
       if (mype==0) print*,'in error_conversion: errmax for wspd=',errmax

       amax=maxval(wspderr)
       call mpi_allreduce(amax,wspderrmax,1,mpi_real,mpi_max,mpi_comm_world,ierror)

       wspderr=wspderr/wspderrmax*errmax

       perr=perr*1000.  !convert to Pa
       qerr=qerr*qs0    !from pseudo-rh to sphum

       if ( ivis > 0 ) viserr=viserr*(log(10.)*vis0) !from log(10) vis to vis

       deallocate(field)
!==========================================================================================
!==> derive dewpoint temperature error

       allocate(field(nx,ny))   !allocate as a global field

       open (53,file='siganl',form='unformatted')

       read(53)
       read(53) rtime,nlon,nlat,nsig
       read(53) field     !full record contains glat,dx
       read(53) field     !full record contains glon,dy


       read(53) field     !panl
       panl0=.5*(maxval(field)+minval(field))

       read(53) field     !terrain
       read(53) field     !temperture
       read(53) field     !qanl
       close(53)
       qanl0=.5*(maxval(abs(field))+minval(abs(field)))

       if (mype==0) print*,'in error_conversion: qanl0,panl0=',qanl0,panl0

       do j=j1,j2
       do i=i1,i2
         qv=qanl0/(1.-qanl0)
         dqv=qerr(i,j)/(1.-qanl0)**2
         e=panl0/100.*qv/(eps+qv)
         de=e/panl0*perr(i,j) + &
            eps*panl0/100./(eps+qv)**2*dqv
         dlne=de/e
         tderr(i,j)=(a*b-alpha)/(b-alog(e))**2*dlne
         tderr(i,j)=abs(tderr(i,j))
       enddo
       enddo

       amin=minval(tderr)
       call mpi_allreduce(amin,tderrmin,1,mpi_real,mpi_min,mpi_comm_world,ierror)

       tderr=tderr-tderrmin+hoberrtd
       errmax=max(min(tdcverr,erruppertd),0.) !replace 0. with mininum of bckg error for td once available
       if (mype==0) print*,'in error_conversion: errmax for td=',errmax

       amax=maxval(tderr)
       call mpi_allreduce(amax,tderrmax,1,mpi_real,mpi_max,mpi_comm_world,ierror)

       tderr=tderr/tderrmax*errmax

       !just in case ...
       perr=max(tiny_single,perr)
       terr=max(tiny_single,terr)
       tderr=max(tiny_single,tderr)
       uerr=max(tiny_single,uerr)
       verr=max(tiny_single,verr)
       qerr=max(tiny_single,qerr)

       if (cvbasedrecalibration) then
           filename='errfield.dat_precval'
           call writeout_errfields(filename,i1,i2,j1,j2,nx,ny, & 
                                   igust,ivis,ipblh,idist,mype,npe) !writes out global fields of 
                                                                    !perr,terr,tderr,uerr,verr 
                                                                    !qerr,wdirerr2,wspderr,uerr2,
                                                                    !verr2,wdirerr,gusterr,viserr,
                                                                    !pblherr,disterr

           call cvbasedrecal(cgrid,cvbasedcmodel0,i1,i2,j1,j2,nx,ny, &
                   igust,ivis,ipblh,idist,mype,npe)

       endif


       deallocate(psiberr)
       deallocate(chiberr)
       deallocate(tberr)
       deallocate(qberr)
       deallocate(pberr)
       deallocate(uberr)
       deallocate(vberr)
       deallocate(wdirberr)

       if (igust > 0 ) deallocate(gustberr)
       if (ivis  > 0 ) deallocate(visberr)
       if (ipblh > 0 ) deallocate(pblhberr)
       if (idist > 0 ) deallocate(distberr)

       deallocate(field)
       deallocate(tempa)

end subroutine error_conversion
!******************************************************************************
!******************************************************************************
       subroutine error_conversion_ops(cgrid,cvbasedrecalibration,cvbasedcmodel0, &
                        igust,ivis,ipblh,idist,ista,iend,jsta,jend,nx,ny,mype,npe)

       use mpi
       use constants, only: tiny_single

       use errs_common, only: psierr,chierr,uerr,verr,uerr2,verr2, &
                              wspderr,wdirerr,wdirerr2,terr,tderr, &
                              qerr,perr,gusterr,viserr,pblherr,disterr

       implicit none

!declare passed variables

       character(60),intent(in):: cgrid
       integer(4),intent(in):: ista,iend,jsta,jend
       integer(4),intent(in):: nx,ny
       integer(4),intent(in):: mype,npe
       integer(4),intent(in):: igust,ivis,ipblh,idist
       logical,intent(in):: cvbasedrecalibration
       character(60),intent(in):: cvbasedcmodel0

!declare local parameters
       integer(4),parameter::nbckgfls=9
       integer(4),parameter::nerrfls=15
       real(4),parameter:: wdirberr0=15. !(dg)
       real(4),parameter::a=243.5
       real(4),parameter::alpha=440.8
       real(4),parameter::b=19.48
       real(4),parameter::eps=0.62197

       real(4),parameter::fraction=0.65
       real(4),parameter::uerrm=4.0
       real(4),parameter::verrm=4.0
       real(4),parameter::wdirerrm=8.0
       real(4),parameter::terrm=4.0
       real(4),parameter::perrm=300.
       real(4),parameter::qerrm=6.e-3

!declare local variables
       integer(4) i,j,im,ip,jm,jp
       integer(4) i1,i2,j1,j2
       integer(4) n,ierror
       integer(4) rtime(6),nlon,nlat,nsig

       real(4) amin,amax,aave
       real(4) amin2,amax2,aave2
       real(4) amin3,amax3,aave3
       real(4) tempamin,tempamax,fieldmin
       real(4) wspderrmin,wspderrmax
       real(4) uberrmin,vberrmin
       real(4) tderrmin,tderrmax

       character(60) bckgfname(nbckgfls)
       real(4),allocatable,dimension(:,:)::psiberr
       real(4),allocatable,dimension(:,:)::chiberr
       real(4),allocatable,dimension(:,:)::tberr
       real(4),allocatable,dimension(:,:)::qberr
       real(4),allocatable,dimension(:,:)::pberr
       real(4),allocatable,dimension(:,:)::uberr
       real(4),allocatable,dimension(:,:)::vberr
       real(4),allocatable,dimension(:,:)::wdirberr
       real(4),allocatable,dimension(:,:)::gustberr
       real(4),allocatable,dimension(:,:)::visberr
       real(4),allocatable,dimension(:,:)::pblhberr
       real(4),allocatable,dimension(:,:)::distberr

       real(4),allocatable,dimension(:,:)::field
       real(4),allocatable,dimension(:,:)::tempa

       real(4) qs0,qanl0,panl0
       real(4) qv,dqv,e,de,dlne,terrmax

       real(4) vis0

       logical fexist

       real(4) psicverr,chicverr,tcverr,pcverr,qcverr,ucverr,vcverr, & 
               wdircverr,u2cverr,v2cverr,wdir2cverr, & 
               wspdcverr,tdcverr, & 
               gustcverr,viscverr,pblhcverr,distcverr

       character(7) cfldname(nerrfls) !error field names
       real(4) hoberr(nerrfls)        !roughly half the observation error
       real(4) avgcverr(nerrfls)      !sqrt of the cross-validation rmse
       real(4) errupper(nerrfls)      !the maximum allowable error
       real(4) hoberrwspd,errupperwspd
       real(4) hoberrtd,erruppertd
       real(4) errmax

       character(60) filename

       namelist/cverrorupdate/psicverr,chicverr,tcverr,pcverr, & 
                qcverr,ucverr,vcverr,wdircverr,u2cverr,v2cverr,wdir2cverr, & 
                wspdcverr,tdcverr, & 
                gustcverr,viscverr,pblhcverr,distcverr

       data cfldname(1)  /'psi'/    ;  data psicverr   /-999./
       data cfldname(2)  /'chi'/    ;  data chicverr   /-999./
       data cfldname(3)  /'t'/      ;  data tcverr     /2.48/
       data cfldname(4)  /'p'/      ;  data pcverr     /2.12/   !ps in hPa
       data cfldname(5)  /'q'/      ;  data qcverr     /0.0011/ !this is for true sphum (in kg/kg)
       data cfldname(6)  /'u'/      ;  data ucverr     /2.26/
       data cfldname(7)  /'v'/      ;  data vcverr     /3.25/
       data cfldname(8)  /'wdir'/   ;  data wdircverr  /20./
       data cfldname(9)  /'u2'/     ;  data u2cverr    /2.26/
       data cfldname(10) /'v2'/     ;  data v2cverr    /3.25/
       data cfldname(11) /'wdir2'/  ;  data wdir2cverr /20./
       data cfldname(12) /'gust'/   ;  data gustcverr  /3.25/
       data cfldname(13) /'vis'/    ;  data viscverr   /3000./
       data cfldname(14) /'pblh'/   ;  data pblhcverr  /-999./
       data cfldname(15) /'dist'/   ;  data distcverr  /-999./
                                       data wspdcverr  /2.60/
                                       data tdcverr    /5.0/

       data hoberr(1)  /-999./  ; data errupper(1)  /-999./
       data hoberr(2)  /-999./  ; data errupper(2)  /-999./
       data hoberr(3)  /0.5/    ; data errupper(3)  /4.0/
       data hoberr(4)  /0.5/    ; data errupper(4)  /3./      !ps in hPa
       data hoberr(5)  /0.0005/ ; data errupper(5)  /0.004/   !this is for true sphum in (kg/kg)
       data hoberr(6)  /0.5/    ; data errupper(6)  /3.0/
       data hoberr(7)  /0.5/    ; data errupper(7)  /3.0/
       data hoberr(8)  /5./     ; data errupper(8)  /30./
       data hoberr(9)  /0.5/    ; data errupper(9)  /3.0/
       data hoberr(10) /0.5/    ; data errupper(10) /3.0/
       data hoberr(11) /0.5/    ; data errupper(11) /30./
       data hoberr(12) /1.0/    ; data errupper(12) /5.0/
       data hoberr(13) /100./   ; data errupper(13) /4500./
       data hoberr(14) /-999./  ; data errupper(14) /-999./
       data hoberr(15) /-999./  ; data errupper(15) /-999./
       data hoberrwspd /0.5/    ; data errupperwspd /3.0/
       data hoberrtd   /1.0/    ; data erruppertd   /6.0/


       data bckgfname(1)/'bckgvar.dat_ps'/
       data bckgfname(2)/'bckgvar.dat_t'/
       data bckgfname(3)/'bckgvar.dat_q'/
       data bckgfname(4)/'bckgvar.dat_u'/
       data bckgfname(5)/'bckgvar.dat_v'/
       data bckgfname(6)/'bckgvar.dat_gust'/
       data bckgfname(7)/'bckgvar.dat_vis'/
       data bckgfname(8)/'bckgvar.dat_pblh'/
       data bckgfname(9)/'bckgvar.dat_dist'/
!==========================================================================================
!==> read in the background error standard deviations. they either come from derive_xbar
!==========================================================================================
       i1=ista
       i2=iend
       j1=jsta
       j2=jend

       allocate ( psiberr  (i1:i2,j1:j2) )
       allocate ( chiberr  (i1:i2,j1:j2) )
       allocate ( tberr    (i1:i2,j1:j2) )
       allocate ( qberr    (i1:i2,j1:j2) )
       allocate ( pberr    (i1:i2,j1:j2) )
       allocate ( uberr    (i1:i2,j1:j2) )
       allocate ( vberr    (i1:i2,j1:j2) )
       allocate ( wdirberr (i1:i2,j1:j2) )

       if (igust > 0 ) allocate ( gustberr (i1:i2,j1:j2) )
       if (ivis  > 0 ) allocate ( visberr  (i1:i2,j1:j2) )
       if (ipblh > 0 ) allocate ( pblhberr (i1:i2,j1:j2) )
       if (idist > 0 ) allocate ( distberr (i1:i2,j1:j2) )

       allocate(field(ny,nx))        !transposed

!==> psi
       open (94,file='bckgvar.dat_psi',form='unformatted')        
       read(94) field
       do j=j1,j2
       do i=i1,i2
          psiberr(i,j)=field(j,i)
       enddo
       enddo
       close(94)

       amin2=minval(field)
       amax2=maxval(field)
       aave2=sum(field)/float(nx*ny)

!==> chi
       open (94,file='bckgvar.dat_chi',form='unformatted')        
       read(94) field
       do j=j1,j2
       do i=i1,i2
          chiberr(i,j)=field(j,i)
       enddo
       enddo
       close(94)

       amin3=minval(field)
       amax3=maxval(field)
       aave3=sum(field)/float(nx*ny)

       deallocate(field)           !deallocate transposed field
       allocate(field(nx,ny))      !allocate non-transposed field

!============================================================================
!==>dump bckg errors for u,v,t,q,gust,vis,pblh,vis on disc. also dump
!   background satutation specific humidity qs
!============================================================================
!

       call derive_xbvar_mpi(cgrid,nx,ny,igust,ivis,ipblh,idist,mype,npe)

!
!============================================================================
!==>read in the background satutation specific humidity qs and compute the 
!   domain averaged value qs0. it is used to convert the bckg error for    
!   specific humidity into the background error for pseudo-relative humidty, 
!   which is the true analysis variable
!============================================================================
       open (94,file='bckg_qsat.dat',form='unformatted')
       read(94) field
       close(94)

       qs0=sum(field(:,:))/float(nx*ny)   !0.01
       if (mype==0) print*,'in error_conversion: qs0=',qs0

!============================================================================
!==>read in the bckg errors for u,v,t,q,gust,vis,pblh, and vis. also, 
!   assume constant value of backg error for wind direction
!============================================================================

       if (mype==0)  print*,'============= in error_conversion ================'
       if (mype==0) print*,'psiberr,min,max,avg='   , amin2,amax2,aave2
       if (mype==0) print*,'chiberr,min,max,avg='   , amin3,amax3,aave3

       wdirberr=wdirberr0 !assume const background error

       do n=1,nbckgfls
          if (n == 6 .and. igust <= 0 ) cycle
          if (n == 7 .and. ivis  <= 0 ) cycle
          if (n == 8 .and. ipblh <= 0 ) cycle
          if (n == 9 .and. idist <= 0 ) cycle

          if (mype==0) print*,'in error_conversion: n,bckgfname=', & 
                                                    n,trim(bckgfname(n))

          open (94,file=trim(bckgfname(n)),form='unformatted')
          read(94) field
          close(94)

          if (n==1) field=field/1000. !get in cbars
          if (n==3) field=field/qs0

          if (n==1) pberr    (i1:i2,j1:j2) = field (i1:i2,j1:j2)
          if (n==2) tberr    (i1:i2,j1:j2) = field (i1:i2,j1:j2)
          if (n==3) qberr    (i1:i2,j1:j2) = field (i1:i2,j1:j2)
          if (n==4) uberr    (i1:i2,j1:j2) = field (i1:i2,j1:j2) 
          if (n==5) vberr    (i1:i2,j1:j2) = field (i1:i2,j1:j2)
          if (n==6) gustberr (i1:i2,j1:j2) = field (i1:i2,j1:j2)
          if (n==7) visberr  (i1:i2,j1:j2) = field (i1:i2,j1:j2)
          if (n==8) pblhberr (i1:i2,j1:j2) = field (i1:i2,j1:j2)
          if (n==9) distberr (i1:i2,j1:j2) = field (i1:i2,j1:j2)

          amin=minval(field)
          amax=maxval(field)
          aave=sum(field)/float(nx*ny)

          if (mype==0)  then
             if (n==1) print*,'pberr,min,max,avg='   , amin,amax,aave
             if (n==2) print*,'tberr,min,max,avg='   , amin,amax,aave 
             if (n==3) print*,'qberr,min,max,avg='   , amin,amax,aave
             if (n==4) print*,'uberr,min,max,avg='   , amin,amax,aave
             if (n==5) then 
                 print*,'vberr,min,max,avg='   , amin,amax,aave
                 print*,'wdirberr,min,max,avg=',wdirberr0,wdirberr0,wdirberr0
             endif
             if (n==6) print*,'gustberr,min,max,avg=' , amin,amax,aave
             if (n==7) print*,'visberr,min,max,avg='  , amin,amax,aave
             if (n==8) print*,'pblhberr,min,max,avg=' , amin,amax,aave
             if (n==9) print*,'distberr,min,max,avg=' , amin,amax,aave
          endif

       enddo
       deallocate(field)

       if (mype==0) print*,'==================================================='

!==========================================================================================
!==> read in x the domain-averaged cross-validation standard deviation error info through
!    namelist if available. otherwise, use 'climatological' info from data statements above.
!    it's used to rescale the analysis uncertainty.
!==========================================================================================
       inquire(file='cverrorupdate_input',exist=fexist)
       if (fexist) then
           open (55,file='cverrorupdate_input',form='formatted')
           read(55,cverrorupdate)
           close(55)
       endif

       avgcverr(1)  = psicverr
       avgcverr(2)  = chicverr
       avgcverr(3)  = tcverr
       avgcverr(4)  = pcverr
       avgcverr(5)  = qcverr
       avgcverr(6)  = ucverr
       avgcverr(7)  = vcverr
       avgcverr(8)  = wdircverr
       avgcverr(9)  = u2cverr
       avgcverr(10) = v2cverr
       avgcverr(11) = wdir2cverr
       avgcverr(12) = gustcverr
       avgcverr(13) = viscverr
       avgcverr(14) = pblhcverr
       avgcverr(15) = distcverr

       !convert from mb to cbar
       hoberr(4)=hoberr(4)/10.
       errupper(4)=errupper(4)/10.
       avgcverr(4)=avgcverr(4)/10.

       !convert from sphm error to pseudo-rh error
       hoberr(5)=hoberr(5)/qs0
       errupper(5)=errupper(5)/qs0
       avgcverr(5)=avgcverr(5)/qs0

       !convert from vis error to log(10)vis
       vis0=1500.
       if (hoberr(13)   > 0.) hoberr(13)=hoberr(13)/(log(10.)*vis0)
       if (errupper(13) > 0.) errupper(13)=errupper(13)/(log(10.)*vis0)
       if (avgcverr(13) > 0.) avgcverr(13)=avgcverr(13)/(log(10.)*vis0)

       do n=1,nerrfls
          if (n == 12 .and. igust <= 0 ) cycle
          if (n == 13 .and. ivis  <= 0 ) cycle
          if (n == 14 .and. ipblh <= 0 ) cycle
          if (n == 15 .and. idist <= 0 ) cycle

          if (mype==0) & 
              print*,' in error_conversion: hoberr,avgcverr,errupper for ',cfldname(n),':', &
              hoberr(n),avgcverr(n),errupper(n)
       enddo

!==========================================================================================
!==> compute the analysis uncertainty
!==========================================================================================
       allocate(tempa(i1:i2,j1:j2))
       allocate(field(i1:i2,j1:j2))   !allocate as a subdomain field

       do n=1,nerrfls
          if (n == 12 .and. igust <= 0 ) cycle
          if (n == 13 .and. ivis  <= 0 ) cycle
          if (n == 14 .and. ipblh <= 0 ) cycle
          if (n == 15 .and. idist <= 0 ) cycle

          if (n==1)  then  ; field=psiberr   ; tempa=psierr   ; endif
          if (n==2)  then  ; field=chiberr   ; tempa=chierr   ; endif
          if (n==3)  then  ; field=tberr     ; tempa=terr     ; endif
          if (n==4)  then  ; field=pberr     ; tempa=perr     ; endif
          if (n==5)  then  ; field=qberr     ; tempa=qerr     ; endif
          if (n==6)  then  ; field=uberr     ; tempa=uerr     ; endif
          if (n==7)  then  ; field=vberr     ; tempa=verr     ; endif
          if (n==8)  then  ; field=wdirberr  ; tempa=wdirerr  ; endif
          if (n==9)  then  ; field=uberr     ; tempa=uerr2    ; endif !assume uberr2=uberr
          if (n==10) then  ; field=vberr     ; tempa=verr2    ; endif !assume vberr2=vberr
          if (n==11) then  ; field=wdirberr  ; tempa=wdirerr2 ; endif !assume wdirberr2=wdirberr
          if (n==12) then  ; field=gustberr  ; tempa=gusterr  ; endif
          if (n==13) then  ; field=visberr   ; tempa=viserr   ; endif
          if (n==14) then  ; field=pblhberr  ; tempa=pblherr  ; endif
          if (n==15) then  ; field=distberr  ; tempa=disterr  ; endif

          if (n==13) then                   !visibility exception. assume anlerr is
             tempa=sqrt(sqrt(sqrt(tempa)))  !this function of the 'std deviation' reduction
             amax=maxval(tempa)
             call mpi_allreduce(amax,tempamax,1,mpi_real,mpi_max,mpi_comm_world,ierror)
             field=tempamax
          endif

          amin=minval(tempa) ; amax=maxval(tempa)

          call mpi_allreduce(amin,tempamin,1,mpi_real,mpi_min,mpi_comm_world,ierror)
          call mpi_allreduce(amax,tempamax,1,mpi_real,mpi_max,mpi_comm_world,ierror)

          if (mype==0) print*,'in error_conversion: n,tempa,min,max=', & 
                               n,tempamin,tempamax


          tempa=field-tempa/tempamax*fraction*field           !current ops way
!         tempa=field-tempa                                   !preferred way
!         tempa=field-sqrt(tempa)                             !alternate way
!         tempa=sqrt(max(0.000001,(field*field-tempa*tempa))) !most correct way

!         amin=minval(tempa)
!         call mpi_allreduce(amin,tempamin,1,mpi_real,mpi_min,mpi_comm_world,ierror)

!         if (hoberr(n) > 0.) then 
!            tempa=tempa-tempamin+hoberr(n) 
!          else
!            tempa(:,:)=max(0.,tempa(:,:))
!         endif
!         
!         amin=minval(field)
!         call mpi_allreduce(amin,fieldmin,1,mpi_real,mpi_min,mpi_comm_world,ierror)

!         if (avgcverr(n) > 0.) then
!            errmax=max(min(avgcverr(n),errupper(n)),fieldmin)
!           else
!            errmax=1.
!         endif
!         if (mype==0) print*,'in error_conversion: n, errmax=',n, errmax

          amax=maxval(tempa)
          call mpi_allreduce(amax,tempamax,1,mpi_real,mpi_max,mpi_comm_world,ierror)

!         tempa=tempa/tempamax*errmax

          if (n==1)  psierr   = tempa
          if (n==2)  chierr   = tempa
          if (n==3)  terr     = tempa*terrm/tempamax
          if (n==4)  perr     = tempa*perrm/tempamax   !*1000. !convert to Pa / perr is already in Pa
          if (n==5)  qerr     = tempa*qerrm/tempamax   !qs0   !from pseudo-rh to sphum /qerr is already in sphm units
          if (n==6)  uerr     = tempa*uerrm/tempamax
          if (n==7)  verr     = tempa*verrm/tempamax
          if (n==8)  wdirerr  = tempa*wdirerrm/tempamax
          if (n==9)  uerr2    = tempa*uerrm/tempamax
          if (n==10) verr2    = tempa*verrm/tempamax
          if (n==11) wdirerr2 = tempa*wdirerrm/tempamax
          if (n==12) gusterr  = tempa 
          if (n==13) viserr   = tempa*(log(10.)*vis0) !from log(10) vis to vis
          if (n==14) pblherr  = tempa
          if (n==15) disterr  = tempa
       enddo

!      amin=minval(uberr)
!      call mpi_allreduce(amin,uberrmin,1,mpi_real,mpi_min,mpi_comm_world,ierror)

!      amin=minval(vberr)
!      call mpi_allreduce(amin,vberrmin,1,mpi_real,mpi_min,mpi_comm_world,ierror)


       wspderr=sqrt(uerr*uerr+verr*verr)            !a very wild approach

!      amin=minval(wspderr)
!      call mpi_allreduce(amin,wspderrmin,1,mpi_real,mpi_min,mpi_comm_world,ierror)

!      wspderr=wspderr-wspderrmin+hoberrwspd
!      errmax=max(min(wspdcverr,errupperwspd),0.5*(uberrmin+vberrmin)) !kinda reasonable
!      if (mype==0) print*,'in error_conversion: errmax for wspd=',errmax

!      amax=maxval(wspderr)
!      call mpi_allreduce(amax,wspderrmax,1,mpi_real,mpi_max,mpi_comm_world,ierror)

!      wspderr=wspderr/wspderrmax*errmax

!      perr=perr*1000.  !convert to Pa
!      qerr=qerr*qs0    !from pseudo-rh to sphum
!      viserr=viserr*(log(10.)*vis0) !from log(10) vis to vis

       deallocate(field)
!==========================================================================================
!==> derive dewpoint temperature error

       allocate(field(nx,ny))   !allocate as a global field

       open (53,file='siganl',form='unformatted')

       read(53) rtime,nlon,nlat,nsig
       read(53) field     !full record contains glat,dx
       read(53) field     !full record contains glon,dy


       read(53) field     !panl
       panl0=.5*(maxval(field)+minval(field))

       read(53) field     !terrain
       read(53) field     !temperture
       read(53) field     !qanl
       close(53)
       qanl0=.5*(maxval(abs(field))+minval(abs(field)))

       if (mype==0) print*,'in error_conversion: qanl0,panl0=',qanl0,panl0

       do j=j1,j2
       do i=i1,i2
         qv=qanl0/(1.-qanl0)
         dqv=qerr(i,j)/(1.-qanl0)**2
         e=panl0/100.*qv/(eps+qv)
         de=e/panl0*perr(i,j) + &
            eps*panl0/100./(eps+qv)**2*dqv
         dlne=de/e
         tderr(i,j)=(a*b-alpha)/(b-alog(e))**2*dlne
         tderr(i,j)=abs(tderr(i,j))
       enddo
       enddo

!      amin=minval(tderr)
!      call mpi_allreduce(amin,tderrmin,1,mpi_real,mpi_min,mpi_comm_world,ierror)

!      tderr=tderr-tderrmin+hoberrtd
!      errmax=max(min(tdcverr,erruppertd),0.) !replace 0. with mininum of bckg error for td once available
!      if (mype==0) print*,'in error_conversion: errmax for td=',errmax

       amax=maxval(tderr)
       call mpi_allreduce(amax,tderrmax,1,mpi_real,mpi_max,mpi_comm_world,ierror)

       tderr=tderr/tderrmax*terrm*1.2

       !just in case ...
       perr=max(tiny_single,perr)
       terr=max(tiny_single,terr)
       tderr=max(tiny_single,tderr)
       uerr=max(tiny_single,uerr)
       verr=max(tiny_single,verr)
       qerr=max(tiny_single,qerr)

       if (cvbasedrecalibration) then
           filename='errfield.dat_precval'
           call writeout_errfields(filename,i1,i2,j1,j2,nx,ny, & 
                                   igust,ivis,ipblh,idist,mype,npe) !writes out global fields of 
                                                                    !perr,terr,tderr,uerr,verr 
                                                                    !qerr,wdirerr2,wspderr,uerr2,
                                                                    !verr2,wdirerr,gusterr,viserr,
                                                                    !pblherr,disterr

           call cvbasedrecal(cgrid,cvbasedcmodel0,i1,i2,j1,j2,nx,ny, &
                   igust,ivis,ipblh,idist,mype,npe)

       endif


       deallocate(psiberr)
       deallocate(chiberr)
       deallocate(tberr)
       deallocate(qberr)
       deallocate(pberr)
       deallocate(uberr)
       deallocate(vberr)
       deallocate(wdirberr)

       if (igust > 0 ) deallocate(gustberr)
       if (ivis  > 0 ) deallocate(visberr)
       if (ipblh > 0 ) deallocate(pblhberr)
       if (idist > 0 ) deallocate(distberr)

       deallocate(field)
       deallocate(tempa)

end subroutine error_conversion_ops
!*********************************************************************
!*********************************************************************
       subroutine writeout_errfields(filename,i1,i2,j1,j2,nx,ny, & 
                                     igust,ivis,ipblh,idist,mype,npe)

       use mpi

       use errs_common, only: uerr,verr,uerr2,verr2, &
                              wspderr,wdirerr,wdirerr2,terr,tderr, &
                              qerr,perr,gusterr,viserr,pblherr,disterr

       implicit none

!declare passed variables

       character(60),intent(in):: filename
       integer(4),intent(in):: mype,npe
       integer(4),intent(in):: i1,i2,j1,j2
       integer(4),intent(in):: igust,ivis,ipblh,idist
       integer(4),intent(in):: nx,ny

!declare local variables
       integer(4) nn,ierror

       real(4),allocatable,dimension(:,:):: slab1,slab2
!---------------------------------------------------------------------
!
       if (mype==0)open (17,file=trim(filename),form='unformatted')

       allocate(slab1(nx,ny))
       allocate(slab2(nx,ny))

       do nn=1,15

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

          slab1(:,:)=0.
          if (nn==1)  slab1 (i1:i2,j1:j2) = perr     (i1:i2,j1:j2)
          if (nn==2)  slab1 (i1:i2,j1:j2) = terr     (i1:i2,j1:j2)
          if (nn==3)  slab1 (i1:i2,j1:j2) = tderr    (i1:i2,j1:j2)
          if (nn==4)  slab1 (i1:i2,j1:j2) = uerr     (i1:i2,j1:j2)
          if (nn==5)  slab1 (i1:i2,j1:j2) = verr     (i1:i2,j1:j2)
          if (nn==6)  slab1 (i1:i2,j1:j2) = qerr     (i1:i2,j1:j2)
          if (nn==7)  slab1 (i1:i2,j1:j2) = wdirerr2 (i1:i2,j1:j2)
          if (nn==8)  slab1 (i1:i2,j1:j2) = wspderr  (i1:i2,j1:j2)
          if (nn==9)  slab1 (i1:i2,j1:j2) = uerr2    (i1:i2,j1:j2)
          if (nn==10) slab1 (i1:i2,j1:j2) = verr2    (i1:i2,j1:j2)
          if (nn==11) slab1 (i1:i2,j1:j2) = wdirerr  (i1:i2,j1:j2)
          if (nn==12) slab1 (i1:i2,j1:j2) = gusterr  (i1:i2,j1:j2)
          if (nn==13) slab1 (i1:i2,j1:j2) = viserr   (i1:i2,j1:j2)
          if (nn==14) slab1 (i1:i2,j1:j2) = pblherr  (i1:i2,j1:j2)
          if (nn==15) slab1( i1:i2,j1:j2) = disterr  (i1:i2,j1:j2)

          call mpi_allreduce(slab1,slab2,nx*ny, &
                   mpi_real,mpi_sum,mpi_comm_world,ierror)

          if (mype==0) write(17) slab2
       enddo
       close(17)

       call mpi_barrier(mpi_comm_world,ierror)

       deallocate(slab1)
       deallocate(slab2)
end subroutine writeout_errfields
!*********************************************************************
!*********************************************************************
!=======================================================================
!-------------------------------------------------------------------------------------
! SUBROUTINE 'EIGEN' FROM THE IBM SCIENTIFIC SUBROUTINE PACKAGE.
!
! NOTE:  TO CONFORM WITH THE FORTRAN77 STANDARD, DUMMY ARRAY DIMENSIONS
!        (1) HAVE BEEN CHANGED TO (*).
!
!     ..................................................................
!
!        SUBROUTINE EIGEN_3
!
!        PURPOSE
!           COMPUTE EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC
!           MATRIX
!
!        USAGE
!           CALL EIGEN_3(A,R,N,MV)
!
!        DESCRIPTION OF PARAMETERS
!           A - ORIGINAL MATRIX (SYMMETRIC), DESTROYED IN COMPUTATION.
!               RESULTANT EIGENVALUES ARE DEVELOPED IN DIAGONAL OF
!               MATRIX A IN DESCENDING ORDER.
!           R - RESULTANT MATRIX OF EIGENVECTORS (STORED COLUMNWISE,
!               IN SAME SEQUENCE AS EIGENVALUES)
!           N - ORDER OF MATRICES A AND R
!           MV- INPUT CODE
!                   0   COMPUTE EIGENVALUES AND EIGENVECTORS
!                   1   COMPUTE EIGENVALUES ONLY (R NEED NOT BE
!                       DIMENSIONED BUT MUST STILL APPEAR IN CALLING
!                       SEQUENCE)
!
!        REMARKS
!           ORIGINAL MATRIX A MUST BE REAL SYMMETRIC (STORAGE MODE=1)
!           MATRIX A CANNOT BE IN THE SAME LOCATION AS MATRIX R
!
!        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
!           NONE
!
!        METHOD
!           DIAGONALIZATION METHOD ORIGINATED BY JACOBI AND ADAPTED
!           BY VON NEUMANN FOR LARGE COMPUTERS AS FOUND IN 'MATHEMATICAL
!           METHODS FOR DIGITAL COMPUTERS', EDITED BY A. RALSTON AND
!           H.S. WILF, JOHN WILEY AND SONS, NEW YORK, 1962, CHAPTER 7
!
!     ..................................................................
!
      SUBROUTINE EIGEN_3(A,R,N,MV)
      use kinds, only: r_kind,i_kind,r_single,r_double,i_long
      integer(i_kind) N,MV
      real(r_kind) A(*),R(*),ANORM,ANRMX,THR,X,Y,SINX,SINX2,COSX,&
                   COSX2,SINCS,RANGE
!
!        ...............................................................
!
!        IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
!        C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
!        STATEMENT WHICH FOLLOWS.
!
!     DOUBLE PRECISION A,R,ANORM,ANRMX,THR,X,Y,SINX,SINX2,COSX,
!    1                 COSX2,SINCS,RANGE
!
!        THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
!        APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
!        ROUTINE.
!
!        THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO
!        CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS.  SQRT IN STATEMENTS
!        40, 68, 75, AND 78 MUST BE CHANGED TO DSQRT.  ABS IN STATEMENT
!        62 MUST BE CHANGED TO DABS. THE CONSTANT IN STATEMENT 5 SHOULD
!        BE CHANGED TO 1.0D-12.
!
!        ...............................................................
!
!        GENERATE IDENTITY MATRIX
!
    5 RANGE=1.0E-12_r_kind
      IF(MV-1) 10,25,10
   10 IQ=-N
      DO 20 J=1,N
      IQ=IQ+N
      DO 20 I=1,N
      IJ=IQ+I
      R(IJ)=0.0
      IF(I-J) 20,15,20
   15 R(IJ)=1.0
   20 CONTINUE
!
!        COMPUTE INITIAL AND FINAL NORMS (ANORM AND ANORMX)
!
   25 ANORM=0.0
      DO 35 I=1,N
      DO 35 J=I,N
      IF(I-J) 30,35,30
   30 IA=I+(J*J-J)/2
      ANORM=ANORM+A(IA)*A(IA)
   35 CONTINUE
      IF(ANORM) 165,165,40
   40 ANORM=1.414*SQRT(ANORM)
      ANRMX=ANORM*RANGE/FLOAT(N)
!
!        INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR
!
      IND=0
      THR=ANORM
   45 THR=THR/FLOAT(N)
   50 L=1
   55 M=L+1
!
!        COMPUTE SIN AND COS
!
   60 MQ=(M*M-M)/2
      LQ=(L*L-L)/2
      LM=L+MQ
   62 IF( ABS(A(LM))-THR) 130,65,65
   65 IND=1
      LL=L+LQ
      MM=M+MQ
      X=0.5*(A(LL)-A(MM))
   68 Y=-A(LM)/ SQRT(A(LM)*A(LM)+X*X)
      IF(X) 70,75,75
   70 Y=-Y
   75 SINX=Y/ SQRT(2.0*(1.0+( SQRT(1.0-Y*Y))))
      SINX2=SINX*SINX
   78 COSX= SQRT(1.0-SINX2)
      COSX2=COSX*COSX
      SINCS =SINX*COSX
!
!        ROTATE L AND M COLUMNS
!
      ILQ=N*(L-1)
      IMQ=N*(M-1)
      DO 125 I=1,N
      IQ=(I*I-I)/2
      IF(I-L) 80,115,80
   80 IF(I-M) 85,115,90
   85 IM=I+MQ
      GO TO 95
   90 IM=M+IQ
   95 IF(I-L) 100,105,105
  100 IL=I+LQ
      GO TO 110
  105 IL=L+IQ
  110 X=A(IL)*COSX-A(IM)*SINX
      A(IM)=A(IL)*SINX+A(IM)*COSX
      A(IL)=X
  115 IF(MV-1) 120,125,120
  120 ILR=ILQ+I
      IMR=IMQ+I
      X=R(ILR)*COSX-R(IMR)*SINX
      R(IMR)=R(ILR)*SINX+R(IMR)*COSX
      R(ILR)=X
  125 CONTINUE
      X=2.0*A(LM)*SINCS
      Y=A(LL)*COSX2+A(MM)*SINX2-X
      X=A(LL)*SINX2+A(MM)*COSX2+X
      A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2)
      A(LL)=Y
      A(MM)=X
!
!        TESTS FOR COMPLETION
!
!        TEST FOR M = LAST COLUMN
!
  130 IF(M-N) 135,140,135
  135 M=M+1
      GO TO 60
!
!        TEST FOR L = SECOND FROM LAST COLUMN
!
  140 IF(L-(N-1)) 145,150,145
  145 L=L+1
      GO TO 55
  150 IF(IND-1) 160,155,160
  155 IND=0
      GO TO 50
!
!        COMPARE THRESHOLD WITH FINAL NORM
!
  160 IF(THR-ANRMX) 165,165,45
!
!        SORT EIGENVALUES AND EIGENVECTORS
!
  165 IQ=-N
      DO 185 I=1,N
      IQ=IQ+N
      LL=I+(I*I-I)/2
      JQ=N*(I-2)
      DO 185 J=I,N
      JQ=JQ+N
      MM=J+(J*J-J)/2
      IF(A(LL)-A(MM)) 170,185,185
  170 X=A(LL)
      A(LL)=A(MM)
      A(MM)=X
      IF(MV-1) 175,185,175
  175 DO 180 K=1,N
      ILR=IQ+K
      IMR=JQ+K
      X=R(ILR)
      R(ILR)=R(IMR)
  180 R(IMR)=X
  185 CONTINUE
      RETURN
      END
!
!=======================================================================
       subroutine para_range(n1,n2,nprocs,irank,ista,iend)
       implicit none

       integer(4) n1,n2,nprocs,irank,ista,iend
       integer(4) iwork1,iwork2
       
       iwork1 = (n2-n1+1) / nprocs
       iwork2 = mod(n2 - n1 + 1 , nprocs)
       ista = irank * iwork1 + n1 + min(irank, iwork2)
       iend = ista + iwork1 -1
       if (iwork2 > irank) iend = iend+1
       return
       end
!------------------------------------------------------
       subroutine ob_loc(xob,yob,nobs,nvar)

       implicit none

       integer(4) nobs,nvar
       real(4) xob(nobs)
       real(4) yob(nobs)

       integer(4) n,lun

       character*8 filename
       character*5 clun1,clun2
       character*7 clun11,clun22
 
       if (nvar.eq.1) filename='u'
       if (nvar.eq.2) filename='v'
       if (nvar.eq.3) filename='t'
       if (nvar.eq.4) filename='p'
       if (nvar.eq.5) filename='q'

       lun=88
       open (lun,file='plot_ob_loc_'//trim(filename)//'.gs', & 
                 form='formatted')
       write(lun,*) "'set display color white'"
       write(lun,*) "'clear'"
       write(lun,*) "'open anlerr.des'"
       write(lun,*) "'set mproj off'"
       if (trim(filename).eq.'u') write(lun,*) "'display uerr'"
       if (trim(filename).eq.'v') write(lun,*) "'display verr'"
       if (trim(filename).eq.'t') write(lun,*) "'display terr'"
       if (trim(filename).eq.'p') write(lun,*) "'display perr'"
       if (trim(filename).eq.'q') write(lun,*) "'display qerr'"
       write(lun,*) 
       write(lun,*) 

       do 100 n=1,nobs
         if (n.le.9) then
           write (clun1,"(i1.1)") n
           write (clun2,"(i1.1)") n
         endif
         if (n.ge.10.and.n.le.99) then
           write (clun1,"(i2.2)") n
           write (clun2,"(i2.2)") n
         endif
         if (n.ge.100.and.n.le.999) then
           write (clun1,"(i3.3)") n
           write (clun2,"(i3.3)") n
         endif
         if (n.ge.1000.and.n.le.9999) then
           write (clun1,"(i4.4)") n
           write (clun2,"(i4.4)") n
         endif
         if (n.ge.10000.and.n.le.99999) then
           write (clun1,"(i5.5)") n
           write (clun2,"(i5.5)") n
         endif
         if (n.ge.100000.and.n.le.999999) then
           write (clun1,"(i6.6)") n
           write (clun2,"(i6.6)") n
         endif
         if (n.ge.1000000.and.n.le.9999999) then
           write (clun1,"(i7.7)") n
           write (clun2,"(i7.7)") n
         endif
         clun11='a.'//trim(clun1)
         clun22='b.'//trim(clun2)
         write(lun,"(1x,a7,3x,'=',3x,f12.6)") clun11,xob(n)
         write(lun,"(1x,a7,3x,'=',3x,f12.6)") clun22,yob(n)
100    continue

       write(lun,*) 
       write(lun,*) 
       write(lun,*) 'll=1'
       write(lun,*) 'while ( ll <',nobs+1,')'
       write(lun,*) "    'q gr2xy 'a.ll' 'b.ll"
       write(lun,*) '    line=sublin(result,1)'
       write(lun,*) '    xval=subwrd(line,3)'
       write(lun,*) '    yval=subwrd(line,6)'
       write(lun,*) "    'draw mark 3 'xval' 'yval' 0.06'"
       write(lun,*) '    ll=ll+1'
       write(lun,*) ' endwhile'

       return
       end
!------------------------------------------------------
subroutine smther_one(g1,is,ie,js,je,ns)

!      apply 1-2-1 smoother in each direction of data slab

! use kinds,only: r_single,i_long
  implicit none

  integer(4)  is, ie, js, je
  integer(4)  i,j,l,ip,im,jp,jm
  integer(4), intent(in) :: ns

  real(4), dimension(is:ie, js:je), intent(inout) :: g1
                                   !  on input: original data slab
                                   !  on ouput: filtered data slab
 

  real(4), allocatable:: g2(:,:)

   allocate(g2(is:ie,js:je))
   do l=1,ns

     do j=js,je
      do i=is,ie
       ip=min(i+1,ie) ; im=max(is,i-1)
         g2(i,j)=.25*(g1(ip,j)+g1(im,j))+.5*g1(i,j)
      end do
     end do

     do i=is,ie
      do j=js,je
       jp=min(j+1,je) ; jm=max(js,j-1)
       g1(i,j)=.25*(g2(i,jp)+g2(i,jm))+.5*g2(i,j)
      end do
     end do

   end do
   deallocate(g2)

   return
end subroutine smther_one
!------------------------------------------------------
       subroutine fetch_anlincs(itime,pinc,tinc,qinc, & 
                                uinc,vinc,tdinc, & 
                                is,ie,js,je,nx,ny,mype)

       use mpi
       implicit none

!Declare passed variables
       integer(4),intent(in):: is,ie,js,je
       integer(4),intent(in):: nx,ny,mype
       integer(4),intent(out):: itime(6)
       real(4),dimension(is:ie,js:je),intent(out):: & 
                         pinc,tinc,qinc,uinc,vinc,tdinc

!Declare local variables
       integer(4) nlon,nlat,nsig
       integer(4) i,j
       real(4),allocatable,dimension(:,:)::field1,field2
       real(4),allocatable,dimension(:,:)::t1,p1,q1,td1
       real(4),allocatable,dimension(:,:)::t2,p2,q2,td2
!--------------------------------------------------------------------

       if (mype==0) print*,'in fetch_anlincs: nx,ny=',nx,ny

       allocate(field1(nx,ny)) 
       allocate(field2(nx,ny)) 

       allocate ( t1  (is:ie,js:je) ) 
       allocate ( p1  (is:ie,js:je) ) 
       allocate ( q1  (is:ie,js:je) ) 
       allocate ( td1 (is:ie,js:je) ) 
       allocate ( t2  (is:ie,js:je) ) 
       allocate ( p2  (is:ie,js:je) ) 
       allocate ( q2  (is:ie,js:je) ) 
       allocate ( td2 (is:ie,js:je) ) 

       open (52,file='sigges',form='unformatted')
       open (53,file='siganl',form='unformatted')

       read(52)
       read(52) itime,nlon,nlat,nsig
       write(*,*) itime,nlon,nlat,nsig
       read(52) field1,field2! dx,dy
       read(52) field1! glat
       read(52) field2! glon

       print*,'in fetch_anlincs / from first guess'
       print*,'in fetch_anlincs / itime=',itime
       print*,'in fetch_anlincs / nlon,nlat,nsig=',nlon,nlat,nsig
       print*,'**********************************************'

       read(53)
       read(53) itime,nlon,nlat,nsig
       read(53) field1,field2! dx,dy
       read(53) field1! glat
       read(53) field2! glon

       print*,'in fetch_anlincs / from analysis'
       print*,'in fetch_anlincs / itime=',itime
       print*,'in fetch_anlincs / nlon,nlat,nsig=',nlon,nlat,nsig
       print*,'**********************************************'

       read(52) field1  !psfc 
       read(53) field2  !psfc 

       do j=js,je
       do i=is,ie
          pinc(i,j)=field2(i,j)-field1(i,j)
          p1(i,j)=field1(i,j)
          p2(i,j)=field2(i,j)
       enddo
       enddo

       read(52) field1  !fis 
       read(53) field2  !fis 

       read(52) field1  !t 
       read(53) field2  !t 

       do j=js,je
       do i=is,ie
          tinc(i,j)=field2(i,j)-field1(i,j)
          t1(i,j)=field1(i,j)
          t2(i,j)=field2(i,j)
       enddo
       enddo

       read(52) field1  !q 
       read(53) field2  !q 

       do j=js,je
       do i=is,ie
          qinc(i,j)=field2(i,j)-field1(i,j)
          q1(i,j)=field1(i,j)
          q2(i,j)=field2(i,j)
       enddo
       enddo

       read(52) field1  !u 
       read(53) field2  !u 

       do j=js,je
       do i=is,ie
          uinc(i,j)=field2(i,j)-field1(i,j)
       enddo
       enddo

       read(52) field1  !v 
       read(53) field2  !v 

       do j=js,je
       do i=is,ie
          vinc(i,j)=field2(i,j)-field1(i,j)
       enddo
       enddo

       call get_dewpt(p1,q1,t1,td1,is,ie,js,je)
       call get_dewpt(p2,q2,t2,td2,is,ie,js,je)

       !dew point increment
       do j=js,je
       do i=is,ie
          tdinc(i,j) = max ( -20., td2(i,j)-td1(i,j) )
       enddo
       enddo

       close(52)
       close(53)

       deallocate(field1)
       deallocate(field2)

       deallocate(t1)
       deallocate(p1)
       deallocate(q1)
       deallocate(td1)
       deallocate(t2)
       deallocate(p2)
       deallocate(q2)
       deallocate(td2)

       end subroutine fetch_anlincs
!------------------------------------------------------
!------------------------------------------------------
      subroutine bilinear_2d0v2(rffcst,ix1,ix2,jx1,jx2,rfobs,xx,yy)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    bilinear_2d0v2
!   prgmmr:
!
! abstract:
!
!
!   input argument list:
!    rffcst               - model grid value
!    ix,jx
!    xx,yy                - define coordinates in grid units
!                         of point for which interpolation is
!                         performed
!
!   output argument list:
!    rfobs                - interpolated value
!
! notes:
!
!     i+1,j |          | i+1,j+1
!         --+----------+---
!           |          | dym
!           |    *     + -
!           |   x,y    | dy
!           |          |
!         --+----+-----+---
!        i,j|<dx>|<dxm>| i,j+1
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block
      implicit none

!declare passed variables
      integer(4),intent(in   ) :: ix1,ix2,jx1,jx2
      real(4) ,intent(in   ) :: rffcst(ix1:ix2,jx1:jx2)
      real(4) ,intent(in   ) :: xx,yy
      real(4) ,intent(  out) :: rfobs

!declare local variables
      integer(4) i,j,ip,jp
      real(4) dx,dy,dxm,dym

      i  = ifix(yy)
      j  = ifix(xx)
      
      dx = xx - float(j)
      dy = yy - float(i)
      dxm= 1.0-dx
      dym= 1.0-dy
 
      i=min(max(ix1,i),ix2) ; j=min(max(jx1,j),jx2)
      ip=min(ix2,i+1)     ; jp=min(jx2,j+1) 

      rfobs=dxm*(dym*rffcst(i,j)+dy*rffcst(ip,j)) &
               + dx *(dym*rffcst(i,jp)+dy*rffcst(ip,jp))

      return
end subroutine bilinear_2d0v2
!------------------------------------------------------
!------------------------------------------------------
 subroutine anl_quality(mype,npe,ista,iend,jsta,jend,rjbuffer_km,ds0, &
            nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist,lobjanl)
  
!********************************************************************
! abstract: use barnes or cressman analysis to obtain gridded field *
!           of analysis uncertainty                                 *
!                                                                   *
! program history log:                                              *
!   2005-10-08  pondeca                                             *
!                                                                   *
!********************************************************************
  use mpi
  use kinds, only: i_kind,r_single
  use cressanl_common, only: nobsmax, xlocs, ylocs, hgt0s, hgts, &
                             hobs, rmuses, oberrs, &
                             dtimes, cstations, & 
                             obstypes, dups, bckgs, xberrs, & 
                             terrain, jpointer, & 
                             kps,kts,kqs,kus,kvs,kugrds,kvgrds,kws,kw2s,kwds, & 
                             ktds,kgusts,kvis,kpblhs,kdists
  implicit none

  include 'param.incl'

! Declare passed variables
  integer(i_kind),intent(in):: mype,npe
  integer(i_kind),intent(in):: ista,iend,jsta,jend
  integer(i_kind),intent(in):: nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist
  real(r_single),intent(in):: rjbuffer_km
  real(r_single),intent(in):: ds0
  logical,intent(in)::lobjanl
  
! Declare local parameters
  real(r_single),parameter :: gravity=9.81
  integer(i_kind),parameter :: nflds=13   !order is ps,t,q,u,v,w,w2,wdir,td,gust,vis,pblh,dist
  integer(i_kind),parameter :: kflds=9    !perform cressman analysis  for 9 fields
                                    !(u,v,t,ps,q,gust,vis,pblh,dist)

  real(r_single),parameter::spval=-99999.
  character(10),parameter::chstar10='**********'
  character(1),parameter::cblank1=' '
 
! Declare local variables
  character(60) cgrid
  integer(i_kind) npass
  integer(i_kind) nx,ny
  integer(i_kind) i,j,m,n,nn,nt,np,k
  integer(i_kind) mlen,nlen
  integer(i_kind) lun1
  integer(i_kind) ierror,itype,jflag
  integer(i_kind) rtime(6),nlon,nlat,nsig
  integer(i_kind) ii,jj,nn0

  real(r_single) rlat,rlon,xx,yy,oberr,ob,ob_model,rinov
  real(r_single) qtflg,dtime,shgt0,hgt0,hgt,slm,rmuse
  real(r_single) rinovmax(kflds),radii(kflds)
  real(r_single) ds
  real(r_single) bb !dummy variable
  real(r_single) w0

  real(r_single) radius_u,rinovmax_u,radius_v,rinovmax_v, &
          radius_t,rinovmax_t,radius_p,rinovmax_p, &
          radius_q,rinovmax_q, & 
          radius_gust,rinovmax_gust,radius_vis,rinovmax_vis, & 
          radius_pblh,rinovmax_pblh,radius_dist,rinovmax_dist

  real(r_single),allocatable,dimension(:,:):: &
                 glon,glat,dx,dy,aux

  real(r_single),allocatable,dimension(:,:)::auxfield

  integer(i_kind),allocatable,dimension(:,:)::iauxfield

  character(8) cstation
  character(8) cprovider,csubprovider

  character(60) tdfname

  logical fexist1

  logical lpadjust
  logical usebckg

  character(2) clun1
  character(8) obstype

  character(60)  cvmodel(20) !names of verification models
  real(r_single)  rmusecv,rmuseb
  integer(i_kind) itotmodel,itotrmuse

  namelist/cress_barnes_anlqlty/cgrid,npass,&
                 radius_u,rinovmax_u,radius_v,rinovmax_v, &
                 radius_t,rinovmax_t,radius_p,rinovmax_p, &
                 radius_q,rinovmax_q, & 
                 radius_gust,rinovmax_gust,radius_vis,rinovmax_vis, &
                 radius_pblh,rinovmax_pblh,radius_dist,rinovmax_dist, & 
                 lpadjust,usebckg, &
                 cvmodel,rmusecv,rmuseb

!=====================================================================
! MPI setup
!  call mpi_init(ierror)
!  call mpi_comm_size(mpi_comm_world,npe,ierror)
!  call mpi_comm_rank(mpi_comm_world,mype,ierror)

   data cgrid/'conus'/
   data npass/1/
   data radius_u/18./ ; data rinovmax_u/10./
   data radius_v/18./ ; data rinovmax_v/10./
   data radius_t/12./ ; data rinovmax_t/7.5/
   data radius_p/11./ ; data rinovmax_p/1000./
   data radius_q/12./ ; data rinovmax_q/5.e-03/

   data radius_gust/18./ ; data rinovmax_gust/10./
   data radius_vis /18./ ; data rinovmax_vis/9000./
   data radius_pblh/18./ ; data rinovmax_pblh/1000./
   data radius_dist/18./ ; data rinovmax_dist/1000./

   data lpadjust/.false./
   data usebckg/.false./
!
!=====================================================================
   data rmuseb /+1./

   do n=1,20
      cvmodel(n)=chstar10
   enddo

   open (55,file='cress_barnesparm.anl',form='formatted')
   read(55,cress_barnes_anlqlty)
   close(55)

   if (mype==0) then
     print*,'========in anl_quality ==========='
     print*,'cgrid=',trim(cgrid)
     print*,'npass=',npass
     print*,'radius_u=',radius_u
     print*,'radius_v=',radius_v
     print*,'radius_t=',radius_t
     print*,'radius_p=',radius_p
     print*,'radius_q=',radius_q
     if (nrf2_gust > 0) print*,'radius_gust =',radius_gust
     if (nrf2_vis  > 0) print*,'radius_vis  =',radius_vis
     if (nrf2_pblh > 0) print*,'radius_pblh =',radius_pblh
     if (nrf2_dist > 0) print*,'radius_dist =',radius_dist

     print*,'in anl_quality:'

     print*,'rinovmax_u=',rinovmax_u
     print*,'rinovmax_v=',rinovmax_v
     print*,'rinovmax_t=',rinovmax_t
     print*,'rinovmax_p=',rinovmax_p
     print*,'rinovmax_q=',rinovmax_q
     if (nrf2_gust > 0) print*,'rinovmax_gust =',rinovmax_gust
     if (nrf2_vis  > 0) print*,'rinovmax_vis  =',rinovmax_vis
     if (nrf2_pblh > 0) print*,'rinovmax_pblh =',rinovmax_pblh
     if (nrf2_dist > 0) print*,'rinovmax_dist =',rinovmax_dist

     print*,'rmusecv=',rmusecv
     print*,'rmuseb=',rmuseb
     print*,'lpadjust=',lpadjust
     print*,'usebckg=',usebckg
   endif

   radii(1)=radius_u
   radii(2)=radius_v
   radii(3)=radius_t
   radii(4)=radius_p
   radii(5)=radius_q
   radii(6)=radius_gust
   radii(7)=radius_vis
   radii(8)=radius_pblh
   radii(9)=radius_dist

   if (usebckg) then
      rinovmax(1)=rinovmax_u
      rinovmax(2)=rinovmax_v
      rinovmax(3)=rinovmax_t
      rinovmax(4)=rinovmax_p
      rinovmax(5)=rinovmax_q
      rinovmax(6)=rinovmax_gust
      rinovmax(7)=rinovmax_vis
      rinovmax(8)=rinovmax_pblh
      rinovmax(9)=rinovmax_dist
     else
      print*,'in anl_quality: reset rinovmax and lpadjust & 
                 & to reflect non-use of first guess fields'
!     rinovmax(:)=1.e+20
      rinovmax(1)=rinovmax_u
      rinovmax(2)=rinovmax_v
      rinovmax(3)=rinovmax_t
      rinovmax(4)=1.e+20!        rinovmax_p
      rinovmax(5)=rinovmax_q
      rinovmax(6)=rinovmax_gust
      rinovmax(7)=rinovmax_vis
      rinovmax(8)=rinovmax_pblh
      rinovmax(9)=rinovmax_dist
      lpadjust=.false.
      if (mype==0) print*,'in anl_quality: lpadjust reset to: ',lpadjust
   endif

   call domain_dims(cgrid,nx,ny,ds)

   if (mype==0) print*,'in anl_quality: nx,ny=', & 
                                        nx,ny

   allocate(aux(nx,ny))
   allocate(auxfield(nx,ny))
   allocate(iauxfield(nx,ny))
!
!=====================================================================
!==>all tasks read in the first guess and the bckg error variances
!=====================================================================
   open (53,file='sigges',form='unformatted')

   read(53) rtime,nlon,nlat,nsig
   read(53) auxfield, aux        !glat,dx
   read(53) auxfield, aux        !glon,dy

   read(53) auxfield 
   do j=jsta,jend
   do i=ista,iend
      bckgs(i,j,4)=auxfield(i,j) !psfc1
   enddo
   enddo


   read(53) auxfield
   terrain=auxfield/gravity !no need for subdomain distribution

   read(53) auxfield
   do j=jsta,jend
   do i=ista,iend
      bckgs(i,j,3)=auxfield(i,j) !t1
   enddo
   enddo


   read(53) auxfield
   do j=jsta,jend
   do i=ista,iend
      bckgs(i,j,5)=auxfield(i,j) !q1
   enddo
   enddo


   read(53) auxfield
   do j=jsta,jend
   do i=ista,iend
      bckgs(i,j,1)=auxfield(i,j) !u1
   enddo
   enddo


   read(53) auxfield
   do j=jsta,jend
   do i=ista,iend
      bckgs(i,j,2)=auxfield(i,j) !v1
   enddo
   enddo


   do n=1,12                     !must jump 12 records to get to gust
      if (n==4 .or. n==5) then   !(see subroutine convert_binary_2d)
         read(53) iauxfield
        else
         read(53) auxfield
       endif
   enddo


   read(53) auxfield
   do j=jsta,jend
   do i=ista,iend
      bckgs(i,j,6)=auxfield(i,j) !gust1
   enddo
   enddo


   read(53) auxfield
   do j=jsta,jend
   do i=ista,iend
      bckgs(i,j,7)=auxfield(i,j) !vis1
   enddo
   enddo


   read(53) auxfield
   do j=jsta,jend
   do i=ista,iend
      bckgs(i,j,8)=auxfield(i,j) !pblh1
   enddo
   enddo


   do j=jsta,jend
   do i=ista,iend
      bckgs(i,j,9)=0.            !dist
   enddo
   enddo

   close(53)

   xberrs=0.
   do nn=1,kflds

      if (nn==6 .and. nrf2_gust <= 0) cycle
      if (nn==7 .and. nrf2_vis  <= 0) cycle
      if (nn==8 .and. nrf2_pblh <= 0) cycle
      if (nn==9 .and. nrf2_dist <= 0) cycle

      if (nn==1) open (94,file='bckgvar.dat_u', form='unformatted')
      if (nn==2) open (94,file='bckgvar.dat_v', form='unformatted')
      if (nn==3) open (94,file='bckgvar.dat_t', form='unformatted')
      if (nn==4) open (94,file='bckgvar.dat_ps',form='unformatted')
      if (nn==5) open (94,file='bckgvar.dat_q', form='unformatted')

      if (nn==6) open (94,file='bckgvar.dat_gust',form='unformatted')
      if (nn==7) open (94,file='bckgvar.dat_vis' ,form='unformatted')
      if (nn==8) open (94,file='bckgvar.dat_pblh',form='unformatted')
      if (nn==9) open (94,file='bckgvar.dat_dist',form='unformatted')

      read(94) auxfield

      do j=jsta,jend
      do i=ista,iend
         xberrs(i,j,nn)=auxfield(i,j)
      enddo
      enddo

      close(94)

      print*,'in anl_quality: mype,nn,xberrs,min,max=', & 
              mype,nn,minval(xberrs(:,:,nn)),maxval(xberrs(:,:,nn))
   enddo
!
!=====================================================================
!==>load array dups, which is used to remove duplicate obs in
!   the cressman analysis
!=====================================================================

   dups=.false.
   do nn=1,kflds
      if (nn==1) nn0=kugrds
      if (nn==2) nn0=kvgrds
      if (nn==3) nn0=kts
      if (nn==4) nn0=kps
      if (nn==5) nn0=kqs
      if (nn==6) nn0=kgusts
      if (nn==7) nn0=kvis
      if (nn==8) nn0=kpblhs
      if (nn==9) nn0=kdists


     do jj=1,nn0
        j=jpointer(jj,nn)
        nlen=0
        do k=1,8
           if (cstations(j)(k:k) /= cblank1) then
              nlen=nlen+1
            else
              exit
           endif
        enddo
 
        do ii=jj+1,nn0
           i=jpointer(ii,nn)
           mlen=0
           do k=1,8
              if (cstations(i)(k:k) /= cblank1) then
                 mlen=mlen+1
               else
                 exit
              endif
           enddo

           if ( mlen==nlen ) then
              if (cstations(i)(1:mlen)==cstations(j)(1:nlen)) then
                  if (abs(dtimes(i)) >= abs(dtimes(j))) dups(i)=.true.
                  if (abs(dtimes(i)) <  abs(dtimes(j))) dups(j)=.true.
              endif
           endif
        enddo
     enddo
   enddo

!=====================================================================
!==> perform cressman (or barnes) analysis
!=====================================================================
   if (lobjanl) & 
   call obj_anl4(rinovmax,lpadjust,usebckg,radii,spval, & 
                 nx,ny,kflds,nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist, & 
                 npass,ista,iend,jsta,jend,rjbuffer_km,ds0,npe,mype)

   if (mype==0) then
      open (54,file='sigcress',form='unformatted')
      open (53,file='sigges',form='unformatted')

      read(53) rtime,nlon,nlat,nsig
      write(54) rtime,nlon,nlat,nsig


      read(53) auxfield, aux        !glat,dx
      write(54) auxfield, aux

      read(53) auxfield, aux        !glon,dy
      write(54) auxfield, aux
 
      close(53)
   endif
  
   do n=1,6    
      if (n==1) nn=4   !pfsc
      if (n==3) nn=3   !temp
      if (n==4) nn=5   !moisture
      if (n==5) nn=1   !uwind
      if (n==6) nn=2   !vwind

      if (n/=2) then 
         aux(:,:)=0.
         do j=jsta,jend
         do i=ista,iend
            aux(i,j)=bckgs(i,j,nn)
         enddo
         enddo
         call mpi_allreduce(aux,auxfield,nx*ny, & 
                  mpi_real4,mpi_sum,mpi_comm_world,ierror)
         if (mype==0) write(54) auxfield
       else
         if (mype==0) write(54) terrain*gravity
       endif
   enddo

   aux(:,:)=0.
   do j=jsta,jend
   do i=ista,iend
      w0=bckgs(i,j,1)*bckgs(i,j,1) + & 
         bckgs(i,j,2)*bckgs(i,j,2)

      if (w0 > 0.) w0=sqrt(w0)

      if ( bckgs(i,j,6) < w0) then 
          aux(i,j)=w0
        else
          aux(i,j)=bckgs(i,j,6)
      endif
   enddo
   enddo
   call mpi_allreduce(aux,auxfield,nx*ny, &            !gust
            mpi_real4,mpi_sum,mpi_comm_world,ierror)

   aux(:,:)=0.
   if (mype==0) then
      do n=1,12
         write(54) aux       !must write 12 records before gust
      enddo                  !(see subroutine convert_binary_2d)
   endif

   if (mype==0) write(54) auxfield        !gust

   do n=1,2    
      if (n==1) nn=7   !vis
      if (n==2) nn=8   !pblh

      aux(:,:)=0.
      do j=jsta,jend
      do i=ista,iend
         aux(i,j)=bckgs(i,j,nn)
      enddo
      enddo
      call mpi_allreduce(aux,auxfield,nx*ny, & 
               mpi_real4,mpi_sum,mpi_comm_world,ierror)
      if (mype==0) write(54) auxfield
   enddo
     
!     aux(:,:)=0.
!     if (mype==0) write(54) aux     !dist   /no need for now
   close(54)

   tdfname='sigcress' 
   call td_flds(nlon,nlat,tdfname,mype,npe)

!==> verification
!    Determine if there is valid verf info:

     do n=1,20
       if (trim(cvmodel(n)) .ne. chstar10) then 
          if (mype==0) print*,' requested verf using model ',trim(cvmodel(n))
        else
          exit
      endif
     enddo
     itotmodel=n-1
     if (mype==0) print*,'itotmodel=',itotmodel

     if (itotmodel > 0) then
         do n=1,itotmodel
            call bias_rmse_cv_mpi(cgrid,ista,iend,jsta,jend,nx,ny, & 
                                  nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_dist, & 
                                  itotmodel,n,cvmodel(n),rmusecv,rmuseb, &
                                  lpadjust,mype,npe,cblank1)

         enddo
     endif

   deallocate(aux)
   deallocate(auxfield)
   deallocate(iauxfield)

end subroutine anl_quality
!------------------------------------------------------
       subroutine obj_anl4(rinovmax,lpadjust,usebckg,radii,spval, & 
                          nx,ny,kflds,igust,ivis,ipblh,idist, & 
                          npass,ista,iend,jsta,jend,rjbuffer_km,ds0,npe,mype)
       use mpi
       use kinds, only: i_kind,r_single
       use cressanl_common, only: nobsmax, xlocs, ylocs, &
                                  hobs, rmuses, oberrs, &
                                  obstypes, bckgs, xberrs, terrain

       implicit none  

       include 'param.incl'

       !==>Declare passed varaibles
       integer(i_kind),intent(in)  :: igust,ivis,ipblh,idist
       integer(i_kind),intent(in)  :: npe,mype
       integer(i_kind),intent(in)  :: ista,iend,jsta,jend
       integer(i_kind),intent(in)  :: nx,ny
       integer(i_kind),intent(in)  :: kflds,npass
       real(r_single),intent(in)   :: rjbuffer_km
       real(r_single),intent(in)   :: ds0
       real(r_single),intent(in)   :: rinovmax(kflds) 
       logical,intent(in)          :: lpadjust
       logical,intent(in)          :: usebckg
       real(r_single),intent(in)   :: radii(kflds)                       
       real(r_single),intent(in)   :: spval

       !==>Declare local parameters
       integer(i_kind),parameter:: ndomains=10  !npe  !ny/(20*nint(radii(nn)))

       !==>Declare local variables
       integer(i_kind) mype1
       integer(i_kind) jbuffer
       integer(i_kind) ierror
       integer(i_kind) iter,nn,ij,npts,npts2
       integer(i_kind) i,j,k,n,m,im0
       integer(i_kind) nobssubmax
       integer(i_kind) nclose
       integer(i_kind),allocatable,dimension(:,:):: nobssub
       integer(i_kind),allocatable,dimension(:,:,:):: npointer
       real(r_single),allocatable,dimension(:,:):: g0,g
       real(r_single),allocatable,dimension(:):: fincrs,qcflgs
       integer(i_kind),allocatable,dimension(:):: j1,j2


       integer(i_kind) nobsmax0,nobssubmax0,iproc
       integer(i_kind) ista0,iend0,jsta0,jend0
       integer(i_kind),allocatable,dimension(:,:):: nobssub0
       integer(i_kind),allocatable,dimension(:,:,:):: npointer0
       real(r_single),allocatable,dimension(:):: xlocs0,ylocs0,fincrs0,qcflgs0

       real(r_single) sumprod,wsum
       real(r_single) yleft,yright
       real(r_single) wij
       real(r_single) xob,yob,elev
       real(r_single) xi,xj,weight4,dist2,radius
       real(r_single) xberr0,oberr0
       real(r_single) u00,v00,t00,p00,q00,gust00,vis00,pblh00,dist00
!******************************************************************
!******************************************************************
       allocate(fincrs(nobsmax))
       allocate(qcflgs(nobsmax))
       allocate(nobssub(ndomains,kflds))
       allocate(j1(ndomains))
       allocate(j2(ndomains))

       mype1=mype+1

       jbuffer=1000._r_single*rjbuffer_km/ds0

       if (mype==0) print*,'in obj_anl4: rjbuffer_km,ds0,jbuffer=', &
                                         rjbuffer_km,ds0,jbuffer


       j1=-999 ; j2=-999

       call para_range3(1,ny,ndomains,j1(1:ndomains),j2(1:ndomains))

       nobssub(:,:)=0
       do n=1,nobsmax
          if     (trim(obstypes(n))=='ugrel-ob') then ;  nn=1
          elseif (trim(obstypes(n))=='vgrel-ob') then ;  nn=2
          elseif (trim(obstypes(n))=='t-ob')     then ;  nn=3
          elseif (trim(obstypes(n))=='p-ob')     then ;  nn=4
          elseif (trim(obstypes(n))=='q-ob')     then ;  nn=5
          elseif (trim(obstypes(n))=='gust-ob')  then ;  nn=6
          elseif (trim(obstypes(n))=='vis-ob')   then ;  nn=7
          elseif (trim(obstypes(n))=='pblh-ob')  then ;  nn=8
          elseif (trim(obstypes(n))=='dist-ob')  then ;  nn=9
          else                                        ; cycle ; endif

          do m=1,ndomains
             yleft  = float(max(1,(j1(m)-jbuffer)))
             yright = float(min(ny,(j2(m)+jbuffer)))

                if (ylocs(n) .ge. yleft .and. ylocs(n) .le. yright) then
                    nobssub(m,nn)=nobssub(m,nn)+1
                 endif
          enddo
       enddo

       do nn=1,kflds
          do m=1,ndomains
             print*,'in obj_anl4: mype,nn,m,nobssub(m,nn)=',mype,nn,m,nobssub(m,nn)
          enddo
       enddo
  
       nobssubmax=maxval(nobssub(:,:)) 

       print*,'in obj_anl4: mype,nobssubmax=',nobssubmax
       
       allocate(npointer(nobssubmax,ndomains,kflds))

       nobssub(:,:)=0
       npointer(:,:,:)=0

       do n=1,nobsmax
          if     (trim(obstypes(n))=='ugrel-ob') then ;  nn=1
          elseif (trim(obstypes(n))=='vgrel-ob') then ;  nn=2
          elseif (trim(obstypes(n))=='t-ob')     then ;  nn=3
          elseif (trim(obstypes(n))=='p-ob')     then ;  nn=4
          elseif (trim(obstypes(n))=='q-ob')     then ;  nn=5
          elseif (trim(obstypes(n))=='gust-ob')  then ;  nn=6
          elseif (trim(obstypes(n))=='vis-ob')   then ;  nn=7
          elseif (trim(obstypes(n))=='pblh-ob')  then ;  nn=8
          elseif (trim(obstypes(n))=='dist-ob')  then ;  nn=9
          else                                        ; cycle ; endif

          do m=1,ndomains
             yleft  = float(max(1,(j1(m)-jbuffer)))
             yright = float(min(ny,(j2(m)+jbuffer)))

                if (ylocs(n) .ge. yleft .and. ylocs(n) .le. yright) then
                    nobssub(m,nn)=nobssub(m,nn)+1
                    im0=nobssub(m,nn)
                    npointer(im0,m,nn)=n
                 endif
          enddo
       enddo

       if (.not.usebckg) then 
!           bckgs=0.                           !original  option
            call avghobs(mype,npe,u00,v00,t00,p00,q00, & 
                         gust00,vis00,pblh00,dist00)

            bckgs(:,:,1)=u00
            bckgs(:,:,2)=v00
            bckgs(:,:,3)=t00
            bckgs(:,:,4)=p00
            bckgs(:,:,5)=q00
            bckgs(:,:,6)=gust00
            bckgs(:,:,7)=vis00
            bckgs(:,:,8)=pblh00
            bckgs(:,:,9)=dist00
            if (mype==0) print*,'in obj_anl4: u00,v00,t00,p00,q00=',&
                                              u00,v00,t00,p00,q00
            if (mype==0) print*,'in obj_anl4: gust00,vis00,pblh00,dist00=',&
                                              gust00,vis00,pblh00,dist00
       endif

       do iter=1,npass
          call computeincrs(rinovmax,lpadjust, &
               fincrs,qcflgs,spval,igust,ivis,ipblh,idist, & 
               ista,iend,jsta,jend,nx,ny,kflds,nobsmax,mype,npe)

          do iproc=1,npe
             nobsmax0=nobsmax
             nobssubmax0=nobssubmax
             ista0=ista
             iend0=iend
             jsta0=jsta
             jend0=jend

             call mpi_bcast (nobsmax0,    1,mpi_integer,iproc-1,mpi_comm_world,ierror)
             call mpi_bcast (nobssubmax0, 1,mpi_integer,iproc-1,mpi_comm_world,ierror)
             call mpi_bcast (ista0,       1,mpi_integer,iproc-1,mpi_comm_world,ierror)
             call mpi_bcast (iend0,       1,mpi_integer,iproc-1,mpi_comm_world,ierror)
             call mpi_bcast (jsta0,       1,mpi_integer,iproc-1,mpi_comm_world,ierror)
             call mpi_bcast (jend0,       1,mpi_integer,iproc-1,mpi_comm_world,ierror)

             allocate(npointer0(nobssubmax0,ndomains,kflds))
             allocate(nobssub0(ndomains,kflds))
             allocate(xlocs0(nobsmax0))
             allocate(ylocs0(nobsmax0))
             allocate(fincrs0(nobsmax0))
             allocate(qcflgs0(nobsmax0))

             if (mype1==iproc) then
                npointer0(:,:,:)=npointer(:,:,:)
                nobssub0(:,:)=nobssub(:,:)
                xlocs0(:)=xlocs(:)
                ylocs0(:)=ylocs(:)
                fincrs0(:)=fincrs(:)
                qcflgs0(:)=qcflgs(:)
             endif
             call mpi_barrier(mpi_comm_world,ierror)

             npts=nobssubmax0*ndomains*kflds
             npts2=ndomains*kflds
             call mpi_bcast (npointer0, npts,     mpi_integer,  iproc-1, mpi_comm_world, ierror)
             call mpi_bcast (nobssub0,  npts2,    mpi_integer,  iproc-1, mpi_comm_world, ierror)
             call mpi_bcast (xlocs0,    nobsmax0, mpi_real,     iproc-1, mpi_comm_world, ierror)
             call mpi_bcast (ylocs0,    nobsmax0, mpi_real,     iproc-1, mpi_comm_world, ierror)
             call mpi_bcast (fincrs0,   nobsmax0, mpi_real,     iproc-1, mpi_comm_world, ierror)
             call mpi_bcast (qcflgs0,   nobsmax0, mpi_real,     iproc-1, mpi_comm_world, ierror)

             allocate(g (ista0:iend0,jsta0:jend0))
             allocate(g0(ista0:iend0,jsta0:jend0))
             
             do nn=1,kflds
                radius=radii(nn)
                g0=0.
                ij=0

                do m=1,ndomains

                   do j=j1(m), j2(m)     !these are the same on all processors
                      xj=float(j)

                      do i=ista0, iend0        !note these. they are the dims on the iproc processor
                         xi=float(i)

                         ij=ij+1
                         if (mype==mod(ij-1,npe)) then

                            sumprod=0.
                            wsum=0.
                            nclose=0

                            do k=1,nobssub0(m,nn)
                               n=npointer0(k,m,nn)

                               if (qcflgs0(n) .lt. 1.) cycle 

                               xob=xlocs0(n)
                               yob=ylocs0(n)

                               dist2=(xi-xob)*(xi-xob)+(xj-yob)*(xj-yob)
                               if (dist2 .le. (2.*radius)**2  ) then
                                  nclose=nclose+1
                                  wij=weight4(xi,xj,xob,yob,radius)
                                  sumprod=sumprod+wij*fincrs0(n)
                                  wsum=wsum+wij
                               endif
                            enddo !k-loop

                            !!!if (wsum.gt.1.e-5) then 
                            !!!   g0(i,j)=sumprod/wsum
                            !!!endif
                            if (nclose .ge. 10) then   !6
                               g0(i,j)=sumprod/wsum
                            endif

                         endif !mype-condition
                      enddo !i-loop
                   enddo !j-loop
                enddo !m-loop

                npts=(iend0-ista0+1)*(jend0-jsta0+1)
                call mpi_allreduce(g0(ista0:iend0,jsta0:jend0),g(ista0:iend0,jsta0:jend0), & 
                                   npts,mpi_real4,mpi_sum,mpi_comm_world,ierror)

                if (mype1==iproc) then
                   do j=jsta0,jend0
                      do i=ista0,iend0
                         bckgs(i,j,nn)=bckgs(i,j,nn)+g(i,j)
                      enddo
                   enddo
                endif
                call mpi_barrier(mpi_comm_world,ierror)

                if (mype==0) then
                   print*,'in obj_anl4: iter,iproc,nn,g,min,max=',iter,iproc,nn, & 
                           minval(g(ista0:iend0,jsta0:jend0)), & 
                           maxval(g(ista0:iend0,jsta0:jend0))
                endif
                call mpi_barrier(mpi_comm_world,ierror)

             enddo !nn-loop

             deallocate(npointer0)
             deallocate(nobssub0)
             deallocate(xlocs0)
             deallocate(ylocs0)
             deallocate(fincrs0)
             deallocate(qcflgs0)
             deallocate(g)
             deallocate(g0)

          enddo !iproc-loop
       enddo  !iter-loop

       deallocate(fincrs)
       deallocate(qcflgs)
       deallocate(nobssub)
       deallocate(npointer)
       deallocate(j1)
       deallocate(j2)

       return
       end

       function weight4(xi,xj,xob,yob,r)

       use kinds, only: r_single
       implicit none

       real(r_single) xi,xj,xob,yob,r,radius2,dist2,weight4

       radius2=r*r
       dist2=(xi-xob)*(xi-xob)+(xj-yob)*(xj-yob)
!      weight4=max(0.,(radius2-dist2)/(radius2+dist2))
       weight4=exp(-dist2/(2.*radius2))
!      weight4=exp(-4.*dist2/radius2) !Barnes type

       return
       end
!------------------------------------------------------
      subroutine avghobs(mype,npe,u00,v00,t00,p00,q00, & 
                         gust00,vis00,pblh00,dist00)

      use mpi
      use kinds, only: i_kind,r_single
      use cressanl_common, only: nobsmax, hobs, rmuses, &
                                 obstypes
      implicit none

      !==>Declare passed variables
      integer(i_kind),intent(in)  :: mype,npe
      real(r_single),intent(out)  :: u00,v00,t00,p00,q00, & 
                                     gust00,vis00,pblh00,dist00

      !Declare local parameters
      integer(i_kind),parameter::nt=9 !u,v,t,p,q,gust,vis,pblh,dist

      !Declare local variables
      integer(i_kind) n,nn,ierror
      integer(i_kind),allocatable::nokays(:),nokays2(:)
      real(r_single),allocatable::h00(:),h002(:)

      allocate(nokays(nt),nokays2(nt))
      allocate(h00(nt),h002(nt))

      u00=0. ; v00=0. ; t00=0. ; p00=0. ; q00=0.
      gust00=0. ; vis00=0. ; pblh00=0. ; dist00=0.
      h00(:)=0.
      nokays(:)=0

      do n=1,nobsmax
          if (abs(rmuses(n)-1.).ge.1.e-03) cycle

          if     (trim(obstypes(n))=='ugrel-ob') then ;  nn=1
          elseif (trim(obstypes(n))=='vgrel-ob') then ;  nn=2
          elseif (trim(obstypes(n))=='t-ob')     then ;  nn=3
          elseif (trim(obstypes(n))=='p-ob')     then ;  nn=4
          elseif (trim(obstypes(n))=='q-ob')     then ;  nn=5
          elseif (trim(obstypes(n))=='gust-ob')  then ;  nn=6
          elseif (trim(obstypes(n))=='vis-ob')   then ;  nn=7
          elseif (trim(obstypes(n))=='pblh-ob')  then ;  nn=8
          elseif (trim(obstypes(n))=='dist-ob')  then ;  nn=9
          else                                        ; cycle ; endif
          h00(nn)=h00(nn)+hobs(n)
          nokays(nn)=nokays(nn)+1
      enddo

      h002(:)=h00
      nokays2(:)=nokays(:)
      call mpi_allreduce(h002,    h00,    nt ,mpi_real4,    mpi_sum, mpi_comm_world,ierror)
      call mpi_allreduce(nokays2, nokays, nt ,mpi_integer, mpi_sum, mpi_comm_world,ierror)

      if (mype==0) then
         do nn=1,nt
            print*,'nn,nokays(nn)=',nn,nokays(nn)
         enddo
      endif

      if (nokays(1) > 0 ) u00=h00(1)/real(nokays(1),r_single)  
      if (nokays(2) > 0 ) v00=h00(2)/real(nokays(2),r_single)  
      if (nokays(3) > 0 ) t00=h00(3)/real(nokays(3),r_single)  
      if (nokays(4) > 0 ) p00=h00(4)/real(nokays(4),r_single)  
      if (nokays(5) > 0 ) q00=h00(5)/real(nokays(5),r_single)  

      if (nokays(6) > 0) gust00=h00(6)/real(nokays(6),r_single)  
      if (nokays(7) > 0) vis00= h00(7)/real(nokays(7),r_single)  
      if (nokays(8) > 0) pblh00=h00(8)/real(nokays(8),r_single)  
      if (nokays(9) > 0) dist00=h00(9)/real(nokays(9),r_single)  

      deallocate(nokays,nokays2)
      deallocate(h00,h002)

      end subroutine avghobs
!------------------------------------------------------
     subroutine computeincrs(rinovmax,lpadjust, &
               fincrs,qcflgs,spval,igust,ivis,ipblh,idist, &
               ista,iend,jsta,jend,nx,ny,kflds,nobsmax,mype,npe)

      use mpi
      use kinds, only: i_kind,r_single
      use cressanl_common, only: xlocs, ylocs, hgt0s, &
                                 hobs, rmuses, oberrs, &
                                 obstypes, insubdoms, dups, bckgs, & 
                                 terrain, jpointer, & 
                                 kts,kqs,kps,kugrds,kvgrds,kgusts, & 
                                 kvis,kpblhs,kdists
      implicit none

 
      !==>Declare passed variables
      integer(i_kind),intent(in)  :: igust,ivis,ipblh,idist
      integer(i_kind),intent(in)  :: ista,iend,jsta,jend
      integer(i_kind),intent(in)  :: nx,ny,kflds
      integer(i_kind),intent(in)  :: mype,npe
      integer(i_kind),intent(in)  :: nobsmax
      real(r_single),intent(in)   :: rinovmax(kflds)                       
      logical,intent(in)          :: lpadjust
      real(r_single),intent(in)   :: spval
      real(r_single),intent(out)  :: qcflgs(nobsmax)
      real(r_single),intent(out)  :: fincrs(nobsmax)

      !Declare local variables
      integer(i_kind) n,nn,nn0,i,j,ii
      integer(i_kind) ierror,mype1
      integer(i_kind) mts,mqs,mps,mus,mvs,mgusts,mvis,mpblhs,mdists
      integer(i_kind) mts2,mqs2,mps2,mus2,mvs2,mgusts2,mvis2,mpblhs2,mdists2
      real(r_single),allocatable,dimension(:,:):: field,field2
      real(r_single),allocatable,dimension(:,:):: bckgtv,tfield
      real(r_single) xx,yy,rmuse,fint
      real(r_single) tges,tob,qob,pges,zsges,hgt0
      real(r_single) rdelz,half,half_tlapse,g_over_rd,rdp,rtvts

      half=0.5
      half_tlapse=0.00325
      g_over_rd=9.81/287.04

      allocate(bckgtv(nx,ny))
      allocate(tfield(nx,ny))
      allocate(field(nx,ny))
      allocate(field2(nx,ny))
  
      field (:,:)=0.
      field2(:,:)=0.

      do j=jsta,jend
      do i=ista,iend
         field (i,j)=bckgs(i,j,3) ! temperature
         field2(i,j)=bckgs(i,j,3)*(1.+0.608*max(0.,bckgs(i,j,5))) !virtual temperature
      enddo
      enddo
      call mpi_allreduce(field,tfield,nx*ny, & 
                  mpi_real4,mpi_sum,mpi_comm_world,ierror)

      call mpi_allreduce(field2,bckgtv,nx*ny, & 
                  mpi_real4,mpi_sum,mpi_comm_world,ierror)

      if (mype==0) then
         print*,'bckg-t, min,max=',minval(tfield),maxval(tfield)
         print*,'bckg-tv, min,max=',minval(bckgtv),maxval(bckgtv)
      endif

      mype1=mype+1

      qcflgs(:)=0.
      fincrs(:)=0.

      mts=0 ; mqs=0 ; mps=0 ; mus=0 ; mvs=0 ; mgusts=0 ; mvis=0
      mpblhs=0 ; mdists=0

      do nn=1,kflds
         if (nn==1) nn0=kugrds
         if (nn==2) nn0=kvgrds
         if (nn==3) nn0=kts
         if (nn==4) nn0=kps
         if (nn==5) nn0=kqs
         if (nn==6) nn0=kgusts
         if (nn==7) nn0=kvis
         if (nn==8) nn0=kpblhs
         if (nn==9) nn0=kdists

         field2(:,:)=0.

         do j=jsta,jend
         do i=ista,iend
            field2(i,j)=bckgs(i,j,nn)
         enddo
         enddo

         call mpi_allreduce(field2,field,nx*ny, & 
                  mpi_real4,mpi_sum,mpi_comm_world,ierror)

         do ii=1,nn0
            n=jpointer(ii,nn)

            xx=xlocs(n)
            yy=ylocs(n)
            rmuse=rmuses(n)

            if (abs(rmuse-1.).gt.1.e-03 .or. dups(n)) cycle

            call bilinear_2d0v2(field,1,nx,1,ny,fint,yy,xx)!Note the reverse order "yy,xx"

            if (nn==4 .and. lpadjust) then
               pges=fint
               call bilinear_2d0v2(terrain,1,nx,1,ny,zsges,yy,xx)

               rdp=0.
               rtvts=1.               !hardwired to be sensible temperature for now
               if (rtvts==0.) then
                  call bilinear_2d0v2(bckgtv,1,nx,1,ny,tges,yy,xx)
                 else
                  call bilinear_2d0v2(tfield,1,nx,1,ny,tges,yy,xx)
               endif

               hgt0=hgt0s(n)
               tob=spval               !ideally would use observed temperature
               rdelz=hgt0-zsges

               if (tob /= spval) then
                  tges = half*(tges+tob)
                 else
                  if(rdelz < 0.)then
                     tges=tges-half_tlapse*rdelz
                  endif
               endif

               rdp = g_over_rd*rdelz/tges

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

!           if (nn==3) then
!              rtvts=1.               !hardwired to be sensible temperature for now
!              if (rtvts==0.) then
!                 call bilinear_2d0v2(bckgtv,1,nx,1,ny,fint,yy,xx)
!              endif
!           endif

            fincrs(n)=hobs(n)-fint
            if (abs(fincrs(n)).le.rinovmax(nn)) then 
                qcflgs(n)=+1.
                if (insubdoms(n)) then
                   if (nn==1) mus=mus+1
                   if (nn==2) mvs=mvs+1
                   if (nn==3) mts=mts+1
                   if (nn==4) mps=mps+1
                   if (nn==5) mqs=mqs+1
                   if (nn==6) mgusts=mgusts+1
                   if (nn==7) mvis=mvis+1 
                   if (nn==8) mpblhs=mpblhs+1 
                   if (nn==9) mdists=mdists+1 
                endif
            endif
         enddo
      enddo

      call mpi_allreduce(mus,    mus2,    1, mpi_integer, mpi_sum, mpi_comm_world, ierror)
      call mpi_allreduce(mvs,    mvs2,    1, mpi_integer, mpi_sum, mpi_comm_world, ierror)
      call mpi_allreduce(mts,    mts2,    1, mpi_integer, mpi_sum, mpi_comm_world, ierror)
      call mpi_allreduce(mps,    mps2,    1, mpi_integer, mpi_sum, mpi_comm_world, ierror)
      call mpi_allreduce(mqs,    mqs2,    1, mpi_integer, mpi_sum, mpi_comm_world, ierror)
      call mpi_allreduce(mgusts, mgusts2, 1, mpi_integer, mpi_sum, mpi_comm_world, ierror)
      call mpi_allreduce(mvis,   mvis2,   1, mpi_integer, mpi_sum, mpi_comm_world, ierror)
      call mpi_allreduce(mpblhs, mpblhs2, 1, mpi_integer, mpi_sum, mpi_comm_world, ierror)
      call mpi_allreduce(mdists, mdists2, 1, mpi_integer, mpi_sum, mpi_comm_world, ierror)

      if (mype==0) then
         print*,'in computeincrs: #of usable u-obs=',mus2
         print*,'in computeincrs: #of usable v-obs=',mvs2
         print*,'in computeincrs: #of usable t-obs=',mts2
         print*,'in computeincrs: #of usable p-obs=',mps2
         print*,'in computeincrs: #of usable q-obs=',mqs2
         print*,'in computeincrs: #of usable guts-obs=',mgusts2
         print*,'in computeincrs: #of usable vis-obs=',mvis2
         print*,'in computeincrs: #of usable pblh-obs=',mpblhs2
         print*,'in computeincrs: #of usable dist-obs=',mdists2
      endif

       deallocate(bckgtv)
       deallocate(tfield)
       deallocate(field)
       deallocate(field2)

       return
       end
!------------------------------------------------------
!------------------------------------------------------
       subroutine cvbasedrecal(cgrid,cmodel0,ista,iend,jsta,jend,nx,ny, & 
                  igust,ivis,ipblh,idist,mype,npe)

       use mpi
       use kinds, only: i_kind,r_single,r_kind
                   
       use errs_common, only: psierr,chierr,uerr,verr,uerr2,verr2, &
                              wspderr,wdirerr,wdirerr2,terr,tderr, &
                              qerr,perr,gusterr,viserr,pblherr,disterr


       implicit none

!declare passed variables
       character(60),intent(in)::cgrid
       character(60),intent(in)::cmodel0
       integer(i_kind),intent(in)::mype,npe
       integer(i_kind),intent(in)::ista,iend,jsta,jend
       integer(i_kind),intent(in)::nx,ny
       integer(i_kind),intent(in)::igust,ivis,ipblh,idist

!declare local parameters
       real(r_single),parameter::epsilon=1.e-06_r_single
       character(1),parameter::ccomma=','
       integer(i_kind),parameter::itotmodel=2
       logical,parameter::lpadjust=.true.

!declare local variables
       real(4) rmusecv,rmuseb
       character(60) cmodel
       integer(i_kind) nflds
       character(60) fnamesuffix
       character(80) rawstatsflname


       integer(i_kind),allocatable:: nobs(:)
       real(r_single),allocatable:: bias(:),abse(:),rmse(:)

       integer(i_kind),allocatable:: nobs2(:)
       real(r_single),allocatable:: bias2(:),abse2(:),rmse2(:)

       real(r_single),allocatable:: field(:,:)

       real(r_single) xx,yy,fint

       real(r_single) rmse0,rmse1,ratio
       logical lres1,lres2,lres3

       integer(i_kind) nobsmax0
       integer(i_kind) n,nn
       logical fexist


       print*,'in cvbasedrecal: cmodel0=',trim(cmodel0)

       if (trim(cmodel0)=='sigfupdate02') then !applicable when miter=2. cv-obs put back at start of 2nd outer loop
           fnamesuffix='_fg02only'  !this is the gsi output at the end of first outer loop
           rawstatsflname='rawstats.dat_for_rmuse_neg2'//trim(fnamesuffix)
           rmusecv=-2.
           rmuseb=+1.
         else if (trim(cmodel0)=='sigfupdate03') then !applicable when miter=3. cv-obs put back at start of 3rd outer loop
           fnamesuffix='_fg03only'  !this is the gsi output at the end of second outer loop
           rawstatsflname='rawstats.dat_for_rmuse_neg3'//trim(fnamesuffix)
           rmusecv=-3.
           rmuseb=+1.
        else
           fnamesuffix='_anlonly'  !applicable when miter=2, and cv-obs aren't put back into the analysis
           rawstatsflname='rawstats.dat_for_rmuse_neg3'//trim(fnamesuffix)
           rmusecv=-3.
           rmuseb=+1.
       endif

       print*,'in cvbasedrecal: fnamesuffix,rmusecv,rmuseb=',trim(fnamesuffix),rmusecv,rmuseb
       print*,'in cvbasedrecal: rawstatsflname=',trim(rawstatsflname)

       do n=1,itotmodel
          if (n==1) cmodel=cmodel0
          if (n==2) cmodel='errfield.dat_precval'
          call bias_rmse_cv_mpi(cgrid,ista,iend,jsta,jend,nx,ny,igust,ivis,ipblh,idist, & 
                                itotmodel,n,cmodel,rmusecv,rmuseb, &
                                lpadjust,mype,npe,fnamesuffix)
       enddo

       inquire(file=trim(rawstatsflname),exist=fexist)
       if (.not.fexist) then
          print*,'in cvbasedrecal: no input file available'
          print*,'apparently no cross-validation was performed to create input file'
        return
       endif

       open (18,file=trim(rawstatsflname),form='unformatted')

       read(18) cmodel !this is either sigfupdate02, sigfupdate03, or siganl
       read(18) nflds,nobsmax0

       allocate(nobs(nflds))
       allocate(bias(nflds))
       allocate(abse(nflds))
       allocate(rmse(nflds))

       read(18) nobs,bias,abse,rmse

       print*,'----------------------------'
       print*,'in cvbasedrecal: first read'
       print*,'cmodel=',trim(cmodel)
       print*,'nflds,nobsmax0=',nflds,nobsmax0
       do nn=1,nflds       !order is ps,t,q,u,v,w,w2,wdir,td,gust,vis,pblh,dist
          print*,'nn,nobs,bias,abse,rmse=',&
                  nn,nobs(nn),bias(nn),abse(nn),rmse(nn)
       enddo
       print*,'----------------------------'


       lres1=any(bias(1:9) > 9998.)   !handle gust & vis separately  !any(bias(1:nflds) > 9998.)
       lres2=any(abse(1:9) > 9998.)   !handle gust & vis separately  !any(abse(1:nflds) > 9998.)
       lres3=any(rmse(1:9) > 9998.)   !handle gust & vis separately  !any(rmse(1:nflds) > 9998.)

       if (lres1 .or. lres2 .or. lres3) then
          print*,'in cvbasedrecal/ first read : apparently no-calidation was done'
          print*,' do nothing and return'
          return
       endif

       open (19,file='cverrorupdate_input_current', form='formatted')

       write(19,'(a14)') '&cverrorupdate'
       write(19,'(a9,f14.6,a1)')  'pcverr = '     ,rmse(1),ccomma
       write(19,'(a9,f14.6,a1)')  'tcverr = '     ,rmse(2),ccomma
       write(19,'(a9,f14.6,a1)')  'qcverr = '     ,rmse(3),ccomma
       write(19,'(a9,f14.6,a1)')  'ucverr = '     ,rmse(4),ccomma
       write(19,'(a10,f14.6,a1)') 'u2cverr = '    ,rmse(4),ccomma
       write(19,'(a9,f14.6,a1)')  'vcverr = '     ,rmse(5),ccomma
       write(19,'(a10,f14.6,a1)') 'v2cverr = '    ,rmse(5),ccomma
       write(19,'(a12,f14.6,a1)') 'wspdcverr = '  ,rmse(6),ccomma
       write(19,'(a12,f14.6,a1)') 'wdircverr = '  ,rmse(8),ccomma
       write(19,'(a9,f14.6,a1)')  'tdcverr= '     ,rmse(9),ccomma
       do nn=10,11
          if (nobs(nn)==0 .or. bias(nn)>9998. .or. abse(nn)>9998. .or. rmse(nn)>9998.) then
             cycle
           else
             if (nn==10) write(19,'(a12,f14.6,a1)')  'gustcverr = '     ,rmse(10),ccomma
             if (nn==11) write(19,'(a11,f14.6,a1)')  'viscverr = '      ,rmse(11),ccomma
          endif
       enddo
       write(19,'(a1)') '/'
       close(19)

       read(18) cmodel  !this is errfield.dat_precval
       read(18) nflds,nobsmax0

       allocate(nobs2(nflds))   !order is ps,t,q,u,v,w,w2,wdir,td,gust,vis,pblh,dist
       allocate(bias2(nflds))
       allocate(abse2(nflds))
       allocate(rmse2(nflds))

       read(18) nobs2,bias2,abse2,rmse2


       print*,'----------------------------'
       print*,'in cvbasedrecal: second read'
       print*,'cmodel=',trim(cmodel)
       print*,'nflds,nobsmax0=',nflds,nobsmax0
       do nn=1,nflds
          print*,'nn,nobs2,bias2,abse2,rmse2=',&
                  nn,nobs2(nn),bias2(nn),abse2(nn),rmse2(nn)
       enddo
       print*,'----------------------------'


       lres1=any(bias2(1:9) > 9998.)   !handle gust & vis separately  !any(bias2(1:nflds) > 9998.)
       lres2=any(abse2(1:9) > 9998.)   !handle gust & vis separately  !any(abse2(1:nflds) > 9998.)
       lres3=any(rmse2(1:9) > 9998.)   !handle gust & vis separately  !any(rmse2(1:nflds) > 9998.)

       if (lres1 .or. lres2 .or. lres3) then
          print*,'in cvbasedrecal/ second read : apparently no-calidation was done'
          print*,' do nothing and return'
          return
       endif

       allocate(field(nx,ny))

       do nn=1,nflds
          if (nn >=12) cycle

          if (nobs(nn) <= 0) cycle

          if (nn==10 .or. nn==11) then 
              if (bias(nn)>9998. .or. abse(nn)>9998. .or. rmse(nn)>9998.) cycle
          endif

          rmse0=max(rmse2(nn),epsilon)
          rmse1=max(rmse(nn),epsilon)

          ratio=rmse1/rmse0            !!note: no sqrt needed
          if (nn==11) ratio=min(ratio,1.33) !avoid unacceptably large analysis errors for vis / 29Jun2012

          print*,'in cvbasedrecal:nn,nobs(nn),rmse1,rmse0,(rmse1/rmse0)=',nn,nobs(nn),rmse1,rmse0,ratio

          if (nn==1) perr=perr*ratio
          if (nn==2) terr=terr*ratio
          if (nn==3) qerr=qerr*ratio

          if (nn==4) then 
             uerr2=uerr2*ratio
             uerr=uerr*ratio                 !reasonable, I think
          endif

          if (nn==5) then 
             verr2=verr2*ratio
             verr=verr*ratio                  !reasonable, I think 
          endif

          if (nn==6) wspderr=wspderr*ratio

!         if (nn==8) then 
!            wdirerr2=wdirerr2*ratio
!            wdirerr=wdirerr*ratio             !reasonable, I think 
!         endif

          if (nn==9)  tderr=tderr*ratio
          if (nn==10) gusterr=gusterr*ratio
          if (nn==11) viserr=viserr*ratio
       enddo
                                                         !note: in future, consider disseminating uerr2,verr2,wdierr2 
                                                         !instead of uerr,verr,wdierr 
       close(18)

 
       deallocate(nobs,nobs2)
       deallocate(bias,bias2)
       deallocate(abse,abse2)
       deallocate(rmse,rmse2)
       deallocate(field)

       return
       end
!---------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------
       subroutine td_flds(nx,ny,fname,mype,npe)
       implicit none

       integer(4),intent(in):: nx,ny,mype,npe
       character(60),intent(in):: fname

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

       real(4),allocatable,dimension(:,:)::field1,field2
       real(4),allocatable,dimension(:,:)::psfc,t,q,td1,td2,td0

       real(4),parameter::qmin=1.e-06
       logical fexist

       allocate(field1(nx,ny))
       allocate(field2(nx,ny))
       allocate(psfc(nx,ny))
       allocate(t(nx,ny))
       allocate(q(nx,ny))
       allocate(td1(nx,ny))
       allocate(td2(nx,ny))
       allocate(td0(nx,ny))

       if (mype==0) print*,'in td_flds: fname=',trim(fname)

       if (trim(fname) /= 'sigges'       .and. & 
           trim(fname) /= 'siganl'       .and. & 
           trim(fname) /= 'sigfupdate02' .and. & 
           trim(fname) /= 'sigfupdate03' .and. & 
           trim(fname) /= 'sigcress'            ) then
           
           if (mype==0) print*,'invalid file name ... returning'

           return
       endif


       ioan1=51

!==> original first guess td. this is the field created by smartinit
       open (ioan1,file='slabs2.dat',form='unformatted')
       do n=1,4
         read(ioan1) td0
       enddo
       close(ioan1)
       if (mype==0) print*,'in td_flds: td0,min,max=',minval(td0),maxval(td0)
 
       if (mype==0) & 
          open (57,file='td.dat_'//trim(fname),form='unformatted')  !output file

       if (trim(fname) == 'sigges') then
          if (mype==0) then
             write(57) td0
             close(57)
          endif
          return
       endif

!==> alternative first guess td as derived from first guess psfc, q, and t.
!    it is used to compute the td increments, which are then added to td0 to
!    yield the total analysis td field

       open (ioan1,file='sigges',form='unformatted')

       read(ioan1) rtime,nlon,nlat,nsig
       read(ioan1) field1,field2 !glat,dx
       read(ioan1) field1,field2 !glon,dy

       if (mype==0) print*,'in td_flds: for sigges: rtime=',rtime
       if (mype==0) print*,'in td_flds: for sigges: nlon,nlat,nsig=',nlon,nlat,nsig

       read(ioan1) psfc
       read(ioan1) field1  !fis
       read(ioan1) t
       read(ioan1) q
       close(ioan1)

       call get_dewpt(psfc,q,t,td1,1,nx,1,ny)
       if (mype==0) print*,'in td_flds: td1,min,max=',minval(td1),maxval(td1)

!==> td for siganl, sigfupdate02, sigfupdate03, or sigcress
       inquire(file=trim(fname),exist=fexist)
       if (fexist) then
          open (ioan1,file=trim(fname),form='unformatted')

          read(ioan1) rtime,nlon,nlat,nsig
          read(ioan1) field1,field2 !glat,dx
          read(ioan1) field1,field2 !glon,dy

          if (mype==0) print*,'in td_flds: for ',trim(fname),': rtime=',rtime
          if (mype==0) print*,'in td_flds: for ',trim(fname),': nlon,nlat,nsig=',nlon,nlat,nsig

          read(ioan1) psfc
          read(ioan1) field1  !fis
          read(ioan1) t
          read(ioan1) q
          close(ioan1)

          q=max(qmin,q)
          call get_dewpt(psfc,q,t,td2,1,nx,1,ny)

          td2   = td2   - td1
          td2(:,:)=max(-20.,td2(:,:))
          call smther_one(td2,1,nx,1,ny,2)

          td2   = td2   + td0       !td1

          if (mype==0) then
             print*,'in td_flds: for ',trim(fname),': td2,min,max=',minval(td2),maxval(td2)
             write(57) td2
             close(57)
          endif
       endif

       deallocate(field1)
       deallocate(field2)
       deallocate(psfc)
       deallocate(t)
       deallocate(q)
       deallocate(td1)
       deallocate(td2)
       deallocate(td0)

       return
       end
!---------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------

       subroutine mp_flush(i)
       implicit none
       integer(4),intent(in) :: i
       return
       end
!---------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------
