module m_berror_stats
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:	 module m_berror_stats
!   prgmmr:	 j guo <jguo@nasa.gov>
!      org:	 NASA/GSFC, Global Modeling and Assimilation Office, 900.3
!     date:	 2010-03-24
!
! abstract:  a module of berror_stats input
!
! program history log:
!   2010-03-24  j guo   - added this document block
!   2011-08-01  lueken  - changed F90 to f90 (no machine logic) and fix indentation
!
!   input argument list: see Fortran 90 style document below
!
!   output argument list: see Fortran 90 style document below
!
! attributes:
!   language: Fortran 90 and/or above
!   machine:
!
!$$$  end subprogram documentation block

! module interface:

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS  !
!BOP -------------------------------------------------------------------
!
! !MODULE: m_berror_stats - a module of berror_stats input
!
! !DESCRIPTION:
!
! !INTERFACE:

      use kinds,only : i_kind
      use constants, only: one
      use control_vectors,only: cvars2d,cvars3d
      use mpeu_util,only: getindex
      use mpeu_util,only: perr,die

      implicit none

      private	! except

        ! reconfigurable parameters, via NAMELIST/setup/
      public :: berror_stats	! reconfigurable filename

        ! interfaces to file berror_stats.
      public :: berror_get_dims	! get dimensions, jfunc::createj_func()
      public :: berror_read_bal	! get cross-cov.stats., balmod::prebal()
      public :: berror_read_wgt	! get auto-cov.stats., prewgt()

      	! external interfaces relating to internal procedures.
      interface berror_get_dims; module procedure get_dims; end interface
      interface berror_read_bal; module procedure read_bal; end interface
      interface berror_read_wgt; module procedure read_wgt; end interface

! !REVISION HISTORY:
! 	30Jul08	- Jing Guo <guo@gmao.gsfc.nasa.gov>
!		- initial prototype/prolog/code to wrap up all file
!		  "berror_stats" related operations.
!       25Feb10 - Zhu
!               - made changes for generalizing control variables
!               - remove berror_nvars
!       14May13 - Jing Guo <jing.guo@nasa.gov>
!               - added I/O messages to aid run-time error diagnosis.
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname='m_berror_stats'
  character(len=256) :: berror_stats_chem = "berror_stats_chem" ! filename

  	! Reconfigurable parameters, vai NAMELISt/setup/
  character(len=256),save :: berror_stats = "berror_stats"	! filename

  integer(i_kind),parameter :: default_unit_ = 22
  integer(i_kind),parameter :: ERRCODE=2
contains
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS  !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: get_dims - get dimensions
!
! !DESCRIPTION:
!
! !INTERFACE:

    subroutine get_dims(msig,mlat,mlon,unit)

      implicit none

      integer(i_kind)         ,intent(  out) :: msig  ! dimension of levels
      integer(i_kind)         ,intent(  out) :: mlat  ! dimension of latitudes
      integer(i_kind),optional,intent(  out) :: mlon  ! dimension of longitudes
      integer(i_kind),optional,intent(in   ) :: unit  ! logical unit [22]

! !REVISION HISTORY:
! 	30Jul08	- Jing Guo <guo@gmao.gsfc.nasa.gov>
!		- the main body of the code is extracted from jfunc.f90
!       18Jun10 - todling - turn mlon into optional
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::get_dims'

  integer(i_kind) :: inerr,mlon_

! Read dimension of stats file
  inerr=default_unit_
  if(present(unit)) inerr = unit
  open(inerr,file=berror_stats,form='unformatted',status='old')
  rewind inerr
  read(inerr) msig,mlat,mlon_
  close(inerr)
  if(present(mlon))then
     mlon=mlon_
  endif
end subroutine get_dims
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS  !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: read_bal - get cross-corr. coefficients
!
! !DESCRIPTION:
!
! !INTERFACE:

    subroutine read_bal(agvin,bvin,wgvin,mype,unit)
      use kinds,only : r_single
      use gridmod,only : nlat,nlon,nsig

      implicit none

      real(r_single),dimension(nlat,nsig,nsig),intent(  out) :: agvin
      real(r_single),dimension(nlat,nsig)     ,intent(  out) :: bvin,wgvin
      integer(i_kind)                         ,intent(in   ) :: mype  ! "my" processor ID
      integer(i_kind),optional                ,intent(in   ) :: unit ! an alternative unit

! !REVISION HISTORY:
! 	30Jul08	- Jing Guo <guo@gmao.gsfc.nasa.gov>
!		- the main body of code for input is extracted from
!		  prebal() in balmod.f90.
!       25Feb10 - Zhu 
!               - change the structure of background error file
!               - read in agvin,wgvin,bvin only
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::read_bal'

!   workspaces/variables for data not returned

  integer(i_kind):: nsigstat,nlatstat
  integer(i_kind):: inerr,ier


