Merge pull request #199 from jderber-NOAA/master

 GitHub Issue #175. Use the global 127L B-Matrix in regional FV3 DA (FV3LAMDA)

src/gsi/balmod.f90

@@ -396,17 +396,6 @@ subroutine prebal_reg(cwcoveqqcov)
 !   Set internal parameters to m_berror_stats
     call berror_set_reg('cwcoveqqcov',cwcoveqqcov)
 
-!   Read dimension of stats file
-    inerr=22
-    call berror_get_dims_reg(msig,mlat)
-
-!   Allocate arrays in stats file
-    allocate ( agvi(0:mlat+1,1:nsig,1:nsig) )
-    allocate ( bvi(0:mlat+1,1:nsig),wgvi(0:mlat+1,1:nsig) )
-
-!   Read in background error stats and interpolate in vertical to that specified in namelist
-    call berror_read_bal_reg(msig,mlat,agvi,bvi,wgvi,mype,inerr)
-
 !   ke_vp used to project SF to balanced VP
 !   below sigma level 0.8
 
@@ -419,12 +408,30 @@ subroutine prebal_reg(cwcoveqqcov)
        endif
     enddo j_loop
 
-    agvk=zero
-    bvk=zero
-    wgvk=zero
     ke_vp=ke-1
     if (twodvar_regional) ke_vp=ke
-    if (.not.twodvar_regional) then
+
+!   Read dimension of stats file
+    inerr=22
+    call berror_get_dims_reg(msig,mlat)
+
+!   Allocate arrays in stats file
+    allocate ( agvi(0:mlat+1,1:nsig,1:nsig) )
+    allocate ( bvi(0:mlat+1,1:nsig),wgvi(0:mlat+1,1:nsig) )
+
+!   Read in background error stats and interpolate in vertical to that specified in namelist
+    call berror_read_bal_reg(msig,mlat,agvi,bvi,wgvi,mype,inerr)
+
+!   Alternatively, zero out all balance correlation matrices
+!   for univariate surface analysis
+    if (twodvar_regional .or. lnobalance) then
+       if(mype==0) write(6,*)"***WARNING*** running univariate analysis." 
+       agvk=zero
+       bvk=zero
+       wgvk=zero
+       if(lnobalance) agvk_lm(:,:)=zero
+    else
+
        do k=1,ke_vp
           do j=1,lon2
              do i=1,lat2
@@ -471,20 +478,9 @@ subroutine prebal_reg(cwcoveqqcov)
           end do
        end do
     endif
-
-
-!   Alternatively, zero out all balance correlation matrices
-!   for univariate surface analysis
-    if (twodvar_regional .or. lnobalance) then
-       if(mype==0) write(6,*)"***WARNING*** running univariate analysis." 
-       bvk(:,:,:)=zero
-       agvk(:,:,:,:)=zero
-       wgvk(:,:,:)=zero
-       if(lnobalance) agvk_lm(:,:)=zero
-    endif
-
     deallocate (agvi,bvi,wgvi)
 
+
     return
   end subroutine prebal_reg
 
@@ -903,6 +899,7 @@ subroutine locatelat_reg(mype)
 !$$$
     use kinds, only: r_single
     use gridmod, only: nlon,nlat,lat2,lon2,istart,jstart,region_lat
+    use m_berror_stats, only: berror_stats
     use constants, only: deg2rad,one
     implicit none
 
@@ -918,10 +915,9 @@ subroutine locatelat_reg(mype)
 
 !   Read in dim of stats file
     lunin=22
-    open(lunin,file='berror_stats',form='unformatted')
+    open(lunin,file=berror_stats,form='unformatted')
     rewind lunin
     read(lunin)msig,mlat
-
 
 !   Allocate and read in lat array in stats file
     allocate( clat_avn(mlat), clat_avn4(mlat) )

src/gsi/glbsoi.f90

@@ -161,6 +161,7 @@ subroutine glbsoi
   use m_prad, only: prad_updatePredx    ! was -- prad_bias()
   use m_obsdiags, only: obsdiags_write
   use gsi_io,only: verbose
+  use m_berror_stats,only: inquire_berror
 
   implicit none
 
@@ -256,6 +257,8 @@ subroutine glbsoi
         end if
      end if
   else
+     lunit=22
+     call inquire_berror(lunit,mype)
      call create_balance_vars
      if(anisotropic) then
         call create_anberror_vars(mype)

src/gsi/gsi_rfv3io_mod.f90

@@ -653,14 +653,14 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin)
        call die('read_fv3_netcdf_guess','not enough PEs to read in fv3 fields' )
     endif
     mype_u=0           
-    mype_v=1
-    mype_t=2
-    mype_p=3
-    mype_q=4
-    mype_ql=5
-    mype_oz=6
-    mype_2d=7 
-    mype_delz=8
+    mype_v=mod(1,npe)
+    mype_t=mod(2,npe)
+    mype_p=mod(3,npe)
+    mype_q=mod(4,npe)
+    mype_ql=mod(5,npe)
+    mype_oz=mod(6,npe)
+    mype_2d=mod(7,npe)
+    mype_delz=mod(8,npe)
 
     allocate(ijns(npe),ijns2d(npe),ijnz(npe) )
     allocate(displss(npe),displss2d(npe),displsz_g(npe) )
