From b846b29ef16136828b1e08504f45bb1f05398706 Mon Sep 17 00:00:00 2001 From: "john.derber" Date: Fri, 13 Aug 2021 14:52:26 +0000 Subject: [PATCH] GitHub Issue NOAA-EMC/GSI#175. Use the global 127L B-Matrix in regional FV3 DA (FV3LAMDA) --- src/gsi/balmod.f90 | 54 ++--- src/gsi/glbsoi.f90 | 3 + src/gsi/gsi_rfv3io_mod.f90 | 23 +- src/gsi/gsimod.F90 | 5 +- src/gsi/m_berror_stats.f90 | 68 +++++- src/gsi/m_berror_stats_reg.f90 | 415 ++++++++++++++++++--------------- src/gsi/prewgt_reg.f90 | 18 +- 7 files changed, 329 insertions(+), 257 deletions(-) diff --git a/src/gsi/balmod.f90 b/src/gsi/balmod.f90 index 59ec5a5b9c..3686efe286 100644 --- a/src/gsi/balmod.f90 +++ b/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) ) diff --git a/src/gsi/glbsoi.f90 b/src/gsi/glbsoi.f90 index d9dbee7b3b..a210ba3258 100644 --- a/src/gsi/glbsoi.f90 +++ b/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) diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90 index fe6a59cb85..7393699ff3 100644 --- a/src/gsi/gsi_rfv3io_mod.f90 +++ b/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 diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 78bebc5140..2081b91f54 100644 --- a/src/gsi/gsimod.F90 +++ b/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 diff --git a/src/gsi/m_berror_stats.f90 b/src/gsi/m_berror_stats.f90 index 82d4cf0d6c..808aee9954 100644 --- a/src/gsi/m_berror_stats.f90 +++ b/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 diff --git a/src/gsi/m_berror_stats_reg.f90 b/src/gsi/m_berror_stats_reg.f90 index 896a43da4a..fe47947078 100644 --- a/src/gsi/m_berror_stats_reg.f90 +++ b/src/gsi/m_berror_stats_reg.f90 @@ -10,10 +10,10 @@ module m_berror_stats_reg use kinds,only : i_kind,r_kind - use constants, only: zero,one,max_varname_length + use constants, only: zero,one,max_varname_length,half use gridmod, only: nsig use chemmod, only : berror_chem,upper2lower,lower2upper - use m_berror_stats, only: berror_stats + use m_berror_stats, only: usenewgfsberror,berror_stats implicit none @@ -24,9 +24,9 @@ module m_berror_stats_reg ! interfaces to file berror_stats. public :: berror_set_reg ! set internal parameters - public :: berror_get_dims_reg ! get dimensions, jfunc::createj_func() - public :: berror_read_bal_reg ! get cross-cov.stats., balmod::prebal() - public :: berror_read_wgt_reg ! get auto-cov.stats., prewgt() + public :: berror_get_dims_reg ! get dimensions, jfunc::createj_func() + public :: berror_read_bal_reg ! get cross-cov.stats., balmod::prebal() + public :: berror_read_wgt_reg ! get auto-cov.stats., prewgt() ! !REVISION HISTORY: ! 25Feb10 - Zhu - adopt code format from m_berror_stats @@ -112,7 +112,6 @@ subroutine berror_set_reg(opt,value) character(len=*),parameter :: myname_=myname//'::berror_set_reg' logical found - found=.false. if(trim(opt)=='cwcoveqqcov') then cwcoveqqcov_=value @@ -153,19 +152,19 @@ subroutine berror_read_bal_reg(msig,mlat,agvi,bvi,wgvi,mype,unit) ! - make changes for generalized control variables !EOP ___________________________________________________________________ - character(len=*),parameter :: myname_=myname//'::berror_read_bal_reg' + character(len=*),parameter :: myname_=myname//'::berror_read_bal_reg' -! workspaces/variables for data not returned +! workspaces/variables for data not returned - integer(i_kind) k,i,m,j,m1,l1,l - integer(i_kind):: nsigstat,nlatstat - integer(i_kind):: inerr + integer(i_kind) k,i,m,j,m1,l1,l + integer(i_kind):: nsigstat,nlatstat + integer(i_kind):: inerr - real(r_kind),dimension(nsig) :: rlsig - real(r_single),dimension(:),allocatable:: clat_avn,sigma_avn - real(r_single),dimension(:,:),allocatable:: bv_avn,wgv_avn - real(r_single),dimension(:,:,:),allocatable:: agv_avn - real(r_kind),dimension(:),allocatable:: rlsigo + real(r_kind),dimension(nsig) :: rlsig + real(r_single),dimension(:),allocatable:: clat_avn,sigma_avn + real(r_single),dimension(:,:),allocatable:: bv_avn,wgv_avn + real(r_single),dimension(:,:,:),allocatable:: agv_avn + real(r_kind),dimension(:),allocatable:: rlsigo ! Open background error statistics file inerr=default_unit_ @@ -186,11 +185,18 @@ subroutine berror_read_bal_reg(msig,mlat,agvi,bvi,wgvi,mype,unit) allocate ( clat_avn(mlat) ) allocate ( sigma_avn(1:msig) ) allocate ( rlsigo(1:msig) ) - allocate ( agv_avn(0:mlat+1,1:msig,1:msig) ) - allocate ( bv_avn(0:mlat+1,1:msig),wgv_avn(0:mlat+1,1:msig) ) + if(usenewgfsberror)then + allocate ( agv_avn(mlat,1:msig,1:msig) ) + allocate ( bv_avn(mlat,1:msig),wgv_avn(mlat,1:msig) ) + else + allocate ( agv_avn(0:mlat+1,1:msig,1:msig) ) + allocate ( bv_avn(0:mlat+1,1:msig),wgv_avn(0:mlat+1,1:msig) ) + end if ! Read background error file to get balance variables read(inerr)clat_avn,(sigma_avn(k),k=1,msig) + + read(inerr)agv_avn,bv_avn,wgv_avn close(inerr) @@ -254,6 +260,8 @@ subroutine berror_read_bal_reg(msig,mlat,agvi,bvi,wgvi,mype,unit) enddo enddo + deallocate (agv_avn,bv_avn,wgv_avn,clat_avn,sigma_avn,rlsigo) + agvi(0,:,:)=agvi(1,:,:) wgvi(0,:)=wgvi(1,:) bvi(0,:)=bvi(1,:) @@ -261,7 +269,6 @@ subroutine berror_read_bal_reg(msig,mlat,agvi,bvi,wgvi,mype,unit) wgvi(mlat+1,:)=wgvi(mlat,:) bvi(mlat+1,:)=bvi(mlat,:) - deallocate (agv_avn,bv_avn,wgv_avn,clat_avn,sigma_avn,rlsigo) return end subroutine berror_read_bal_reg !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -277,7 +284,7 @@ end subroutine berror_read_bal_reg subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qoption,varcw,cwoption,mype,unit) use kinds,only : r_single,r_kind - use gridmod,only : nsig + use gridmod,only : nsig, twodvar_regional use gsi_io, only : verbose use control_vectors,only: nrf,nc2d,nc3d,mvars,nvars use control_vectors,only: cvars => nrf_var @@ -342,9 +349,7 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt real(r_kind),parameter:: vz_oz = 0.53333333_r_kind ! workspace variables not returned - 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:: corqq_avn real(r_single),dimension(:,:),allocatable:: corz_avn,hwll_avn,vztdq_avn real(r_single),dimension(1:mlat,msig,nrf):: corz_tmp @@ -364,10 +369,11 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt integer(i_kind) :: nrf3_dbz integer(i_kind) :: inerr,istat integer(i_kind) :: nsigstat,nlatstat,isig - integer(i_kind) :: loc,m1,m,i,n,j,k,n0,ivar,ic + integer(i_kind) :: loc,m1,m,i,n,j,k,ivar,ic integer(i_kind),allocatable,dimension(:) :: nrf2_loc,nrf3_loc,nmotl_loc real(r_kind) :: factoz real(r_kind) :: raux + real(r_kind),dimension(nsig):: dlsig ! corz = sqrt(corz) real(r_kind), parameter :: corz_default=one,hwll_default=100000_r_kind,vz_default=one @@ -377,13 +383,13 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt ! character(256) :: filename ! filename = 'howv_var_berr.bin' - allocate ( clat_avn(mlat) ) - allocate ( sigma_avn(1:msig) ) - allocate ( agv_avn(0:mlat+1,1:msig,1:msig) ) - allocate ( bv_avn(0:mlat+1,1:msig),wgv_avn(0:mlat+1,1:msig) ) - print_verbose=.false. if(verbose)print_verbose=.true. + + do k=1,nsig + rlsig(k)=log(ges_prslavg(k)/ges_psfcavg) + enddo + ! Open background error statistics file inerr=default_unit_ if(present(unit)) inerr=unit @@ -393,14 +399,6 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt rewind inerr read(inerr) nsigstat,nlatstat -! Read background error file to get balance variables - read(inerr)clat_avn,(sigma_avn(k),k=1,msig) - read(inerr)agv_avn,bv_avn,wgv_avn - -! compute vertical(pressure) interpolation index and weight - do k=1,nsig - rlsig(k)=log(ges_prslavg(k)/ges_psfcavg) - enddo if(mype==0) then write(6,*) myname_,'(PREWGT_REG): read error amplitudes ', & @@ -409,16 +407,20 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt mype,nsigstat,nlatstat end if - allocate(nrf3_loc(nc3d),nrf2_loc(nc2d),nmotl_loc(mvars)) - do n=1,nc3d - nrf3_loc(n)=getindex(cvars,cvars3d(n)) - enddo - do n=1,nc2d - nrf2_loc(n)=getindex(cvars,cvars2d(n)) - enddo - do n=1,mvars - nmotl_loc(n)=getindex(cvars,cvarsmd(n)) - enddo +! Read background error file to get past alance variables + read(inerr) + read(inerr) + + allocate(nrf3_loc(nc3d),nrf2_loc(nc2d),nmotl_loc(mvars)) + do n=1,nc3d + nrf3_loc(n)=getindex(cvars,cvars3d(n)) + enddo + do n=1,nc2d + nrf2_loc(n)=getindex(cvars,cvars2d(n)) + enddo + do n=1,mvars + nmotl_loc(n)=getindex(cvars,cvarsmd(n)) + enddo ! Read amplitudes @@ -427,44 +429,56 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt if (berror_chem) then read(inerr,iostat=istat) varshort,isig var=upper2lower(varshort) - if (var == 'pm25') var = 'pm2_5' + if (trim(var) == 'pm25') var = 'pm2_5' else read(inerr,iostat=istat) varshort, isig var=varshort endif + if (istat /= 0) exit read + do n=1,nrf + if (trim(var)==cvars(n)) then + nrf_err(n)=.true. + loc=n + exit + else + loc=-999 + end if + end do + if(loc == -999)then + if(mype == 0)write(6,*) 'variable in input file, but not used in analysis ',var,isig + read(inerr) + read(inerr) + if(isig > 1)read(inerr) + cycle read + end if - if (istat /= 0) exit allocate ( corz_avn(1:mlat,1:isig) ) - allocate ( hwll_avn(0:mlat+1,1:isig) ) - allocate ( vztdq_avn(1:isig,0:mlat+1) ) - if (var/='q' .or. (var=='cw' .and. cwoption==2)) then + if (trim(var)/='q' .or. (trim(var)=='cw' .and. cwoption==2)) then read(inerr) corz_avn else allocate ( corqq_avn(1:mlat,1:isig) ) read(inerr) corz_avn,corqq_avn end if - read(inerr) hwll_avn - if (isig>1) then - read(inerr) vztdq_avn + if(usenewgfsberror)then + allocate ( hwll_avn(mlat,1:isig) ) + else + allocate ( hwll_avn(0:mlat+1,1:isig) ) end if + read(inerr) hwll_avn -! load the variances - do n=1,nrf - if (var==cvars(n)) then - nrf_err(n)=.true. - loc=n - exit - else - loc=-999 - end if - end do if (isig==msig) then + if(usenewgfsberror)then + allocate ( vztdq_avn(1:isig,mlat) ) + else + allocate ( vztdq_avn(1:isig,0:mlat+1) ) + end if + read(inerr) vztdq_avn do n=1,nc3d if (nrf3_loc(n)==loc) then - if ((var=='q' .and. qoption==2) .or. (var=='cw' .and. cwoption==2)) then + if ((trim(var)=='q' .and. qoption==2) .or. (trim(var)=='cw' .and. cwoption==2)) then ! choose which q stat to use do k=1,msig do i=1,mlat @@ -479,43 +493,55 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt end do end if do k=1,msig - do i=0,mlat+1 - hwll_tmp(i,k,n)=hwll_avn(i,k) - vz_tmp(k,i,n)=vztdq_avn(k,i) - end do + if(usenewgfsberror)then + do i=1,mlat + hwll_tmp(i,k,n)=hwll_avn(i,k) + vz_tmp(k,i,n)=vztdq_avn(k,i) + end do + hwll_tmp(0,k,n)=hwll_avn(1,k) + hwll_tmp(mlat+1,k,n)=hwll_avn(mlat,k) + vz_tmp(k,0,n)=vztdq_avn(k,1) + vz_tmp(k,mlat+1,n)=vztdq_avn(k,mlat) + else + do i=0,mlat+1 + hwll_tmp(i,k,n)=hwll_avn(i,k) + vz_tmp(k,i,n)=vztdq_avn(k,i) + end do + end if end do exit end if end do - end if + deallocate ( vztdq_avn ) - if (isig==1) then + else if(isig == 1)then do n=1,nc2d if (nrf2_loc(n)==loc) then do i=1,mlat corp(i,n)=corz_avn(i,1) hwllp(i,n)=hwll_avn(i,1) end do - hwllp(0,n)=hwll_avn(0,1) - hwllp(mlat+1,n)=hwll_avn(mlat+1,1) + if(usenewgfsberror)then + hwllp(0,n)=hwll_avn(1,1) + hwllp(mlat+1,n)=hwll_avn(mlat,1) + else + hwllp(0,n)=hwll_avn(0,1) + hwllp(mlat+1,n)=hwll_avn(mlat+1,1) + end if exit end if end do end if - deallocate ( corz_avn ) deallocate ( hwll_avn ) - deallocate ( vztdq_avn ) - if (var=='q' .or. var=='cw') deallocate ( corqq_avn ) + if(allocated(corqq_avn)) deallocate ( corqq_avn ) enddo read close(inerr) - deallocate(clat_avn,sigma_avn) - deallocate(agv_avn,bv_avn,wgv_avn) - ! 3d variable do n=1,nc3d loc=nrf3_loc(n) + if (loc <= 0)cycle if (nrf_err(loc)) then do k=1,nsig m=lsig(k) @@ -599,39 +625,38 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt vz(:,:,nrf3_oz)=vz_oz end if - if (cwcoveqqcov_ .and. nrf3_cw>0) then - corz(:,:,nrf3_cw)=corz(:,:,nrf3_q) - hwll(:,:,nrf3_cw)=hwll(:,:,nrf3_q) - vz(:,:,nrf3_cw)=vz(:,:,nrf3_q) - end if + if(nrf3_cw > 0)then + if (cwcoveqqcov_ ) then + corz(:,:,nrf3_cw)=corz(:,:,nrf3_q) + hwll(:,:,nrf3_cw)=hwll(:,:,nrf3_q) + vz(:,:,nrf3_cw)=vz(:,:,nrf3_q) - if ((.not. cwcoveqqcov_) .and. nrf3_cw>0) then - corz(:,:,nrf3_cw)=zero - if (cwoption==2) then - do k=1,nsig - if (ges_prslavg(k)>15.0_r_kind) then - do j=1,mlat - varcw(j,k)=max(real(corz(j,k,nrf3_cw),r_kind),zero) - corz(j,k,nrf3_cw)=one - enddo - end if - enddo - end if + else + corz(:,:,nrf3_cw)=zero + if (cwoption==2) then + do k=1,nsig + if (ges_prslavg(k)>15.0_r_kind) then + do j=1,mlat + varcw(j,k)=max(real(corz(j,k,nrf3_cw),r_kind),zero) + corz(j,k,nrf3_cw)=one + enddo + end if + enddo + end if - if (cwoption==1 .or. cwoption==3) then - do k=1,nsig - if (ges_prslavg(k)>15.0_r_kind) then - do j=1,mlat - corz(j,k,nrf3_cw)=one - end do - end if - end do - hwll(:,:,nrf3_cw)=0.5_r_kind*hwll(:,:,nrf3_q) - vz(:,:,nrf3_cw)=0.5_r_kind*vz(:,:,nrf3_q) + if (cwoption==1 .or. cwoption==3) then + do k=1,nsig + if (ges_prslavg(k)>15.0_r_kind) then + do j=1,mlat + corz(j,k,nrf3_cw)=one + end do + end if + end do + hwll(:,:,nrf3_cw)=0.5_r_kind*hwll(:,:,nrf3_q) + vz(:,:,nrf3_cw)=0.5_r_kind*vz(:,:,nrf3_q) + end if end if - end if - - if (icloud_cv .and. n_clouds_fwd>0 .and. nrf3_cw<=0 .and. cwoption==3) then + else if (icloud_cv .and. n_clouds_fwd>0 .and. cwoption==3) then do n=1,size(cvars3d) do ic=1,n_clouds_fwd if(trim(cvars3d(n))==trim(cloud_names_fwd(ic))) then @@ -665,50 +690,46 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt ! 2d variable do n=1,nc2d loc=nrf2_loc(n) - if (nrf_err(loc)) cycle + if(loc <= 0)cycle +! If we want to use the sst in global file uncomment the following if statement if (n==nrf2_sst) then - do i=1,mlat - corp(i,n)=zero_3 - end do - do i=0,mlat+1 - hwllp(i,n)=hwll(i,1,nrf3_sf) - hwllp(i,nc2d+1)=hwll(i,1,nrf3_sf) !not very nice, since it assumes that stl and sti - hwllp(i,nc2d+2)=hwll(i,1,nrf3_sf) !are always the first motley variables in convinfo - end do - end if - if (n==nrf2_gust) then +! if(.not. usenewgfsberror)then + do i=1,mlat + corp(i,n)=zero_3 + end do + do i=0,mlat+1 + hwllp(i,n)=hwll(i,1,nrf3_sf) + end do +! end if + else if (n==nrf2_gust) then do i=1,mlat corp(i,n)=three end do do i=0,mlat+1 hwllp(i,n)=hwll(i,1,nrf3_q) end do - end if - if (n==nrf2_vis) then + else if (n==nrf2_vis) then do i=1,mlat corp(i,n)=3.0_r_kind end do do i=0,mlat+1 hwllp(i,n)=hwll(i,1,nrf3_t) end do - end if - if (n==nrf2_pblh) then + else if (n==nrf2_pblh) then do i=1,mlat corp(i,n)=500.0_r_kind end do do i=0,mlat+1 hwllp(i,n)=three*hwll(i,1,nrf3_t) end do - end if - if (n==nrf2_wspd10m) then + else if (n==nrf2_wspd10m) then do i=1,mlat corp(i,n)=three end do do i=0,mlat+1 hwllp(i,n)=hwll(i,1,nrf3_q) end do - end if - if (n==nrf2_td2m) then + else if (n==nrf2_td2m) then raux=maxval(corz(1:mlat,1,nrf3_q)) do i=1,mlat corp(i,n)=(corz(i,1,nrf3_q)/raux)*three !tentatively @@ -716,32 +737,28 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt do i=0,mlat+1 hwllp(i,n)=hwll(i,1,nrf3_q) !tentatively end do - end if - if (n==nrf2_mxtm) then + else if (n==nrf2_mxtm) then do i=1,mlat corp(i,n)=corz(i,1,nrf3_t) end do do i=0,mlat+1 hwllp(i,n)=hwll(i,1,nrf3_t) end do - end if - if (n==nrf2_mitm) then + else if (n==nrf2_mitm) then do i=1,mlat corp(i,n)=corz(i,1,nrf3_t) end do do i=0,mlat+1 hwllp(i,n)=hwll(i,1,nrf3_t) end do - end if - if (n==nrf2_pmsl) then + else if (n==nrf2_pmsl) then do i=1,mlat corp(i,n)=corp(i,nrf2_ps) end do do i=0,mlat+1 hwllp(i,n)=hwllp(i,nrf2_ps) end do - end if - if (n==nrf2_howv) then + else if (n==nrf2_howv) then call read_howv_stats(mlat,1,2,cov_dum) do i=1,mlat corp(i,n)=cov_dum(i,1,1) !#ww3 @@ -759,16 +776,14 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt ! do i=0,mlat+1 ! hwllp(i,n)=hwll(i,1,nrf3_sf) !tentatively !#ww3 hwllp(i,n)=150000_r_kind ! ! end do - end if - if (n==nrf2_tcamt) then + else if (n==nrf2_tcamt) then do i=1,mlat corp(i,n)=50.0_r_kind end do do i=0,mlat+1 hwllp(i,n)=hwll(i,1,nrf3_t) end do - end if - if (n==nrf2_lcbas) then + else if (n==nrf2_lcbas) then do i=1,mlat corp(i,n)=40000.0_r_kind end do @@ -776,8 +791,7 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt hwllp(i,n)=hwll(i,1,nrf3_t) end do if(print_verbose)print*, 'm_berror_reg: maxhwllp_lcbas=',maxval(hwllp(:,n)) - end if - if (n==nrf2_cldch) then + else if (n==nrf2_cldch) then do i=1,mlat corp(i,n)=3.0_r_kind end do @@ -785,16 +799,14 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt hwllp(i,n)=hwll(i,1,nrf3_t) end do if(print_verbose)print*, 'm_berror_reg: maxhwllp_cldch=',maxval(hwllp(:,n)) - end if - if (n==nrf2_uwnd10m) then + else if (n==nrf2_uwnd10m) then do i=1,mlat corp(i,n)=three end do do i=0,mlat+1 hwllp(i,n)=hwll(i,1,nrf3_q) end do - end if - if (n==nrf2_vwnd10m) then + else if (n==nrf2_vwnd10m) then do i=1,mlat corp(i,n)=three end do @@ -807,103 +819,126 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt ! motley variable - n0=nc2d do n=1,mvars - if (cvarsmd(n)=='stl') cycle - if (cvarsmd(n)=='sti') cycle + loc=nmotl_loc(n) + if(loc <=0)cycle + loc=n+nc2d + if (cvarsmd(n)=='sti' .or. cvarsmd(n)=='stl') then + do i=1,mlat + corp(i,loc)=corz(i,1,nrf3_sf) + end do + do i=0,mlat+1 + hwllp(i,loc)=hwll(i,1,nrf3_sf) + end do - if (cvarsmd(n)=='pswter') then + else if (cvarsmd(n)=='pswter') then do i=1,mlat - corp(i,n0+n)=corp(i,nrf2_ps) + corp(i,loc)=corp(i,nrf2_ps) end do do i=0,mlat+1 - hwllp(i,n0+n)=hwllp(i,nrf2_ps) + hwllp(i,loc)=hwllp(i,nrf2_ps) end do - endif - if (cvarsmd(n)=='twter') then + else if (cvarsmd(n)=='twter') then do i=1,mlat - corp(i,n0+n)=corz(i,1,nrf3_t) + corp(i,loc)=corz(i,1,nrf3_t) end do do i=0,mlat+1 - hwllp(i,n0+n)=hwll(i,1,nrf3_t) + hwllp(i,loc)=hwll(i,1,nrf3_t) end do - endif - if (cvarsmd(n)=='qwter') then + else if (cvarsmd(n)=='qwter') then do i=1,mlat - corp(i,n0+n)=corz(i,1,nrf3_q) + corp(i,loc)=corz(i,1,nrf3_q) end do do i=0,mlat+1 - hwllp(i,n0+n)=hwll(i,1,nrf3_q) + hwllp(i,loc)=hwll(i,1,nrf3_q) end do - endif - if (cvarsmd(n)=='gustwter') then + else if (cvarsmd(n)=='gustwter') then do i=1,mlat - corp(i,n0+n)=corp(i,nrf2_gust) + corp(i,loc)=corp(i,nrf2_gust) end do do i=0,mlat+1 - hwllp(i,n0+n)=hwllp(i,nrf2_gust) + hwllp(i,loc)=hwllp(i,nrf2_gust) end do - endif - if (cvarsmd(n)=='wspd10mwter') then + else if (cvarsmd(n)=='wspd10mwter') then do i=1,mlat - corp(i,n0+n)=corp(i,nrf2_wspd10m) + corp(i,loc)=corp(i,nrf2_wspd10m) end do do i=0,mlat+1 - hwllp(i,n0+n)=hwllp(i,nrf2_wspd10m) + hwllp(i,loc)=hwllp(i,nrf2_wspd10m) end do - endif - if (cvarsmd(n)=='td2mwter') then + else if (cvarsmd(n)=='td2mwter') then do i=1,mlat - corp(i,n0+n)=corp(i,nrf2_td2m) + corp(i,loc)=corp(i,nrf2_td2m) end do do i=0,mlat+1 - hwllp(i,n0+n)=hwllp(i,nrf2_td2m) + hwllp(i,loc)=hwllp(i,nrf2_td2m) end do - endif - if (cvarsmd(n)=='mxtmwter') then + else if (cvarsmd(n)=='mxtmwter') then do i=1,mlat - corp(i,n0+n)=corp(i,nrf2_mxtm) + corp(i,loc)=corp(i,nrf2_mxtm) end do do i=0,mlat+1 - hwllp(i,n0+n)=hwllp(i,nrf2_mxtm) + hwllp(i,loc)=hwllp(i,nrf2_mxtm) end do - endif - if (cvarsmd(n)=='mitmwter') then + else if (cvarsmd(n)=='mitmwter') then do i=1,mlat - corp(i,n0+n)=corp(i,nrf2_mitm) + corp(i,loc)=corp(i,nrf2_mitm) end do do i=0,mlat+1 - hwllp(i,n0+n)=hwllp(i,nrf2_mitm) + hwllp(i,loc)=hwllp(i,nrf2_mitm) end do - endif - if (cvarsmd(n)=='uwnd10mwter') then + else if (cvarsmd(n)=='uwnd10mwter') then do i=1,mlat - corp(i,n0+n)=corp(i,nrf2_uwnd10m) + corp(i,loc)=corp(i,nrf2_uwnd10m) end do do i=0,mlat+1 - hwllp(i,n0+n)=hwllp(i,nrf2_uwnd10m) + hwllp(i,loc)=hwllp(i,nrf2_uwnd10m) end do - endif - if (cvarsmd(n)=='vwnd10mwter') then + else if (cvarsmd(n)=='vwnd10mwter') then do i=1,mlat - corp(i,n0+n)=corp(i,nrf2_vwnd10m) + corp(i,loc)=corp(i,nrf2_vwnd10m) end do do i=0,mlat+1 - hwllp(i,n0+n)=hwllp(i,nrf2_vwnd10m) + hwllp(i,loc)=hwllp(i,nrf2_vwnd10m) end do - endif +! if not found use default + else + do i=1,mlat + corp(i,loc)=corz_default + end do + do i=0,mlat+1 + hwllp(i,loc)=hwll_default + end do + end if enddo deallocate(nrf3_loc,nrf2_loc,nmotl_loc) +! Normalize vz with del sigmma and convert to vertical grid units! + if(.not. twodvar_regional)then + dlsig(1)=rlsig(1)-rlsig(2) + do k=2,nsig-1 + dlsig(k)=half*(rlsig(k-1)-rlsig(k+1)) + enddo + dlsig(nsig)=rlsig(nsig-1)-rlsig(nsig) + + do n=1,nc3d + do j=0,mlat+1 + do k=1,nsig + vz(k,j,n)=vz(k,j,n)*dlsig(k) + end do + end do + end do + end if + return end subroutine berror_read_wgt_reg diff --git a/src/gsi/prewgt_reg.f90 b/src/gsi/prewgt_reg.f90 index 5682a38896..0e2804dea1 100644 --- a/src/gsi/prewgt_reg.f90 +++ b/src/gsi/prewgt_reg.f90 @@ -117,7 +117,7 @@ subroutine prewgt_reg(mype) real(r_kind) samp2,dl1,dl2,d real(r_kind) samp,hwl,cc - real(r_kind),dimension(nsig):: rate,dlsig,rlsig + real(r_kind),dimension(nsig):: rate,rlsig real(r_kind),dimension(nsig,nsig):: turn real(r_kind),dimension(ny,nx)::sl real(r_kind) fact,psfc015 @@ -161,7 +161,7 @@ subroutine prewgt_reg(mype) ! Allocate arrays in stats file allocate ( corz(1:mlat,1:nsig,1:nc3d) ) - allocate ( corp(1:mlat,nc2d) ) + allocate ( corp(1:mlat,nvars-nc3d) ) allocate ( hwll(0:mlat+1,1:nsig,1:nc3d),hwllp(0:mlat+1,nvars-nc3d) ) allocate ( vz(1:nsig,0:mlat+1,1:nc3d) ) @@ -223,20 +223,6 @@ subroutine prewgt_reg(mype) enddo enddo endif ! regional_ozone -! Normalize vz with del sigmma and convert to vertical grid units! - dlsig(1)=rlsig(1)-rlsig(2) - do k=2,nsig-1 - dlsig(k)=half*(rlsig(k-1)-rlsig(k+1)) - enddo - dlsig(nsig)=rlsig(nsig-1)-rlsig(nsig) - - do n=1,nc3d - do j=0,mlat+1 - do k=1,nsig - vz(k,j,n)=vz(k,j,n)*dlsig(k) - end do - end do - end do ! As used in the code, the horizontal length scale ! parameters are used in an inverted form. Invert