!   Open background error statistics file
    inerr=default_unit_
    if(present(unit)) inerr=unit
    open(inerr,file=berror_stats,form='unformatted',status='old',iostat=ier)
    if(ier/=0) call die(myname_, &
       'open("'//trim(berror_stats)//'") error, iostat =',ier)

!   Read header.  Ensure that vertical resolution is consistent
!   with that specified via the user namelist

    rewind inerr
    read(inerr,iostat=ier) nsigstat,nlatstat
    if(ier/=0) call die(myname_, &
       'read("'//trim(berror_stats)//'") for (nsigstat,nlatstat) error, iostat =',ier)

    if(mype==0) then
       if (nsig/=nsigstat .or. nlat/=nlatstat) then
          write(6,*) myname_,'(PREBAL):  ***ERROR*** resolution of ', &
             '"',trim(berror_stats),'"', &
             'incompatiable with guess'
          write(6,*) myname_,'(PREBAL):  ***ERROR*** nsigstat,nlatstat=', &
             nsigstat,nlatstat
          write(6,*) myname_,'(PREBAL):  ***ERROR*** expects nsig,nlat=', &
             nsig,nlat
          call stop2(ERRCODE)
       end if

       write(6,*) myname_,'(PREBAL):  get balance variables', &
          '"',trim(berror_stats),'".  ', &
          'mype,nsigstat,nlatstat =', &
           mype,nsigstat,nlatstat
    end if

!   Read background error file to get balance variables
    read(inerr,iostat=ier) agvin,bvin,wgvin
    if(ier/=0) call die(myname_, &
       'read("'//trim(berror_stats)//'") for (agvin,bvin,wgvin) error, iostat =',ier)
    close(inerr,iostat=ier)
    if(ier/=0) call die(myname_, &
       'close("'//trim(berror_stats)//'") error, iostat =',ier)

    return
end subroutine read_bal
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS  !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: read_wgt - read auto-corr. coeffs.
!
! !DESCRIPTION:
!
! !INTERFACE:

    subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,mype,unit)

      use kinds,only : r_single,r_kind
      use gridmod,only : nlat,nlon,nsig
      use jfunc,only: varq,qoption
      use gridmod, only : raqms ! ajl
!     aod code
      use control_vectors,only: nrf
      use chemmod, only : berror_chem,upper2lower,lower2upper,&
           hzscl_fraction_chem,vs_fraction_chem,aero_ratios
      use gsi_chemguess_mod, only: gsi_chemguess_get
      use constants, only : max_varname_length
!     end aod code

      implicit none

      real(r_single),dimension(:,:,:),intent(inout) :: corz 
      real(r_single),dimension(:,:)  ,intent(inout) :: corp  
!     note test three all have a dimsion with nlat in them
      real(r_single),dimension(:,:,:),intent(inout) :: hwll
      real(r_single),dimension(:,:)  ,intent(inout) :: hwllp
      real(r_single),dimension(:,:,:),intent(inout) :: vz

      real(r_single),dimension(nlat,nlon),intent(out) :: corsst
      real(r_single),dimension(nlat,nlon),intent(out) :: hsst

      integer(i_kind)                    ,intent(in   ) :: mype  ! "my" processor ID
      integer(i_kind),optional           ,intent(in   ) :: unit ! an alternative unit
!      integer(i_kind) :: naer

! !REVISION HISTORY:
! 	30Jul08	- Jing Guo <guo@gmao.gsfc.nasa.gov>
!		- the main body of the code for input is extracted from
!		  prewgt() in prewgt.f90.
!       25Feb10  - Zhu - change the structure of background error file
!                      - make changes for generalizing control variables
!                      - move varq here from prewgt
!       28May10 - Todling - Obtain variable id's on the fly (add getindex) 
!                         - simpler logics to associate cv w/ berrors
!       14Jun10 - Todling - Allow any 3d berror not in file to be templated 
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::read_wgt'

!  workspace variables not returned
  real(r_single),dimension(nlat,nsig,nsig):: agvin
  real(r_single),dimension(nlat,nsig) :: wgvin,bvin
! aod add
  real(r_single),dimension(:),allocatable::  clat_avn,sigma_avn
  real(r_single),dimension(:,:),allocatable::  bv_avn,wgv_avn,corqq_avn
  real(r_single),dimension(:,:,:),allocatable:: agv_avn
  real(r_single),dimension(:,:),allocatable:: corz_avn,hwll_avn,vztdq_avn


  real(r_single),allocatable,dimension(:,:,:):: corz_tmp_chem
  real(r_single),allocatable,dimension(:,:,:):: hwll_tmp_chem
  real(r_single),allocatable,dimension(:,:,:):: vz_tmp_chem

  real(r_single),allocatable,dimension(:,:):: corz_tmp_chem_tot
  real(r_single),allocatable,dimension(:,:):: hwll_tmp_chem_tot
  real(r_single),allocatable,dimension(:,:):: vz_tmp_chem_tot

  character*5 :: varshort
  character(len=max_varname_length) :: var_chem
  logical,dimension(nrf):: nrf_err
  character(len=5), allocatable, dimension(:) :: cvars_gocart
  integer(i_kind) :: n_gocart_var
  logical :: vargocart,varaerotot
  real(r_single) :: fact_hs,fact_vs
  integer(i_kind) :: msig_chem,mlat_chem,j,m1


! end aod add
 
  integer(i_kind) :: i,n,k
  integer(i_kind) :: inerr,istat,ier
  integer(i_kind) :: nsigstat,nlatstat
  integer(i_kind) :: loc,nn,isig
  real(r_kind) :: corq2x
  character*5 var
  logical,allocatable,dimension(:) :: found3d
  logical,allocatable,dimension(:) :: found2d

  real(r_single),allocatable,dimension(:,:):: hwllin
  real(r_single),allocatable,dimension(:,:):: corzin
  real(r_single),allocatable,dimension(:,:):: corq2
  real(r_single),allocatable,dimension(:,:):: vscalesin

! Open background error statistics file
  inerr=default_unit_
  if(present(unit)) inerr=unit
!  write(6,*)mype,'file berror_stats',berror_stats
  open(inerr,file=berror_stats,form='unformatted',status='old',iostat=ier)
  if(ier/=0) call die(myname_, &
     'open("'//trim(berror_stats)//'") error, iostat =',ier)

! Read header.  Ensure that vertical resolution is consistent
! with that specified via the user namelist

  rewind inerr
  read(inerr,iostat=ier)nsigstat,nlatstat
  if(ier/=0) call die(myname_, &
     'read("'//trim(berror_stats)//'") for (nsigstat,nlatstat) error, iostat =',ier)

  if(mype==0) then
     if(nsigstat/=nsig .or. nlatstat/=nlat) then
        write(6,*)'PREBAL: **ERROR** resolution of berror_stats incompatiable with GSI'
        write(6,*)'PREBAL:  berror nsigstat,nlatstat=', nsigstat,nlatstat, &
           ' -vs- GSI nsig,nlat=',nsig,nlat
        call stop2(101)
     end if

     write(6,*) myname_,'(PREWGT):  read error amplitudes ', &
        '"',trim(berror_stats),'".  ', &
        'mype,nsigstat,nlatstat =', &
         mype,nsigstat,nlatstat
  end if
  read(inerr,iostat=ier) agvin,bvin,wgvin
  if(ier/=0) call die(myname_, &
     'read("'//trim(berror_stats)//'") for (agvin,bvin,wgvin) error, iostat =',ier)

! Read amplitudes
  allocate(found3d(size(cvars3d)),found2d(size(cvars2d)))
!  if(mype.eq.0)then
!     write(6,*)'cvars3d',cvars3d
!     write(6,*)'size',size(cvars3d)
!  endif
  found3d=.false.
  found2d=.false.
  read: do
     read(inerr,iostat=istat) var, isig
!     write(6,*)mype,'var',var,'isig',isig,'istat',istat
!     call flush(6)
     if (istat/=0) exit

     allocate ( corzin(nlat,isig) )
     if (var=='q') allocate ( corq2(nlat,isig) )
     allocate ( hwllin(nlat,isig) )
     if (isig>1) allocate ( vscalesin(nlat,isig) )

     if (var/='sst') then
        if (var=='q' .or. var=='Q') then
           read(inerr,iostat=ier) corzin,corq2
           if(ier/=0) call die(myname_, &
              'read("'//trim(berror_stats)//'") for (corzin,corq2) error, iostat =',ier)
        else
           read(inerr,iostat=ier) corzin
! ajl add for raqms
           if(raqms)then
             if (var == 'oz')then
!              note corzin is sqrt(variance) so in normal quantity units
               corzin=corzin*(47.98/28.97) ! convert vmr to mixing ration ! ajl
!             leave as ppv
!             elseif (var == 'co')then
!               corzin=corzin*1.e6 ! convert from vmr to ppmv ajl
             elseif (var == 'no2')then
               corzin=corzin*(46.0055/28.97) ! convert vmr to mixing ration ! ajl
!               write(6,*)'corzin',maxval(corzin),minval(corzin)
             endif
             if(var.eq.'co')then
!               do k=1,isig
                 do i=1,nlat
                   corzin(i,24)=max(3.e-9,corzin(i,24))
                   corzin(i,25)=max(3.e-9,corzin(i,25))
                   corzin(i,26)=max(3.6e-9,corzin(i,26))
                   corzin(i,27)=max(4.2e-9,corzin(i,27))
                   corzin(i,28)=max(4.8e-9,corzin(i,28))
                   corzin(i,29)=max(4.8e-9,corzin(i,29))
                   corzin(i,30)=max(5.e-9,corzin(i,30))
                   corzin(i,31)=max(5.2e-9,corzin(i,31))
                   corzin(i,32)=max(5.4e-9,corzin(i,32))
                   corzin(i,33)=max(6.4e-9,corzin(i,33))
                   corzin(i,34)=max(7.4e-9,corzin(i,34))
                   corzin(i,35)=max(11.e-9,corzin(i,35))
                 end do
!               end do
             endif
             if(mype.eq.0)then
               write(6,*)'var',var,nlat
               do k=1,isig
                 write(6,*)'corzin conus',k,maxval(corzin(116:146,k)),minval(corzin(116:146,k))
               end do
             endif
           endif
!            write(6,*)'corzin',maxval(corzin),minval(corzin)
! ajl end add for raqms
           if(ier/=0) call die(myname_, &
              'read("'//trim(berror_stats)//'") for (corzin) error, iostat =',ier)
        end if
        read(inerr,iostat=ier) hwllin
        if(ier/=0) call die(myname_, &
           'read("'//trim(berror_stats)//'") for (hwllin) error, iostat =',ier)
        if (isig>1) then
           read(inerr,iostat=ier) vscalesin
!           write(6,*)'ajl vscalesin',maxval(vscalesin),minval(vscalesin)
           if(ier/=0) call die(myname_, &
              'read("'//trim(berror_stats)//'") for (vscalesin) error, iostat =',ier)
        endif
     else
        read(inerr,iostat=ier) corsst
        if(ier/=0) call die(myname_, &
           'read("'//trim(berror_stats)//'") for (corsst) error, iostat =',ier)
        read(inerr,iostat=ier) hsst
        if(ier/=0) call die(myname_, &
           'read("'//trim(berror_stats)//'") for (hsst) error, iostat =',ier)
     end if
!     write(6,*)mype,'size corz',size(corz) ! ajl
     if (isig>1) then
!        write(6,*)'cvars3d',cvars3d
        n=getindex(cvars3d,var)
!        if(mype.eq.0)then
!         write(6,*)'n',n,'var',var,'cvars3d',cvars3d ! ajl
!        endif
        if(n>0)then
           found3d(n)=.true.
!          if(raqms)then
!          check if aersol           
!          endif
           do k=1,isig
              do i=1,nlat
                 corz(i,k,n)=corzin(i,k)
                 vz(k,i,n)=vscalesin(i,k)
              end do
              if(mype.eq.0)write(6,*)var,'vz k ',k,maxval(vz(k,:,n)),minval(vz(k,:,n)) !             ! ajl
              if(mype.eq.0)write(6,*)'corz k ',k,maxval(corz(:,k,n)),minval(corz(:,k,n)) !             ! ajl
           end do
           if (var=='q' .and. qoption==2)then
              do k=1,isig
                 do i=1,nlat
                    corq2x=corq2(i,k)
                    varq(i,k)=min(max(corq2x,0.00015_r_kind),one)
                 enddo
              enddo
              do k=1,isig
                 do i=1,nlat
                    corz(i,k,n)=one
                 end do
              end do
           end if
           do k=1,isig
              do i=1,nlat
                 hwll(i,k,n)=hwllin(i,k)
              end do
!              if(mype.eq.0)write(6,*)'hwll k ',k,maxval(hwll(k,:,:)),minval(hwll(k,:,:)) !             ! ajl
           end do
        endif ! n>0
     end if ! end of isig

     if (isig==1) then
        n= getindex(cvars2d,var)
        if (n>0.and.var/='sst') then
           found2d(n)=.true.
           do i=1,nlat
              corp(i,n)=corzin(i,1)
              hwllp(i,n)=hwllin(i,1)
           end do
        end if ! n>0
     end if ! isig=1

     deallocate(corzin,hwllin)
     if (isig>1) deallocate(vscalesin)
     if (var=='q') deallocate(corq2)
  enddo read 
  close(inerr)
! add aod
  !if(mype.eq.0)write(6,*)'ajl berror_chem',berror_chem
  IF(berror_chem) THEN

    CALL gsi_chemguess_get ( 'aerosols::3d', n_gocart_var, ier )
    IF (ier/=0) CALL stop2(260)
    ALLOCATE(cvars_gocart(n_gocart_var))
    CALL gsi_chemguess_get ('aerosols::3d',cvars_gocart,ier)
    IF (ier/=0) CALL stop2(261)

    fact_hs=hzscl_fraction_chem
    fact_vs=1.0_r_kind/vs_fraction_chem
    varaerotot=.false.

! open background error statistics file
    inerr=default_unit_
    if(present(unit)) inerr=unit
    !if(mype.eq.0)write(6,*)'ajl open',berror_stats_chem
    open(inerr,file=berror_stats_chem,form='unformatted',status='old')

! read header.
    rewind inerr
    read(inerr) nsigstat,nlatstat
    !write(6,*)'nsigstat',nsigstat,nlatstat,'nsig',nsig
    call flush(6)

    msig_chem=nsigstat
    mlat_chem=nlatstat

    IF (msig_chem/=nsig) THEN
      WRITE(6,*)'msig_chem=',msig_chem,' nsig=',nsig
      call stop2(262)
    ENDIF
    if (raqms )then
! **** ajl raqms


!      not sure non raqms code works correctly
!      raqms has input using nlat which equals nlat_chem  output arrays are
!      dimesioned nlat for global
       read(inerr,iostat=ier) agvin,bvin,wgvin
! read amplitudes

     ALLOCATE ( corz_avn(mlat_chem,1:msig_chem) )
     ALLOCATE ( hwll_avn(mlat_chem,1:msig_chem) )
     ALLOCATE ( vztdq_avn(1:msig_chem,mlat_chem) )
     read_chem_raqms: DO
        read(inerr,iostat=istat) varshort, isig
!        if(mype.eq.0)then
!        write(6,*)'varshort',varshort,isig
!        call flush(6)
!        endif

        IF (istat /= 0) EXIT
        
        var=upper2lower(varshort)
!        if(mype.eq.0)then
!        write(6,*)'var',var,'aero_ratios',aero_ratios
!        call flush(6)
!        endif

        IF (.NOT. aero_ratios) THEN

           vargocart=.FALSE.

           DO i=1,n_gocart_var
!        if(mype.eq.0)then
!              write(6,*)'var',var,'cvars_gocart',i,cvars_gocart(i)
!              call flush(6)
!        endif
              IF (var==cvars_gocart(i)) THEN 
                 vargocart=.TRUE.
                 exit
              ENDIF
           ENDDO

           IF (.NOT. vargocart) THEN 
              READ(inerr) 
              READ(inerr) 
              if (isig > 1) READ(inerr) 
              CYCLE read_chem_raqms
           ENDIF

!           IF(mype==0) THEN
!              WRITE(6,*)'Re-assigning bkg statistics for ',var
!           ENDIF

        ELSE

           vargocart=.FALSE.
           varaerotot=.FALSE.

           IF (varshort=='PMTOT') THEN
!this is a temporary name for name where total aero stats are written out
              IF(mype==0) THEN
                 WRITE(6,*)'Re-assigning and scaling bkg statistics for',var
              ENDIF
              varaerotot=.TRUE.
           ELSE
              DO i=1,n_gocart_var
                 IF (var==cvars_gocart(i)) THEN 
                    vargocart=.TRUE.
                    EXIT
                 ENDIF
              ENDDO
              
              IF (.NOT. vargocart) THEN 
                 READ(inerr) 
                 READ(inerr) 
                 IF (isig > 1) READ(inerr) 
                 CYCLE read_chem_raqms
              ENDIF
              
           ENDIF
        ENDIF

        IF (isig/=msig_chem) THEN
           IF (mype==0) THEN
              WRITE(6,*)'Incompatible berror files msig_chem/=isig',&
                   msig_chem,isig
           ENDIF
           CALL stop2(263)
        ENDIF
!        if(mype.eq.0)then
!        write(6,*)'read corz',shape(corz_avn)
!        call flush(6)
!        endif
      

        READ(inerr) corz_avn
!        if(mype.eq.0)then
!          write(6,*)'corz_avn',maxval(corz_avn),minval(corz_avn)
!        write(6,*)'read hwll',shape(hwll_avn)
!        call flush(6)
!        endif
        READ(inerr) hwll_avn
!        if(mype.eq.0)then
!          write(6,*)'hwll_avn',maxval(hwll_avn),minval(hwll_avn)
!        write(6,*)'read vzt',shape(vztdq_avn)
!        call flush(6)
!        endif
        READ(inerr) vztdq_avn


        n=getindex(cvars3d,var)
!        if(mype.eq.0)then
!          write(6,*)'vz ',maxval(vztdq_avn),minval(vztdq_avn)
!          write(6,*)'n',n,'cvars3d',cvars3d,'var',var
!          call flush(6)
!        endif
        
        IF (n > 0) THEN
!           need to convert variance to micro gram/kilogram just like wrf-chem
!           now have micrograms/kilogram
!            var=upper2lower(var)
!            naer=getindex(aeronames_lower,var)
            DO k=1,msig_chem
              DO j=1,mlat_chem
!                 corz(j,k,n)=corz_avn(j,k)*convert_raqms_microgram_kilogram(naer)
!                 fixed back ground file so don't need to convert
                 corz(j,k,n)=corz_avn(j,k)
              END DO
            END DO 
!            if(mype.eq.0)write(6,*)'corz eq',corz(nlat/2+1,:,n)

           IF (.NOT. aero_ratios) THEN
           
              DO k=1,msig_chem
                 DO j=1,mlat_chem
                    hwll(j,k,n)=hwll_avn(j,k)*fact_hs
                    vz(k,j,n)=vztdq_avn(k,j)*fact_vs
                 END DO
              END DO
!              if(mype.eq.0)then
!              write(6,*)'hwll eq',hwll(nlat/2+1,:,n)
!              write(6,*)'vz eq ',vz(:,nlat/2+1,n)
!              endif

           ENDIF
        endif
      end do read_chem_raqms
      deallocate (corz_avn,hwll_avn,vztdq_avn)
    else

     allocate ( clat_avn(mlat_chem) )
     allocate ( sigma_avn(1:msig_chem) )
     allocate ( agv_avn(0:mlat_chem+1,1:msig_chem,1:msig_chem) )
     allocate ( bv_avn(0:mlat_chem+1,1:msig_chem),wgv_avn(0:mlat_chem+1,1:msig_chem) )


     ALLOCATE(corz_tmp_chem(1:mlat_chem,msig_chem,nrf))
     ALLOCATE(hwll_tmp_chem(0:mlat_chem+1,msig_chem,nrf))
     ALLOCATE(vz_tmp_chem(msig_chem,0:mlat_chem+1,nrf))

     IF (aero_ratios) THEN
        ALLOCATE(corz_tmp_chem_tot(1:mlat_chem,msig_chem))
        ALLOCATE(hwll_tmp_chem_tot(0:mlat_chem+1,msig_chem))
        ALLOCATE(vz_tmp_chem_tot(msig_chem,0:mlat_chem+1))
     ENDIF

     ALLOCATE ( corz_avn(1:mlat_chem,1:msig_chem) )
     ALLOCATE ( hwll_avn(0:mlat_chem+1,1:msig_chem) )
     ALLOCATE ( vztdq_avn(1:msig_chem,0:mlat_chem+1) )

! read background error file to get balance variables
     if(.not.raqms)then
!      ajl raqms berror file does not have this just a dummy read
       read(inerr)clat_avn,(sigma_avn(k),k=1,msig_chem)
     endif
 
     call flush(6)
!     write(6,*)'read agv_avn',shape(agv_avn),shape(bv_avn),shape(wgv_avn)
     read(inerr)agv_avn,bv_avn,wgv_avn
!     write(6,*)'did read agv'
!     call flush(6)


! read amplitudes

     read_chem: DO
        write(6,*)'read varshoort'
        call flush(6)
        read(inerr,iostat=istat) varshort, isig
!        write(6,*)'varshort',varshort,isig
!        call flush(6)

        IF (istat /= 0) EXIT
        
        var=upper2lower(varshort)
!        write(6,*)'var',var,'aero_ratios',aero_ratios
!        call flush(6)

        IF (.NOT. aero_ratios) THEN

           vargocart=.FALSE.

           DO i=1,n_gocart_var
!              write(6,*)'var',var,'cvars_gocart',i,cvars_gocart(i)
!              call flush(6)
              IF (var==cvars_gocart(i)) THEN 
                 vargocart=.TRUE.
                 exit
              ENDIF
           ENDDO

           IF (.NOT. vargocart) THEN 
              READ(inerr) 
              READ(inerr) 
              if (isig > 1) READ(inerr) 
              CYCLE read_chem
           ENDIF

           IF(mype==0) THEN
              WRITE(6,*)'Re-assigning bkg statistics for ',var
           ENDIF

        ELSE

           vargocart=.FALSE.
           varaerotot=.FALSE.

           IF (varshort=='PMTOT') THEN
!this is a temporary name for name where total aero stats are written out
              IF(mype==0) THEN
                 WRITE(6,*)'Re-assigning and scaling bkg statistics for',var
              ENDIF
              varaerotot=.TRUE.
           ELSE
              DO i=1,n_gocart_var
                 IF (var==cvars_gocart(i)) THEN 
                    vargocart=.TRUE.
                    EXIT
                 ENDIF
              ENDDO
              
              IF (.NOT. vargocart) THEN 
                 READ(inerr) 
                 READ(inerr) 
                 IF (isig > 1) READ(inerr) 
                 CYCLE read_chem
              ENDIF
              
           ENDIF
        ENDIF

        IF (isig/=msig_chem) THEN
           IF (mype==0) THEN
              WRITE(6,*)'Incompatible berror files msig_chem/=isig',&
                   msig_chem,isig
           ENDIF
           CALL stop2(263)
        ENDIF
!        write(6,*)'read corz',shape(corz_avn)
!        call flush(6)
      

        READ(inerr) corz_avn
!        write(6,*)'read hwll',shape(hwll_avn)
!        call flush(6)
        READ(inerr) hwll_avn
!        write(6,*)'read vzt',shape(vztdq_avn)
!        call flush(6)
        READ(inerr) vztdq_avn


        n=getindex(cvars3d,var)
        if(mype.eq.0)then
        write(6,*)'n',n,'cvars3d',cvars3d,'var',var
        call flush(6)
        endif
        
        IF (n > 0) THEN
!          ajl fixed units so dont need to convert already micro gram/kilogram
            DO k=1,msig_chem
              DO j=1,mlat_chem
                 corz_tmp_chem(j,k,n)=corz_avn(j,k)
              END DO
            END DO

           IF (.NOT. aero_ratios) THEN
           
              DO k=1,msig_chem
                 DO j=0,mlat_chem+1
                    hwll_tmp_chem(j,k,n)=hwll_avn(j,k)*fact_hs
                    vz_tmp_chem(k,j,n)=vztdq_avn(k,j)*fact_vs
                 END DO
              END DO

           ENDIF

        ELSE IF (varaerotot) THEN

           DO k=1,msig_chem
              DO j=1,mlat_chem
                 corz_tmp_chem_tot(j,k)=corz_avn(j,k)
              END DO
           END DO

           DO k=1,msig_chem
              DO j=0,mlat_chem+1
                 hwll_tmp_chem_tot(j,k)=hwll_avn(j,k)*fact_hs
                 vz_tmp_chem_tot(k,j)=vztdq_avn(k,j)*fact_vs
              END DO
           END DO

        ELSE
           
           IF (mype==0) WRITE(6,*)'Should not be here, stopping'
           
           CALL stop2(264)           

        ENDIF

     enddo read_chem

     close(inerr)

!assign total variance whihc is scaled in prewg_reg.f90 based on
!mass ratio

     IF (aero_ratios) THEN

        DO k=1,msig_chem

           DO j=1,mlat_chem

              DO i=1,n_gocart_var
                 var=cvars_gocart(i)
                 n=getindex(cvars3d,var)
                 corz_tmp_chem(j,k,n)=corz_tmp_chem_tot(j,k)
              END DO
           ENDDO

        END DO

     ENDIF

     m1=MIN(mlat_chem/2+1,mlat_chem)

     DO i=1,n_gocart_var

        var=cvars_gocart(i)

        n=getindex(cvars3d,var)

        IF (n < 1) CALL stop2(265)

        IF (.NOT. aero_ratios) THEN

           DO k=1,nsig

!              DO j=0,mlat+1
              DO j=0,nlat+1
                 hwll(j,k,n)=hwll_tmp_chem(m1,k,n)
                 vz(k,j,n)=vz_tmp_chem(k,m1,n)
              END DO
              
!              DO j=1,mlat
              DO j=1,nlat
                 corz(j,k,n)=corz_tmp_chem(m1,k,n)
              END DO

           ENDDO

        ELSE

           DO k=1,nsig
!              DO j=0,mlat+1
              DO j=0,nlat+1
                 hwll(j,k,n)=hwll_tmp_chem_tot(m1,k)
                 vz(k,j,n)=vz_tmp_chem_tot(k,m1)
              END DO
!              DO j=1,mlat
              DO j=1,nlat
                 corz(j,k,n)=corz_tmp_chem(m1,k,n)
              END DO
           ENDDO
           
        ENDIF
        
     ENDDO

     deallocate(clat_avn,sigma_avn)
     deallocate(agv_avn,bv_avn,wgv_avn)

     DEALLOCATE(corz_tmp_chem,hwll_tmp_chem,vz_tmp_chem)

     IF (aero_ratios) DEALLOCATE(hwll_tmp_chem_tot,vz_tmp_chem_tot,&
          corz_tmp_chem_tot)
   endif ! ajl raqms
  ENDIF
  IF (berror_chem)  DEALLOCATE(cvars_gocart)
! end add aod

! corz, hwll & vz for undefined 3d variables
  do n=1,size(cvars3d)
     if(.not.found3d(n)) then
        if(n>0) then
           write(6,*)'not found3d cvars3d',n,cvars3d(n)
           if(cvars3d(n)=='oz') then

              call setcoroz_(corz(:,:,n),mype)
           else
              call setcorchem_(cvars3d(n),corz(:,:,n),ier)
              if(ier/=0) cycle ! if this happens, code will crash later
           endif
           call sethwlloz_(hwll(:,:,n),mype)
           call setvscalesoz_(vz(:,:,n))
        endif
        if(mype==0) write(6,*) myname_, ': WARNING, using general Berror template for ', cvars3d(n)
     endif
  enddo

! need simliar general template for undefined 2d variables ...

  deallocate(found3d,found2d)

  return
end subroutine read_wgt

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS  !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: setcoroz_ - a modeled corr.coeffs. of ozone
!
! !DESCRIPTION:
!
! !INTERFACE:

    subroutine setcoroz_(coroz,mype)
      use kinds,    only: r_single,r_kind
      use constants,only: zero,rozcon,one
      use mpimod,   only: npe,mpi_rtype,mpi_sum,mpi_comm_world

      use gridmod,  only: nlat,nsig
      use gridmod,  only: lon1,lat1

      use guess_grids,only: ntguessig
      use guess_grids,only: ges_oz   ! ozone fields
      use guess_grids,only: ges_prsi ! interface pressures (kPa)

      implicit none

      real(r_single),dimension(nlat,nsig),intent(  out) :: coroz ! of ozone
      integer(i_kind)              ,intent(in   ) :: mype     ! ID of this processor

! !REVISION HISTORY:
! 	31Jul08	- Jing Guo <guo@gmao.gsfc.nasa.gov>
!		- adopted from PREWGT of previous version
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::setcoroz_'
  real(r_kind),parameter:: r25 = one/25.0_r_kind

!! -- workspace and working variables

    real(r_kind),dimension(nsig+1,npe) :: work_oz,work_oz1
    real(r_kind),dimension(nsig) :: ozmz
    real(r_kind) :: asum,bsum

    integer(i_kind) :: mlat,msig
    integer(i_kind) :: i,j,k,n,mm1
    integer(i_kind) :: ierror

!! -- synity check
    if(mype==0) then
       write(6,*) myname_,'(PREWGT): mype = ',mype
    endif

    mlat=size(coroz,1)
    msig=size(coroz,2)
    if(mlat/=nlat .or. msig/=nsig) then
       write(6,*) myname_,'(PREWGT): shape mismatching on PE ',mype
       write(6,*) myname_,'(PREWGT): shape(coroz) = ',shape(coroz)
       write(6,*) myname_,'(PREWGT): while expecting nlat = ',nlat
       write(6,*) myname_,'(PREWGT): while expecting nsig = ',nsig
       call stop2(ERRCODE)
    endif

!! -- The first part is taken from read_guess().

! Calculate global means for ozone
! Calculate sums for ozone to estimate variance.
  mm1=mype+1
  work_oz = zero
  do k = 1,nsig
     do j = 2,lon1+1
        do i = 2,lat1+1
           work_oz(k,mm1) = work_oz(k,mm1) + ges_oz(i,j,k,ntguessig)* &
              rozcon*(ges_prsi(i,j,k,ntguessig)-ges_prsi(i,j,k+1,ntguessig))
        end do
     end do
  end do
  work_oz(nsig+1,mm1)=float(lon1*lat1)

  call mpi_allreduce(work_oz,work_oz1,(nsig+1)*npe,mpi_rtype,mpi_sum,&
       mpi_comm_world,ierror)
  if(ierror/=0) then
     write(6,*) myname_,'(PREWGT): MPI_allreduce() error on PE ',mype
     call stop2(ierror)
  endif

!! -- All it does above, through mm1 plus mpi_allreduce() to work_oz1[],
!! seems no more than a mpi_allgatherv() to me.  The "reduce" part is
!! actually done below ...

  bsum=zero
  do n=1,npe
     bsum=bsum+work_oz1(nsig+1,n)
  end do
  do k=1,nsig
     ozmz(k)=zero
     asum=zero
     do n=1,npe
        asum=asum+work_oz1(k,n)
     end do
     if (bsum>zero) ozmz(k)=asum/bsum
  enddo

!! -- now this part is taken from prewgt().

!   load variances onto subdomains
  do k=1,nsig
     coroz(:,k) = max(ozmz(k),0.0002_r_kind)*r25
  enddo

end subroutine setcoroz_
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS  !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: sethwlloz_ - a modeled hwll of ozone
!
! !DESCRIPTION:
!
! !INTERFACE:

    subroutine sethwlloz_(hwlloz,mype)
      use kinds,   only: r_single,r_kind
      use mpimod,  only: levs_id
      use gridmod, only: nnnn1o,nsig,nlon,nlat
      use constants,only: two,three,pi,rearth_equator
      implicit none

      real(r_single),dimension(nlat,nsig),intent(  out) :: hwlloz
      integer(i_kind)              ,intent(in   ) :: mype ! ID of this processor

! !REVISION HISTORY:
! 	31Jul08	- Jing Guo <guo@gmao.gsfc.nasa.gov>
!		- initial prototype/prolog/code
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::sethwlloz_'

  real(r_kind),parameter :: r400=400._r_kind
  real(r_kind),parameter :: r800=800._r_kind
  real(r_kind),parameter :: r40000=40000._r_kind

  integer(i_kind) :: k,k1
  real(r_kind) :: fact
  real(r_kind) :: s2u
    
  if(mype==0) then
     write(6,*) myname_,'(PREWGT): mype = ',mype
  endif

  s2u=(two*pi*rearth_equator)/nlon
  do k=1,nnnn1o
     k1=levs_id(k)
     if(k1>0) then
     write(6,*) myname_,'(PREWGT): mype = ',mype, k1
        if(k1<=nsig*3/4)then
        !  fact=1./hwl
           fact=r40000/(r400*nlon)
        else
           fact=r40000/(nlon*(r800-r400*(nsig-k1)/(nsig-nsig*3/4)))
        endif
        fact=fact*three
        hwlloz(:,k1)=s2u/fact
     endif
  enddo


  if(mype==0) then
     write(6,*) myname_,'(PREWGT): mype = ',mype, 'finish sethwlloz_'
  endif


end subroutine sethwlloz_
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS  !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: setvscalesoz_ - a modeled vscales for ozone
!
! !DESCRIPTION:
!
! !INTERFACE:

    subroutine setvscalesoz_(vscalesoz)
      use gridmod,only : nlat,nlon,nsig
      use kinds,only: r_single,r_kind
      implicit none

      real(r_single),dimension(nsig,nlat),intent(  out) :: vscalesoz

! !REVISION HISTORY:
! 	31Jul08	- Jing Guo <guo@gmao.gsfc.nasa.gov>
!		- initial prototype/prolog/code
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::setvscalesoz_'
  real(r_kind),parameter:: eight_tenths = 0.8_r_kind

  	! a fixed value is used.
  vscalesoz(:,:)=eight_tenths

end subroutine setvscalesoz_

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS  !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: setcorchem_ - a modeled corr.coeffs. of chemistry
!
! !DESCRIPTION:
!
! !INTERFACE:

    subroutine setcorchem_(cname,corchem,rc)
      use kinds,    only: r_single,r_kind
      use mpimod,   only: mype
      use constants,only: zero,one
      use mpimod,   only: npe,mpi_rtype,mpi_sum,mpi_comm_world

      use gridmod,  only: nlat,nsig
      use gridmod,  only: lon1,lat1

      use guess_grids,only: ntguessig
      use guess_grids,only: ges_prsi ! interface pressures (kPa)

      use gsi_chemguess_mod, only: gsi_chemguess_bundle
      use gsi_bundlemod,     only: gsi_bundlegetpointer

      implicit none

      character(len=*)                   ,intent(in   ) :: cname   ! constituent name
      real(r_single),dimension(nlat,nsig),intent(  out) :: corchem ! constituent correlations
      integer(i_kind)                    ,intent(  out) :: rc      ! return error code

! !REVISION HISTORY:
!    15Jul20010 - Todling - created from Guo's OZ routine
!
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::setcorchem_'
  real(r_kind),parameter:: r25 = one/25.0_r_kind

!! -- workspace and working variables

    real(r_kind),dimension(nsig+1,npe) :: work_chem,work_chem1
    real(r_kind),dimension(nsig) :: chemz
    real(r_kind) :: asum,bsum

    integer(i_kind) :: mlat,msig
    integer(i_kind) :: i,j,k,n,iptr,mm1
    integer(i_kind) :: ierror

    rc=0

!! -- synity check
    if(mype==0) then
       write(6,*) myname_,'(PREWGT): mype = ',mype
    endif

!   Get information for how to use CO2
    iptr=-1
    if(size(gsi_chemguess_bundle)>0) then ! check to see if bundle's allocated
       call gsi_bundlegetpointer(gsi_chemguess_bundle(1),cname,iptr,ierror)
       if(ierror/=0)then
          rc=-2  ! field not found
          return 
       endif
    else
       rc=-1     ! chem not allocated
       return
    endif

    mlat=size(corchem,1)
    msig=size(corchem,2)
    if(mlat/=nlat .or. msig/=nsig) then
       write(6,*) myname_,'(PREWGT): shape mismatching on PE ',mype
       write(6,*) myname_,'(PREWGT): shape(corchem',trim(cname),') = ',shape(corchem)
       write(6,*) myname_,'(PREWGT): while expecting nlat = ',nlat
       write(6,*) myname_,'(PREWGT): while expecting nsig = ',nsig
       call stop2(ERRCODE)
    endif

! -- The first part is taken from read_guess().

!   Calculate global means for constituent
!   Calculate sums for constituent to estimate variance.
    mm1=mype+1
    work_chem = zero
    do k = 1,nsig
       do j = 2,lon1+1
          do i = 2,lat1+1
             work_chem(k,mm1) = work_chem(k,mm1) + gsi_chemguess_bundle(ntguessig)%r3(iptr)%q(i,j,k)* &
                (ges_prsi(i,j,k,ntguessig)-ges_prsi(i,j,k+1,ntguessig))
!_RT not sure yet how to handle scaling factor (rozcon) in general
!_RT            rozcon*(ges_prsi(i,j,k,ntguessig)-ges_prsi(i,j,k+1,ntguessig))
          end do
       end do
    end do
    work_chem(nsig+1,mm1)=float(lon1*lat1)
  
    call mpi_allreduce(work_chem,work_chem1,(nsig+1)*npe,mpi_rtype,mpi_sum,&
         mpi_comm_world,ierror)
    if(ierror/=0) then
       write(6,*) myname_,'(PREWGT): MPI_allreduce() error on PE ',mype
       call stop2(ierror)
    endif

!   -- All it does above, through mm1 plus mpi_allreduce() to work_chem1[],
!   seems no more than a mpi_allgatherv() to me.  The "reduce" part is
!   actually done below ...

    bsum=zero
    do n=1,npe
       bsum=bsum+work_chem1(nsig+1,n)
    end do
    do k=1,nsig
       chemz(k)=zero
       asum=zero
       do n=1,npe
          asum=asum+work_chem1(k,n)
       end do
       if (bsum>zero) chemz(k)=asum/bsum
    enddo

! -- now this part is taken from prewgt().

!   load variances onto subdomains
    do k=1,nsig
       corchem(:,k) = max(chemz(k),0.0002_r_kind)*r25
    enddo

    if (mype==0) then
       write(6,*) myname_, ': Defined general B cov for: ', trim(cname)
    endif

end subroutine setcorchem_
end module m_berror_stats