@@ -978,7 +978,6 @@ subroutine gsi_fv3ncdf2d_read_v1(filenamein,varname,varname2,work_sub,mype_io)
     real(r_kind)   ,intent(out  ) :: work_sub(lat2,lon2) 
     integer(i_kind)   ,intent(in   ) :: mype_io
     real(r_kind),allocatable,dimension(:,:,:):: uu
-    integer(i_kind),allocatable,dimension(:):: dim_id,dim
     real(r_kind),allocatable,dimension(:):: work
     real(r_kind),allocatable,dimension(:,:):: a
 
@@ -1000,7 +999,6 @@ subroutine gsi_fv3ncdf2d_read_v1(filenamein,varname,varname2,work_sub,mype_io)
        endif
 
        iret=nf90_inquire(gfile_loc,ndimensions,nvariables,nattributes,unlimiteddimid)
-       allocate(dim(ndimensions))
        allocate(a(nlat,nlon))
 
        iret=nf90_inq_varid(gfile_loc,trim(adjustl(varname)),var_id)
@@ -1012,9 +1010,6 @@ subroutine gsi_fv3ncdf2d_read_v1(filenamein,varname,varname2,work_sub,mype_io)
        endif
 
        iret=nf90_inquire_variable(gfile_loc,var_id,ndims=ndim)
-       if(allocated(dim_id    )) deallocate(dim_id    )
-       allocate(dim_id(ndim))
-       iret=nf90_inquire_variable(gfile_loc,var_id,dimids=dim_id)
        if(allocated(uu       )) deallocate(uu       )
        allocate(uu(nx,ny,1))
        iret=nf90_get_var(gfile_loc,var_id,uu)
@@ -1030,7 +1025,7 @@ subroutine gsi_fv3ncdf2d_read_v1(filenamein,varname,varname2,work_sub,mype_io)
           end do
 
        iret=nf90_close(gfile_loc)
-       deallocate (uu,a,dim,dim_id)
+       deallocate (uu,a)
 
     endif !mype
 

src/gsi/gsimod.F90

@@ -59,7 +59,7 @@ module gsimod
                        l4dvar,nhr_obsbin,nhr_subwin,nwrvecs,iorthomax,&
                        lbicg,lsqrtb,lcongrad,lbfgsmin,ltlint,ladtest,ladtest_obs, lgrtest,&
                        idmodel,clean_4dvar,iwrtinc,lanczosave,jsiga,ltcost,liauon, &
-		       l4densvar,ens_nstarthr,lnested_loops,lwrite4danl,nhr_anal,thin4d,tau_fcst,efsoi_order
+                       l4densvar,ens_nstarthr,lnested_loops,lwrite4danl,nhr_anal,thin4d,tau_fcst,efsoi_order
   use gsi_4dvar, only: mPEs_observer
   use m_obsdiags, only: alwaysLocal => obsdiags_alwaysLocal
   use obs_ferrscale, only: lferrscale
@@ -103,6 +103,7 @@ module gsimod
   use derivsmod, only: init_anadv
   use berror, only: norh,ndeg,vs,bw,init_berror,hzscl,hswgt,pert_berr,pert_berr_fct,&
      bkgv_flowdep,bkgv_rewgtfct,bkgv_write,fpsproj,nhscrf,adjustozvar,fut2ps,cwcoveqqcov
+  use m_berror_stats, only: usenewgfsberror
   use anberror, only: anisotropic,ancovmdl,init_anberror,npass,ifilt_ord,triad4, &
      binom,normal,ngauss,rgauss,anhswgt,an_vs,&
      grid_ratio,grid_ratio_p,an_flen_u,an_flen_t,an_flen_z, &
@@ -768,7 +769,7 @@ module gsimod
 !     cwcoveqqcov  - sets cw Bcov to be the same as B-cov(q) (presently glb default)
 
   namelist/bkgerr/vs,nhscrf,hzscl,hswgt,norh,ndeg,noq,bw,norsp,fstat,pert_berr,pert_berr_fct, &
-	bkgv_flowdep,bkgv_rewgtfct,bkgv_write,fpsproj,adjustozvar,fut2ps,cwcoveqqcov
+      bkgv_flowdep,bkgv_rewgtfct,bkgv_write,fpsproj,adjustozvar,fut2ps,cwcoveqqcov,usenewgfsberror
 
 ! ANBKGERR (anisotropic background error related variables):
 !     anisotropic - if true, then use anisotropic background error

src/gsi/m_berror_stats.f90

@@ -43,9 +43,13 @@ module m_berror_stats
 
    private    ! except
 
-   ! reconfigurable parameters, via NAMELIST/setup/
-   public :: berror_stats    ! reconfigurable filename
+!   def usenewgfsberror - use modified gfs berror stats for global and regional.
+!                        for global skips extra record
+!                        for regional properly defines array sizes, etc.
+!   def berror_stats  - reconfigurable filename via NAMELIST/setup/
 
+   public :: usenewgfsberror
+   public :: berror_stats,inquire_berror
    ! interfaces to file berror_stats.
    public :: berror_get_dims ! get dimensions, jfunc::createj_func()
    public :: berror_read_bal ! get cross-cov.stats., balmod::prebal()
@@ -77,8 +81,10 @@ module m_berror_stats
    integer(i_kind),parameter :: default_unit_ = 22
    integer(i_kind),parameter :: default_rc_   = 2
 
-   character(len=256),save :: berror_stats = "berror_stats"      ! filename
    logical,save :: cwcoveqqcov_
+   logical usenewgfsberror
+
+   character(len=256),save :: berror_stats = "berror_stats"      ! filename
 
 contains
 
@@ -125,6 +131,52 @@ subroutine get_dims(msig,mlat,mlon,lunit)
    return
 end subroutine get_dims
 
+subroutine inquire_berror(lunit,mype)
+
+   use kinds,    only: r_single
+
+   implicit none
+
+   integer(i_kind),intent(in   ) :: lunit ! logical unit [22]
+   integer(i_kind),intent(in   ) :: mype
+
+   character(len=*),parameter :: myname_=myname//'::inquire_berror'
+   integer(i_kind) :: inerr,msig,mlat,mlon_,ier,errtot,i
+   real(r_single),dimension(:),allocatable::  clat_avn,sigma_avn
+
+   ! Read dimension of stats file
+   inerr = lunit
+   open(inerr,file=berror_stats,form='unformatted',status='old',iostat=ier)
+   call check_iostat(ier,myname_,'open('//trim(berror_stats)//')')
+   rewind inerr
+   read(inerr,iostat=ier) msig,mlat,mlon_
+   call check_iostat(ier,myname_,'read header')
+   errtot=ier
+
+   errtot=0
+   allocate ( clat_avn(mlat) )
+   allocate ( sigma_avn(1:msig) )
+   read(inerr,iostat=ier)clat_avn,sigma_avn
+!  Checking to see if sigma_avn fits required format if so newgfsberror file.
+   do i=1,msig-1
+     if(sigma_avn(i) < sigma_avn(i+1) .or. sigma_avn(i) < zero .or. sigma_avn(i) > one)then
+         errtot=1
+     end if
+   end do
+   deallocate(clat_avn,sigma_avn)
+   if(errtot /= 0)then
+     usenewgfsberror=.false.
+     if(mype == 0)write(6,*) 'usenewgfsberror = .false. old format file '
+   else
+     usenewgfsberror=.true.
+     if(mype == 0)write(6,*) 'usenewgfsberror = .true. new format file '
+   end if
+   close(inerr,iostat=ier)
+   call check_iostat(ier,myname_,'close('//trim(berror_stats)//')')
+
+   return
+end subroutine inquire_berror
+
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS  !
 !BOP -------------------------------------------------------------------
@@ -215,6 +267,8 @@ subroutine read_bal(agvin,bvin,wgvin,pputin,fut2ps,mype,lunit)
    read(inerr,iostat=ier) nsigstat,nlatstat
    call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (nsigstat,nlatstat)')
 
+!  dummy read to skip lats,sigma
+   if(usenewgfsberror)read(inerr)
    if ( mype==0 ) then
       if ( nsig/=nsigstat .or. nlat/=nlatstat ) then
          write(6,*) myname_,'(PREBAL):  ***ERROR*** resolution of ', &
@@ -261,7 +315,7 @@ subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwopt
 
    use kinds,only : r_single,r_kind
    use gridmod,only : nlat,nlon,nsig
-      use radiance_mod, only: n_clouds_fwd,cloud_names_fwd
+   use radiance_mod, only: n_clouds_fwd,cloud_names_fwd
 
    implicit none
 
@@ -309,13 +363,14 @@ subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwopt
 
    integer(i_kind) :: i,n,k,iq,icw,ivar,ic
    integer(i_kind) :: inerr,istat,ier
-   integer(i_kind) :: nsigstat,nlatstat
+   integer(i_kind) :: nsigstat,nlatstat,mlon_
    integer(i_kind) :: isig
    real(r_kind) :: corq2x
    character(len=5) :: var
    logical,allocatable,dimension(:) :: found3d
    logical,allocatable,dimension(:) :: found2d
 
+!  real(r_single),allocatable,dimension(:)  :: clat,sigma
    real(r_single),allocatable,dimension(:,:):: hwllin
    real(r_single),allocatable,dimension(:,:):: corzin
    real(r_single),allocatable,dimension(:,:):: corq2
@@ -331,8 +386,9 @@ subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwopt
    ! with that specified via the user namelist
 
    rewind inerr
-   read(inerr,iostat=ier)nsigstat,nlatstat
+   read(inerr,iostat=ier)nsigstat,nlatstat,mlon_
    call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (nsigstat,nlatstat)')
+   if(usenewgfsberror)read(inerr,iostat=ier)
 
    if ( mype==0 ) then
       if ( nsigstat/=nsig .or. nlatstat/=nlat ) then