From b846b29ef16136828b1e08504f45bb1f05398706 Mon Sep 17 00:00:00 2001 From: "john.derber" Date: Fri, 13 Aug 2021 14:52:26 +0000 Subject: [PATCH 1/5] 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 From f938842c186e615fcf711f286b04ba0abb64c73b Mon Sep 17 00:00:00 2001 From: "john.derber" Date: Fri, 22 Oct 2021 15:09:39 +0000 Subject: [PATCH 2/5] GitHub Issue NOAA-EMC/GSI#219 Improve Minimization and fix bug in vqc --- regression/regression_var.sh | 2 +- src/gsi/berror.f90 | 143 ++++++--------------------- src/gsi/compute_derived.f90 | 5 +- src/gsi/compute_qvar3d.f90 | 4 +- src/gsi/constants.f90 | 1 - src/gsi/crtm_interface.f90 | 31 ++---- src/gsi/evalqlim.f90 | 8 +- src/gsi/genqsat.f90 | 3 +- src/gsi/gsimod.F90 | 8 +- src/gsi/intjcmod.f90 | 6 +- src/gsi/jfunc.f90 | 8 +- src/gsi/observer.F90 | 8 +- src/gsi/pcgsoi.f90 | 110 +++++++++------------ src/gsi/prt_guess.f90 | 8 +- src/gsi/read_guess.F90 | 4 +- src/gsi/setupcldtot.F90 | 1 - src/gsi/setupq.f90 | 15 ++- src/gsi/setuprhsall.f90 | 5 - src/gsi/setupw.f90 | 29 +++--- src/gsi/stpcalc.f90 | 185 ++++++++++++++++++----------------- src/gsi/stpjcmod.f90 | 14 +-- src/gsi/stpw.f90 | 38 ++++--- src/gsi/update_guess.f90 | 4 +- 23 files changed, 268 insertions(+), 372 deletions(-) diff --git a/regression/regression_var.sh b/regression/regression_var.sh index 37096edde1..92af951596 100755 --- a/regression/regression_var.sh +++ b/regression/regression_var.sh @@ -131,7 +131,7 @@ case $machine in # On Hera, there are no scrubbers to remove old contents from stmp* directories. # After completion of regression tests, will remove the regression test subdirecories - export clean=".true." + export clean=".false." ;; Jet) diff --git a/src/gsi/berror.f90 b/src/gsi/berror.f90 index 71bd4f1ea8..cc8ae09208 100644 --- a/src/gsi/berror.f90 +++ b/src/gsi/berror.f90 @@ -33,7 +33,6 @@ module berror ! 2011-04-07 todling - move newpc4pred to radinfo ! 2012-10-09 Gu - add fut2ps to project unbalanced temp to surface pressure in static B modeling ! 2013-05-27 zhu - add background error variances for aircraft temperature bias correction coefficients -! 2013-10-02 zhu - add reset_predictors_var ! ! subroutines included: ! sub init_berror - initialize background error related variables @@ -120,7 +119,6 @@ module berror public :: create_berror_vars public :: destroy_berror_vars public :: set_predictors_var - public :: reset_predictors_var public :: init_rftable public :: initable public :: create_berror_vars_reg @@ -385,7 +383,6 @@ subroutine set_predictors_var integer(i_kind) i,j,ii real(r_kind) stndev - real(r_kind) obs_count logical new_tail stndev = one/biaspredvar @@ -407,17 +404,12 @@ subroutine set_predictors_var do j=1,npred ii=ii+1 if (inew_rad(i)) then - varprd(ii)=10000.0_r_kind + varA(j,i)=r10 else - if (ostats(i)<=20.0_r_kind) then - varA(j,i)=two*varA(j,i)+1.0e-6_r_kind - varprd(ii)=varA(j,i) - else - varprd(ii)=1.1_r_kind*varA(j,i)+1.0e-6_r_kind - end if - if (varprd(ii)>r10) varprd(ii)=r10 - if (varA(j,i)>10000.0_r_kind) varA(j,i)=10000.0_r_kind + varA(j,i)=1.1_r_kind*varA(j,i)+1.0e-6_r_kind + varA(j,i)= min(r10,varA(j,i)) end if + varprd(ii)=varA(j,i) end do end do @@ -427,50 +419,33 @@ subroutine set_predictors_var do j=1,npredt ii=ii+1 - if (aircraft_t_bc_pof) then - obs_count = ostats_t(j,i) - new_tail = varA_t(j,i)==zero - end if + new_tail = varA_t(j,i)==zero if (aircraft_t_bc) then - obs_count = ostats_t(1,i) new_tail = .true. if (any(varA_t(:,i)/=zero)) new_tail = .false. end if if (new_tail) then - varprd(ii)=one_tenth*one_tenth - if (aircraft_t_bc .and. j==2) varprd(ii)=1.0e-4_r_kind - if (aircraft_t_bc .and. j==3) varprd(ii)=1.0e-5_r_kind - else - if (obs_count<=10.0_r_kind) then - if (aircraft_t_bc .and. j==2) then - varA_t(j,i)=1.01_r_kind*varA_t(j,i)+1.0e-6_r_kind - else if (aircraft_t_bc .and. j==3) then - varA_t(j,i)=1.01_r_kind*varA_t(j,i)+1.0e-7_r_kind - else - varA_t(j,i)=1.01_r_kind*varA_t(j,i)+1.0e-5_r_kind - end if - varprd(ii)=varA_t(j,i) + if (aircraft_t_bc .and. j==2) then + varA_t(j,i)=1.0e-4_r_kind + else if (aircraft_t_bc .and. j==3) then + varA_t(j,i)=1.0e-5_r_kind else - if (aircraft_t_bc .and. j==2) then - varprd(ii)=1.005_r_kind*varA_t(j,i)+1.0e-6_r_kind - else if (aircraft_t_bc .and. j==3) then - varprd(ii)=1.005_r_kind*varA_t(j,i)+1.0e-7_r_kind - else - varprd(ii)=1.005_r_kind*varA_t(j,i)+1.0e-5_r_kind - end if + varA_t(j,i)=one_tenth*one_tenth end if - if (varprd(ii)>one_tenth) varprd(ii)=one_tenth - if (varA_t(j,i)>one_tenth) varA_t(j,i)=one_tenth + else if (aircraft_t_bc .and. j==2) then - if (varprd(ii)>1.0e-3_r_kind) varprd(ii)=1.0e-3_r_kind - if (varA_t(j,i)>1.0e-3_r_kind) varA_t(j,i)=1.0e-3_r_kind - end if - if (aircraft_t_bc .and. j==3) then - if (varprd(ii)>1.0e-4_r_kind) varprd(ii)=1.0e-4_r_kind - if (varA_t(j,i)>1.0e-4_r_kind) varA_t(j,i)=1.0e-4_r_kind + varA_t(j,i)=1.005_r_kind*varA_t(j,i)+1.0e-6_r_kind + varA_t(j,i)=min(varA_t(j,i),1.0e-3_r_kind) + else if (aircraft_t_bc .and. j==3) then + varA_t(j,i)=1.005_r_kind*varA_t(j,i)+1.0e-7_r_kind + varA_t(j,i)=min(varA_t(j,i),1.0e-4_r_kind) + else + varA_t(j,i)=1.005_r_kind*varA_t(j,i)+1.0e-5_r_kind + varA_t(j,i)=min(varA_t(j,i),one_tenth) end if end if + varprd(ii)=varA_t(j,i) end do end do end if @@ -479,70 +454,6 @@ subroutine set_predictors_var return end subroutine set_predictors_var - - subroutine reset_predictors_var -!$$$ subprogram documentation block -! . . . . -! subprogram: reset_predictors_var sets variances for bias correction predictors -! prgmmr: yanqiu org: np20 date: 2013-10-01 -! -! abstract: resets variances for bias correction predictors -! -! program history log: -! output argument list: -! 2013-10-01 zhu -! -! attributes: -! language: f90 -! machine: ibm RS/6000 SP -! -!$$$ - use constants, only: one,one_tenth - use radinfo, only: newpc4pred,jpch_rad,npred,ostats,inew_rad,iuse_rad - use aircraftinfo, only: aircraft_t_bc_pof,aircraft_t_bc,biaspredt,ntail,npredt,ostats_t - use gridmod, only: twodvar_regional - use jfunc, only: nrclen, ntclen - implicit none - - integer(i_kind) i,j,ii,obs_count - real(r_kind) stndev - - stndev = one/biaspredt - -! reset variances for bias predictor coeff. based on current data count - if (.not. twodvar_regional .and. newpc4pred) then - ii=0 - do i=1,jpch_rad - do j=1,npred - ii=ii+1 - if (.not.inew_rad(i) .and. iuse_rad(i)>0 .and. ostats(i)<=20.0_r_kind) then - varprd(ii)=1.0e-6_r_kind - end if - end do - end do - - if ((aircraft_t_bc_pof .or. aircraft_t_bc) .and. ntclen>0) then - ii=nrclen-ntclen - do i=1,ntail - do j=1,npredt - ii=ii+1 - - if (aircraft_t_bc_pof) obs_count = ostats_t(j,i) - if (aircraft_t_bc) obs_count = ostats_t(1,i) - - if (obs_count<=10.0_r_kind .and. varprd(ii)>stndev) then - varprd(ii)=stndev - if (aircraft_t_bc .and. j==2) varprd(ii)=one_tenth*stndev - if (aircraft_t_bc .and. j==3) varprd(ii)=one_tenth*one_tenth*stndev - end if - end do - end do - end if - end if - - return - end subroutine reset_predictors_var - subroutine pcinfo !$$$ subprogram documentation block ! . . . . @@ -590,13 +501,12 @@ subroutine pcinfo do i=1,jpch_rad do j=1,npred ii=ii+1 -! if (ostats(i)>zero) vprecond(nclen1+ii)=vprecond(nclen1+ii)/(one+rstats(j,i)*varprd(ii)) if (ostats(i)>zero) vprecond(nclen1+ii)=one/(one+rstats(j,i)*varprd(ii)) if (ostats(i)>20.0_r_kind) then if (rstats(j,i)>zero) then varA(j,i)=one/(one/varprd(ii)+rstats(j,i)) else - varA(j,i)=10000.0_r_kind + if(varA(j,i) <= zero)varA(j,i)=10000.0_r_kind end if end if end do @@ -612,13 +522,18 @@ subroutine pcinfo ii=ii+1 jj=jj+1 - if (aircraft_t_bc_pof) obs_count = ostats_t(j,i) - if (aircraft_t_bc) obs_count = ostats_t(1,i) + obs_count=0 + if (aircraft_t_bc_pof) then + obs_count = ostats_t(j,i) + else if (aircraft_t_bc) then + obs_count = ostats_t(1,i) + end if -! if (obs_count>zero) vprecond(nclen1+ii)=vprecond(nclen1+ii)/(one+rstats_t(j,i)*varprd(jj)) if (obs_count>zero) vprecond(nclen1+ii)=one/(one+rstats_t(j,i)*varprd(jj)) if (obs_count>3.0_r_kind) then varA_t(j,i)=one/(one/varprd(jj)+rstats_t(j,i)) + else + if(varA_t(j,i) <= zero)varA_t(j,i)=10000.0_r_kind end if end do end do diff --git a/src/gsi/compute_derived.f90 b/src/gsi/compute_derived.f90 index a33dbe1f41..deb9200da8 100644 --- a/src/gsi/compute_derived.f90 +++ b/src/gsi/compute_derived.f90 @@ -89,7 +89,7 @@ subroutine compute_derived(mype,init_pass) use kinds, only: r_kind,i_kind use jfunc, only: jiter,jiterstart,& qoption,switch_on_derivatives,& - tendsflag,clip_supersaturation + tendsflag,superfact,clip_supersaturation use control_vectors, only: cvars3d use control_vectors, only: nrf_var use control_vectors, only: an_amp0 @@ -178,7 +178,6 @@ subroutine compute_derived(mype,init_pass) if(init_pass .and. (ntguessig<1 .or. ntguessig>nfldsig)) & call die(myname,'invalid init_pass, ntguessig =',ntguessig) - ! Get required indexes from control vector names nrf3_q=getindex(cvars3d,'q') iq_loc=getindex(nrf_var,'q') @@ -207,7 +206,7 @@ subroutine compute_derived(mype,init_pass) end do end do end do - + ! Load guess cw for use in inner loop ! Get pointer to cloud water mixing ratio it=ntguessig diff --git a/src/gsi/compute_qvar3d.f90 b/src/gsi/compute_qvar3d.f90 index 54aa7721bf..851212d4f5 100644 --- a/src/gsi/compute_qvar3d.f90 +++ b/src/gsi/compute_qvar3d.f90 @@ -37,7 +37,7 @@ subroutine compute_qvar3d !$$$ use kinds, only: r_kind,i_kind,r_single use berror, only: dssv - use jfunc, only: varq,qoption,varcw,cwoption,clip_supersaturation + use jfunc, only: varq,qoption,varcw,cwoption,clip_supersaturation,superfact use derivsmod, only: qsatg,qgues use control_vectors, only: cvars3d use gridmod, only: lat2,lon2,nsig @@ -94,7 +94,7 @@ subroutine compute_qvar3d ! Limit q to be >= qmin ges_q(i,j,k)=max(ges_q(i,j,k),qmin) ! Limit q to be <= ges_qsat - if(clip_supersaturation) ges_q(i,j,k)=min(ges_q(i,j,k),ges_qsat(i,j,k,it)) + if(clip_supersaturation) ges_q(i,j,k)=min(ges_q(i,j,k),superfact*ges_qsat(i,j,k,it)) end do end do end do diff --git a/src/gsi/constants.f90 b/src/gsi/constants.f90 index da2368d70f..77aa027061 100644 --- a/src/gsi/constants.f90 +++ b/src/gsi/constants.f90 @@ -238,7 +238,6 @@ module constants real(r_kind),parameter:: ke2 = 0.00002_r_kind real(r_kind),parameter:: row = r1000 real(r_kind),parameter:: rrow = one/row -! real(r_kind),parameter:: qmin = 1.e-7_r_kind !lower bound on ges_q ! Constant used to process ozone real(r_kind),parameter:: constoz = 604229.0_r_kind diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index d03715059f..1b90b5b279 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -171,8 +171,6 @@ module crtm_interface real(r_kind) , save ,allocatable,dimension(:,:) :: cloud_efr ! effective radius of cloud type in CRTM real(r_kind) , save ,allocatable,dimension(:) :: cf ! effective radius of cloud type in CRTM real(r_kind) , save ,allocatable,dimension(:) :: hwp_guess ! column total for each hydrometeor - - real(r_kind) , save ,allocatable,dimension(:,:,:,:) :: gesqsat ! qsat to calc rh for aero particle size estimate real(r_kind) , save ,allocatable,dimension(:) :: table,table2,tablew ! GFDL saturation water vapor pressure tables real(r_kind) , save ,allocatable,dimension(:) :: des2,desw ! GFDL saturation water vapor presure real(r_kind) , save ,allocatable,dimension(:) :: lcloud4crtm_wk ! cloud info usage index for each channel @@ -840,16 +838,6 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo endif ! nvege_type endif ! regional or IGBP -! Calculate RH when aerosols are present and/or cloud-fraction used - if (n_actual_aerosols_wk>0 .or. n_clouds_fwd_wk>0) then - allocate(gesqsat(lat2,lon2,nsig,nfldsig)) - ice=.true. - iderivative=0 - do ii=1,nfldsig - call genqsat(gesqsat(1,1,1,ii),ges_tsen(1,1,1,ii),ges_prsl(1,1,1,ii),lat2,lon2,nsig,ice,iderivative) - end do - endif - ! Initial GFDL saturation water vapor pressure tables if (n_actual_aerosols_wk>0 .or. n_clouds_fwd_wk>0 .and. imp_physics==11) then @@ -905,7 +893,6 @@ subroutine destroy_crtm if (error_status /= success) & write(6,*)myname_,': ***ERROR*** error_status=',error_status if (n_actual_aerosols_wk>0 .or. n_clouds_fwd_wk>0) then - deallocate(gesqsat) if (imp_physics==11) then deallocate(table) deallocate(table2) @@ -1048,7 +1035,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & use radinfo, only: nsigradjac use gsi_nstcouplermod, only: nst_gsi use guess_grids, only: ges_tsen,& - ges_prsl,ges_prsi,tropprs,dsfct,add_rtm_layers, & + ges_prsl,ges_prsi,ges_qsat,tropprs,dsfct,add_rtm_layers, & hrdifsig,nfldsig,hrdifsfc,nfldsfc,ntguessfc,isli2,sno2, & hrdifaer,nfldaer ! for separate aerosol input file use cloud_efr_mod, only: efr_ql,efr_qi,efr_qr,efr_qs,efr_qg,efr_qh @@ -1908,14 +1895,14 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & end if ! lread_ext_aerosol end if ! n_actual_aerosols_wk > 0 do k=1,nsig - qs(k) = (gesqsat(ix ,iy ,k,itsig )*w00+ & - gesqsat(ixp,iy ,k,itsig )*w10+ & - gesqsat(ix ,iyp,k,itsig )*w01+ & - gesqsat(ixp,iyp,k,itsig )*w11)*dtsig + & - (gesqsat(ix ,iy ,k,itsigp)*w00+ & - gesqsat(ixp,iy ,k,itsigp)*w10+ & - gesqsat(ix ,iyp,k,itsigp)*w01+ & - gesqsat(ixp,iyp,k,itsigp)*w11)*dtsigp + qs(k) = (ges_qsat(ix ,iy ,k,itsig )*w00+ & + ges_qsat(ixp,iy ,k,itsig )*w10+ & + ges_qsat(ix ,iyp,k,itsig )*w01+ & + ges_qsat(ixp,iyp,k,itsig )*w11)*dtsig + & + (ges_qsat(ix ,iy ,k,itsigp)*w00+ & + ges_qsat(ixp,iy ,k,itsigp)*w10+ & + ges_qsat(ix ,iyp,k,itsigp)*w01+ & + ges_qsat(ixp,iyp,k,itsigp)*w11)*dtsigp rh(k) = q(k)/qs(k) end do endif diff --git a/src/gsi/evalqlim.f90 b/src/gsi/evalqlim.f90 index 085fa87724..a5ad216cd8 100644 --- a/src/gsi/evalqlim.f90 +++ b/src/gsi/evalqlim.f90 @@ -33,7 +33,7 @@ subroutine evalqlim(sval,pbc,rval) use kinds, only: r_kind,i_kind,r_quad use constants, only: zero,one,zero_quad use gridmod, only: lat1,lon1,nsig,istart,wgtfactlats - use jfunc, only: factqmin,factqmax + use jfunc, only: factqmin,factqmax,superfact use derivsmod, only: qgues,qsatg use mpl_allreducemod, only: mpl_allreduce use gsi_bundlemod, only: gsi_bundle @@ -83,9 +83,9 @@ subroutine evalqlim(sval,pbc,rval) endif ! Compute penalty for excess q if (q>qsatg(i,j,k)) then - term=(factqmax*wgtfactlats(ii))*(q-qsatg(i,j,k))& - /(qsatg(i,j,k)*qsatg(i,j,k)) - zbc(2) = zbc(2) + term*(q-qsatg(i,j,k)) + term=(factqmax*wgtfactlats(ii))*((q-superfact*qsatg(i,j,k))& + /qsatg(i,j,k))**2 + zbc(2) = zbc(2) + term ! Adjoint rq(i,j,k) = rq(i,j,k) + term endif diff --git a/src/gsi/genqsat.f90 b/src/gsi/genqsat.f90 index 2719ed28f8..ed0eb152e6 100644 --- a/src/gsi/genqsat.f90 +++ b/src/gsi/genqsat.f90 @@ -121,7 +121,6 @@ subroutine genqsat(qsat,tsen,prsl,lat2,lon2,nsig,ice,iderivative) end do do i=1,lat2 tdry = mint(i) - if( abs(tdry) < 1.0e-8_r_kind ) tdry = 1.0e-8_r_kind tr = ttp/tdry if (tdry >= ttp .or. .not. ice) then estmax(i) = psat * (tr**xa) * exp(xb*(one-tr)) @@ -137,7 +136,6 @@ subroutine genqsat(qsat,tsen,prsl,lat2,lon2,nsig,ice,iderivative) do k = 1,nsig do i = 1,lat2 tdry = tsen(i,j,k) - if( abs(tdry) < 1.0e-8_r_kind ) tdry = 1.0e-8_r_kind tr = ttp/tdry if (tdry >= ttp .or. .not. ice) then es = psat * (tr**xa) * exp(xb*(one-tr)) @@ -164,6 +162,7 @@ subroutine genqsat(qsat,tsen,prsl,lat2,lon2,nsig,ice,iderivative) qsat(i,j,k) = max(qmin,qsat(i,j,k)) if(iderivative > 0)then +! if(es <= esmax .and. iderivative == 2 .and. qsat(i,j,k) > qmin )then if(es <= esmax .and. iderivative == 2)then idpupdate=.true. idtupdate=.true. diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 748a3d42dd..bc26650a24 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -93,7 +93,7 @@ module gsimod pvis,pcldch,scale_cv,estvisoe,estcldchoe,vis_thres,cldch_thres,cao_check use qcmod, only: troflg,lat_c,nrand use pcpinfo, only: npredp,diag_pcp,dtphys,deltim,init_pcp - use jfunc, only: iout_iter,iguess,miter,factqmin,factqmax, & + use jfunc, only: iout_iter,iguess,miter,factqmin,factqmax,superfact,limitqobs, & factql,factqi,factqr,factqs,factqg, & factv,factl,factp,factg,factw10m,facthowv,factcldch,niter,niter_no_qc,biascor,& init_jfunc,qoption,cwoption,switch_on_derivatives,tendsflag,jiterstart,jiterend,R_option,& @@ -489,6 +489,8 @@ module gsimod ! gencode - source generation code ! factqmin - weighting factor for negative moisture constraint ! factqmax - weighting factor for supersaturated moisture constraint +! superfact- amount of supersaturation allowed 1.01 = 1% supersaturation +! limitqobs- limit q obs to be <= 100%RH based on model temperatures ! clip_supersaturation - flag to remove supersaturation during each outer loop default=.false. ! deltim - model timestep ! dtphys - physics timestep @@ -686,7 +688,7 @@ module gsimod ! NOTE: for now, if in regional mode, then iguess=-1 is forced internally. ! add use of guess file later for regional mode. - namelist/setup/gencode,factqmin,factqmax,clip_supersaturation, & + namelist/setup/gencode,factqmin,factqmax,superfact,limitqobs,clip_supersaturation, & factql,factqi,factqr,factqs,factqg, & factv,factl,factp,factg,factw10m,facthowv,factcldch,R_option,deltim,dtphys,& biascor,bcoption,diurnalbc,& @@ -1941,6 +1943,8 @@ subroutine gsimain_initialize dmesh=one factqmin=zero factqmax=zero + superfact=1._r_kind + limitqobs=.false. if (hilbert_curve) then write(6,*) 'Disabling hilbert_curve cross validation when oneobtest=.true.' hilbert_curve=.false. diff --git a/src/gsi/intjcmod.f90 b/src/gsi/intjcmod.f90 index d264cd4b5c..c0c23151ee 100644 --- a/src/gsi/intjcmod.f90 +++ b/src/gsi/intjcmod.f90 @@ -71,7 +71,7 @@ subroutine intlimq(rval,sval,itbin) ! !$$$ use gridmod, only: nsig,lat1,lon1,istart,wgtfactlats - use jfunc, only: factqmin,factqmax + use jfunc, only: factqmin,factqmax,superfact use gsi_metguess_mod, only: gsi_metguess_bundle use guess_grids, only: ges_qsat use mpimod, only: mype @@ -116,8 +116,8 @@ subroutine intlimq(rval,sval,itbin) /(ges_qsat(i,j,k,itbin)*ges_qsat(i,j,k,itbin)) ! Upper constraint limit - else if (q > ges_qsat(i,j,k,itbin)) then - rq(i,j,k) = rq(i,j,k) + (factqmax*wgtfactlats(ii))*(q-ges_qsat(i,j,k,itbin))/ & + else if (q > superfact*ges_qsat(i,j,k,itbin)) then + rq(i,j,k) = rq(i,j,k) + (factqmax*wgtfactlats(ii))*(q-superfact*ges_qsat(i,j,k,itbin))/ & (ges_qsat(i,j,k,itbin)*ges_qsat(i,j,k,itbin)) end if diff --git a/src/gsi/jfunc.f90 b/src/gsi/jfunc.f90 index 1f1adc5702..d65b4c8fd3 100644 --- a/src/gsi/jfunc.f90 +++ b/src/gsi/jfunc.f90 @@ -130,7 +130,7 @@ module jfunc public :: switch_on_derivatives,jiterend,jiterstart,jiter,iter,niter,miter public :: diurnalbc,bcoption,biascor,nval2d,xhatsave,first public :: factqmax,factqmin,clip_supersaturation,last,yhatsave,nvals_len,nval_levs,nval_levs_ens,iout_iter,nclen - public :: factql,factqi,factqr,factqs,factqg + public :: factql,factqi,factqr,factqs,factqg,superfact,limitqobs public :: niter_no_qc,print_diag_pcg,penorig,gnormorig,iguess public :: factg,factv,factp,factl,R_option,factw10m,facthowv,factcldch,diag_precon,step_start public :: pseudo_q2 @@ -139,7 +139,7 @@ module jfunc logical first,last,switch_on_derivatives,tendsflag,print_diag_pcg,tsensible,diag_precon logical clip_supersaturation,R_option - logical pseudo_q2 + logical pseudo_q2,limitqobs logical cnvw_option integer(i_kind) iout_iter,miter,iguess,nclen,qoption,cwoption integer(i_kind) jiter,jiterstart,jiterend,iter @@ -150,7 +150,7 @@ module jfunc integer(i_kind),dimension(0:50):: niter,niter_no_qc real(r_kind) factqmax,factqmin,gnormorig,penorig,biascor(2),diurnalbc,factg,factv,factp,factl,& - factw10m,facthowv,factcldch,step_start + factw10m,facthowv,factcldch,step_start,superfact real(r_kind) factql,factqi,factqr,factqs,factqg integer(i_kind) bcoption real(r_kind),allocatable,dimension(:,:):: varq @@ -202,6 +202,8 @@ subroutine init_jfunc factqmin=zero factqmax=zero + superfact=1.00_r_kind + limitqobs=.false. factql=zero factqi=zero factqr=zero diff --git a/src/gsi/observer.F90 b/src/gsi/observer.F90 index 7ac91e5eed..9b37690073 100644 --- a/src/gsi/observer.F90 +++ b/src/gsi/observer.F90 @@ -174,10 +174,10 @@ subroutine guess_init_ endif ! Read output from previous min. - if (l4dvar.and.jiterstart>1) then - else - ! If requested and if available, read guess solution. - endif +! if (l4dvar.and.jiterstart>1) then +! else +! ! If requested and if available, read guess solution. +! endif ! Generate coefficients for compact differencing if(.not.regional)then diff --git a/src/gsi/pcgsoi.f90 b/src/gsi/pcgsoi.f90 index bce3aa0835..bc840bb1ce 100644 --- a/src/gsi/pcgsoi.f90 +++ b/src/gsi/pcgsoi.f90 @@ -170,15 +170,14 @@ subroutine pcgsoi() character(5) step(2) integer(i_kind) i,istep,iobs,ii,nprt real(r_kind) stp,b,converge - real(r_kind) gsave,small_step + real(r_kind) gsave,small_step,aindex real(r_kind) gnormx,penx,penalty,penaltynew - real(r_double) pennorm - real(r_quad) zjo - real(r_quad) :: zdla - real(r_quad),dimension(4):: dprod - real(r_kind),dimension(3):: gnorm real(r_kind) :: zgini,zfini,fjcost(4),fjcostnew(4),zgend,zfend real(r_kind) :: fjcost_e + real(r_kind),dimension(3):: gnorm + real(r_double) pennorm + real(r_quad) :: zdla,zjo + real(r_quad),dimension(4):: dprod type(control_vector) :: gradx,grady,dirx,diry,ydiff,xdiff type(gsi_bundle) :: sval(nobs_bins), rval(nobs_bins) type(gsi_bundle) :: eval(ntlevs_ens) @@ -322,7 +321,7 @@ subroutine pcgsoi() end if ! 2. Multiply by background error - call multb(lanlerr,gradx,grady) + call multb(gradx,grady) if(iorthomax>0) then ! save gradients @@ -338,34 +337,23 @@ subroutine pcgsoi() ! 3. Calculate new norm of gradients and factors going into b calculation dprod(1) = qdot_prod_sub(gradx,grady) - if(iter > 0)then - gsave=gnorm(3) - if (lanlerr) then -! xdiff used as a temporary array - do i=1,nclen - xdiff%values(i)=vprecond(i)*gradx%values(i) - end do - dprod(2) = qdot_prod_sub(xdiff,grady) - call mpl_allreduce(2,qpvals=dprod) - gnorm(2)=dprod(2) - gnorm(3)=dprod(2) - else - do i=1,nclen - xdiff%values(i)=vprecond(i)*(gradx%values(i)-xdiff%values(i)) - ydiff%values(i)=vprecond(i)*(grady%values(i)-ydiff%values(i)) - end do - dprod(2) = qdot_prod_sub(xdiff,grady) - dprod(3) = qdot_prod_sub(ydiff,gradx) + if(iter > 0 .and. .not. lanlerr)then + dprod(3) = qdot_prod_sub(xdiff,grady) + dprod(4) = qdot_prod_sub(ydiff,gradx) ! xdiff used as a temporary array - do i=1,nclen - xdiff%values(i)=vprecond(i)*gradx%values(i) - end do - dprod(4) = qdot_prod_sub(xdiff,grady) - call mpl_allreduce(4,qpvals=dprod) -! Two dot products in gnorm(2) should be same, but are slightly -! different due to round off, so use average. - gnorm(2)=0.5_r_quad*(dprod(2)+dprod(3)) - gnorm(3)=dprod(4) + do i=1,nclen + xdiff%values(i)=vprecond(i)*gradx%values(i) + end do + dprod(2) = qdot_prod_sub(xdiff,grady) + call mpl_allreduce(4,qpvals=dprod) +! Two dot products in dprod(3) and dprod(4) should be same, but are slightly +! different due to round off, so use average. + gnorm(2)=dprod(2)-0.5_r_quad*(dprod(3)+dprod(4)) + gnorm(3)=dprod(2) + if(mype == 0)then + aindex=abs(dprod(3)/dprod(2)) + write(iout_iter,*) 'NL Index ',aindex + if(aindex > 0.5 .or. print_verbose) write(iout_iter,*) 'NL Values ', dprod(3),dprod(2) end if else ! xdiff used as a temporary array @@ -375,38 +363,36 @@ subroutine pcgsoi() dprod(2) = qdot_prod_sub(xdiff,grady) call mpl_allreduce(2,qpvals=dprod) if(print_diag_pcg) call prt_control_norms(grady,'grady') + gnorm(2)=dprod(2) gnorm(3)=dprod(2) end if gnorm(1)=dprod(1) - if(mype == 0)write(iout_iter,*)'Minimization iteration',iter + if(mype == 0)write(iout_iter,*)'Minimization iteration',iter ! 4. Calculate b and new search direction b=zero if (.not. restart .or. iter > 0) then - if (gsave>1.e-16_r_kind .and. iter>0) b=gnorm(2)/gsave - if (b7.0_r_kind) then - if (mype==0) then - if (iout_6) write(6,105) gnorm(2),gsave,b - write(iout_iter,105) gnorm(2),gsave,b + if (iter > 1 .or. .not. read_success)then + if (gsave>1.e-16_r_kind) b=gnorm(2)/gsave + if (b10.0_r_kind) then + if (mype==0) then + if (iout_6) write(6,105) gnorm(2),gsave,b + write(iout_iter,105) gnorm(2),gsave,b + endif + b=zero endif - b=zero - endif - if (mype==0 .and. print_verbose) write(6,888)'pcgsoi: gnorm(1:3),b=',gnorm,b - if(read_success .and. iter == 1)b=zero - - if (.not. lanlerr) then - do i=1,nclen - xdiff%values(i)=gradx%values(i) - ydiff%values(i)=grady%values(i) - end do + if (mype==0 .and. print_verbose) write(6,888)'pcgsoi: gnorm(1:3),b=',gnorm,b end if -! Calculate new search direction + do i=1,nclen - dirx%values(i)=-vprecond(i)*grady%values(i)+b*dirx%values(i) - diry%values(i)=-vprecond(i)*gradx%values(i)+b*diry%values(i) +! Calculate new search direction + ydiff%values(i)=vprecond(i)*grady%values(i) + dirx%values(i)=-ydiff%values(i)+b*dirx%values(i) + xdiff%values(i)=vprecond(i)*gradx%values(i) + diry%values(i)=-xdiff%values(i)+b*diry%values(i) end do else ! If previous solution available, transfer into local arrays. @@ -416,9 +402,11 @@ subroutine pcgsoi() end do call read_guess_solution(diry,mype,read_success) ! Multiply by background error - call multb(lanlerr,diry,dirx) + call multb(diry,dirx) + restart=.false. endif - + gsave=gnorm(3) + ! 5. Calculate stepsize and update solution ! Convert search direction from control space to physical space do ii=1,nobs_bins @@ -568,7 +556,7 @@ subroutine pcgsoi() end do ! Multiply by background error - call multb(lanlerr,gradx,grady) + call multb(gradx,grady) ! Print final Jo table zgend=dot_product(gradx,grady,r_quad) @@ -630,8 +618,7 @@ subroutine pcgsoi() ! if (mype==0) write(6,*)'pcgsoi: Updating guess' if(iwrtinc<=0) call update_guess(sval,sbias) -! cloud analysis after iteration -! if(jiter == miter .and. i_gsdcldanal_type==1) then +! gsd cloud analysis after iteration if(jiter == miter) then if(i_gsdcldanal_type==2) then call gsdcloudanalysis4nmmb(mype) @@ -833,7 +820,7 @@ subroutine periodic_(gradx) end subroutine periodic_ -subroutine multb(lanlerr,vec1,vec2) +subroutine multb(vec1,vec2) !$$$ subprogram documentation block ! . . . . ! subprogram: multb multiply vec1 by background error to equal vec2 @@ -863,7 +850,6 @@ subroutine multb(lanlerr,vec1,vec2) type(control_vector),intent(inout) :: vec1 type(control_vector),intent(inout) :: vec2 - logical ,intent(in ) :: lanlerr if(periodic)call periodic_(vec1) ! start by setting vec2=vec1 and then operate on vec2 (unless gram_schmidt) @@ -892,7 +878,7 @@ end subroutine multb subroutine c2s(hat,val,bias,llprt,ltest) !$$$ subprogram documentation block ! . . . . -! subprogram: multb control2state for all options +! subprogram: c2s control2state for all options ! prgmmr: derber ! ! abstract: generalized control2state @@ -958,7 +944,7 @@ end subroutine c2s subroutine c2s_ad(hat,val,bias,llprt) !$$$ subprogram documentation block ! . . . . -! subprogram: multb control2state for all options +! subprogram: c2s_ad adjoint of control2state for all options ! prgmmr: derber ! ! abstract: generalized control2state diff --git a/src/gsi/prt_guess.f90 b/src/gsi/prt_guess.f90 index ad5cb80025..e1e251195a 100644 --- a/src/gsi/prt_guess.f90 +++ b/src/gsi/prt_guess.f90 @@ -165,7 +165,7 @@ subroutine prt_guess(sgrep) zloc(nvars+7) = minval(ges_cwmr_it(2:lat1+1,2:lon1+1,1:nsig)) zloc(nvars+8) = minval(ges_cf_it(2:lat1+1,2:lon1+1,1:nsig)) zloc(nvars+9) = minval(ges_div_it(2:lat1+1,2:lon1+1,1:nsig)) - zloc(nvars+10) = minval(ges_vor_it(2:lat1+1,2:lon1+1,1:nsig)) + zloc(nvars+10) = minval(ges_vor_it(2:lat1+1,2:lon1+1,1:nsig)) zloc(nvars+11) = minval(ges_prsl (2:lat1+1,2:lon1+1,1:nsig,ntsig)) zloc(nvars+12) = minval(ges_ps_it (2:lat1+1,2:lon1+1 )) zloc(nvars+13) = minval(sfct (2:lat1+1,2:lon1+1, ntsfc)) @@ -178,7 +178,7 @@ subroutine prt_guess(sgrep) zloc(2*nvars+7) = maxval(ges_cwmr_it(2:lat1+1,2:lon1+1,1:nsig)) zloc(2*nvars+8) = maxval(ges_cf_it(2:lat1+1,2:lon1+1,1:nsig)) zloc(2*nvars+9) = maxval(ges_div_it(2:lat1+1,2:lon1+1,1:nsig)) - zloc(2*nvars+10) = maxval(ges_vor_it(2:lat1+1,2:lon1+1,1:nsig)) + zloc(2*nvars+10) = maxval(ges_vor_it(2:lat1+1,2:lon1+1,1:nsig)) zloc(2*nvars+11) = maxval(ges_prsl (2:lat1+1,2:lon1+1,1:nsig,ntsig)) zloc(2*nvars+12) = maxval(ges_ps_it (2:lat1+1,2:lon1+1 )) zloc(2*nvars+13) = maxval(sfct (2:lat1+1,2:lon1+1, ntsfc)) @@ -188,8 +188,8 @@ subroutine prt_guess(sgrep) ! Gather contributions - call mpi_allgather(zloc,3*nvars+3,mpi_rtype, & - & zall,3*nvars+3,mpi_rtype, mpi_comm_world,ierror) + call mpi_gather(zloc,3*nvars+3,mpi_rtype, & + & zall,3*nvars+3,mpi_rtype,0, mpi_comm_world,ierror) if (mype==0) then zmin=zero diff --git a/src/gsi/read_guess.F90 b/src/gsi/read_guess.F90 index 364d4e305b..ecd85aa983 100644 --- a/src/gsi/read_guess.F90 +++ b/src/gsi/read_guess.F90 @@ -89,7 +89,7 @@ subroutine read_guess(iyear,month,idd,mype) !$$$ use kinds, only: r_kind,i_kind - use jfunc, only: bcoption,clip_supersaturation + use jfunc, only: bcoption,clip_supersaturation,superfact use guess_grids, only: nfldsig,ges_tsen,load_prsges,load_geop_hgt,ges_prsl,& ges_tsen1, geop_hgti, ges_geopi, ges_q1 use m_gsiBiases,only : bkg_bias_correction,nbc @@ -250,7 +250,7 @@ subroutine read_guess(iyear,month,idd,mype) do k=1,nsig do j=1,lon2 do i=1,lat2 - satval = min(ges_q(i,j,k),satq(i,j,k)) + satval = min(ges_q(i,j,k),superfact*satq(i,j,k)) satval = max(qmin,satval) ges_q(i,j,k) = satval ges_tsen(i,j,k,it)= ges_tv(i,j,k)/(one+fv*ges_q(i,j,k)) diff --git a/src/gsi/setupcldtot.F90 b/src/gsi/setupcldtot.F90 index fc3978e3d2..565075bec3 100755 --- a/src/gsi/setupcldtot.F90 +++ b/src/gsi/setupcldtot.F90 @@ -109,7 +109,6 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di external:: tintrp2a1,tintrp2a11 external:: tintrp31,tintrp3 external:: grdcrd1 - external:: genqsat external:: stop2 ! Declare local variables diff --git a/src/gsi/setupq.f90 b/src/gsi/setupq.f90 index 77c43ba121..8189b0bcb5 100644 --- a/src/gsi/setupq.f90 +++ b/src/gsi/setupq.f90 @@ -149,13 +149,13 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav use nc_diag_read_mod, only: nc_diag_read_init, nc_diag_read_get_dim, nc_diag_read_close use gsi_4dvar, only: nobs_bins,hr_obsbin,min_offset use oneobmod, only: oneobtest,maginnov,magoberr - use guess_grids, only: ges_lnprsl,hrdifsig,nfldsig,ges_tsen,ges_prsl,pbl_height + use guess_grids, only: ges_lnprsl,hrdifsig,nfldsig,ges_tsen,ges_prsl,pbl_height,ges_qsat use gridmod, only: lat2,lon2,nsig,get_ijk,twodvar_regional use constants, only: zero,one,r1000,r10,r100 use constants, only: huge_single,wgtlim,three use constants, only: tiny_r_kind,five,half,two,huge_r_kind,r0_01 use qcmod, only: npres_print,ptopq,pbotq,dfact,dfact1,njqc,vqc,nvqc - use jfunc, only: jiter,last,jiterstart,miter + use jfunc, only: jiter,last,jiterstart,miter,superfact,limitqobs use convinfo, only: nconvtype,cermin,cermax,cgross,cvar_b,cvar_pg,ictype use convinfo, only: ibeta,ikapa use convinfo, only: icsubtype @@ -417,7 +417,7 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav iderivative=0 do jj=1,nfldsig call genqsat(qg(1,1,1,jj),ges_tsen(1,1,1,jj),ges_prsl(1,1,1,jj),lat2,lon2,nsig,ice,iderivative) - call genqsat(qg2m(1,1,jj),ges_tsen(1,1,1,jj),ges_prsl(1,1,1,jj),lat2,lon2, 1,ice,iderivative) + qg2m(:,:,jj)=qg(:,:,1,jj) end do @@ -516,9 +516,15 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! Scale errors by guess saturation q - call tintrp31(qg,qsges,dlat,dlon,dpres,dtime,hrdifsig,& + qob = data(iqob,i) + if(limitqobs) then + call tintrp31(ges_qsat,qsges,dlat,dlon,dpres,dtime,hrdifsig,& mype,nfldsig) + qob=min(qob,superfact*qsges) + end if + call tintrp31(qg,qsges,dlat,dlon,dpres,dtime,hrdifsig,& + mype,nfldsig) ! Interpolate 2-m qs to obs locations/times if((i_use_2mq4b > 0) .and. ((itype > 179 .and. itype < 190) .or. itype == 199) & .and. .not.twodvar_regional)then @@ -527,7 +533,6 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! Load obs error and value into local variables obserror = max(cermin(ikx)*r0_01,min(cermax(ikx)*r0_01,data(ier,i))) - qob = data(iqob,i) rmaxerr=rmaxerr*qsges rmaxerr=max(small2,rmaxerr) diff --git a/src/gsi/setuprhsall.f90 b/src/gsi/setuprhsall.f90 index 8e34876201..3efcb69859 100644 --- a/src/gsi/setuprhsall.f90 +++ b/src/gsi/setuprhsall.f90 @@ -150,7 +150,6 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass) use lag_fields, only: lag_presetup,lag_state_write,lag_state_read,lag_destroy_uv use mpeu_util, only: getindex use mpl_allreducemod, only: mpl_allreduce - use berror, only: reset_predictors_var use rapidrefresh_cldsurf_mod, only: l_PBL_pseudo_SurfobsT,l_PBL_pseudo_SurfobsQ,& l_PBL_pseudo_SurfobsUV,i_gsdcldanal_type,& i_cloud_q_innovation @@ -579,10 +578,6 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass) call mpl_allreduce(npredt,ntail,rstats_t) end if - if (newpc4pred .or. aircraft_t_bc_pof .or. aircraft_t_bc) then - call reset_predictors_var - end if - ! Collect satellite and precip. statistics call mpi_reduce(aivals,aivals1,size(aivals1),mpi_rtype,mpi_sum,mype_rad, & mpi_comm_world,ierror) diff --git a/src/gsi/setupw.f90 b/src/gsi/setupw.f90 index aa6fe55094..11a5727ce6 100644 --- a/src/gsi/setupw.f90 +++ b/src/gsi/setupw.f90 @@ -253,9 +253,9 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav real(r_double) rstation_id real(r_kind) qcu,qcv,trop5,tfact,fact real(r_kind) scale,ratio,obserror,obserrlm - real(r_kind) residual,ressw,ress,val,vals,val2,valqc2,dudiff,dvdiff - real(r_kind) valqc,valu,valv,dx10,rlow,rhgh,drpx,prsfc,var_jb - real(r_kind) cg_t,cvar,wgt,term,rat_err2,qcgross + real(r_kind) residual,ressw,ress,vals,val2,valqc2,dudiff,dvdiff,rat_err2u + real(r_kind) valqc,valu,valv,dx10,rlow,rhgh,drpx,prsfc,var_jb,rat_err2v + real(r_kind) cg_t,cvar,wgt,term,rat_err2,qcgross,valqcu,valqcv real(r_kind) presw,factw,dpres,ugesin,vgesin,rwgt,dpressave real(r_kind) sfcchk,prsln2,error,dtime,dlon,dlat,r0_001,rsig,thirty,rsigp real(r_kind) ratio_errors,goverrd,spdges,spdob,ten,psges,zsges @@ -1222,8 +1222,6 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ratio_errors=ratio_errors*sqrt(data(ihil,i)) ! hilbert weight ! Compute penalty terms (linear & nonlinear qc). if(luse(i))then - val = valu*valu+valv*valv - vals=sqrt(val) if(vqc) then cg_t=cvar_b(ikx) cvar=cvar_pg(ikx) @@ -1238,17 +1236,23 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ibb=0 ikk=0 endif + vals=valu call vqc_setup(vals,ratio_errors,error,cvar,& - cg_t,ibb,ikk,var_jb,rat_err2,wgt,valqc) - rwgt = wgt/wgtlim + cg_t,ibb,ikk,var_jb,rat_err2u,wgt,valqcu) + rwgt = 0.5_r_kind*wgt/wgtlim + vals=valv + call vqc_setup(vals,ratio_errors,error,cvar,& + cg_t,ibb,ikk,var_jb,rat_err2v,wgt,valqcv) + rwgt = rwgt+0.5_r_kind*wgt/wgtlim + valqc=half*(valqcu+valqcv) ! Accumulate statistics for obs belonging to this task if (muse(i)) then if(rwgt < one) awork(21) = awork(21)+one jsig = dpres jsig=max(1,min(jsig,nsig)) - awork(4*nsig+jsig+100)=awork(4*nsig+jsig+100)+valu*valu*rat_err2 - awork(5*nsig+jsig+100)=awork(5*nsig+jsig+100)+valv*valv*rat_err2 + awork(4*nsig+jsig+100)=awork(4*nsig+jsig+100)+valu*valu*rat_err2u + awork(5*nsig+jsig+100)=awork(5*nsig+jsig+100)+valv*valv*rat_err2v awork(6*nsig+jsig+100)=awork(6*nsig+jsig+100)+one awork(3*nsig+jsig+100)=awork(3*nsig+jsig+100)+valqc end if @@ -1257,8 +1261,7 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! as a function of observation type. ress = scale*sqrt(dudiff**2+dvdiff**2) ressw = ress*ress - val2 = half*(valu*valu+valv*valv) - valqc2 = half*valqc + val2 = half*(valu*valu*rat_err2u+valv*valv*rat_err2v) nn=1 if (.not. muse(i)) then nn=2 @@ -1269,8 +1272,8 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav bwork(k,ikx,1,nn) = bwork(k,ikx,1,nn)+one ! count bwork(k,ikx,2,nn) = bwork(k,ikx,2,nn)+spdb ! speed bias bwork(k,ikx,3,nn) = bwork(k,ikx,3,nn)+ressw ! (o-g)**2 - bwork(k,ikx,4,nn) = bwork(k,ikx,4,nn)+val2*rat_err2 ! penalty - bwork(k,ikx,5,nn) = bwork(k,ikx,5,nn)+valqc2 ! nonlin qc penalty + bwork(k,ikx,4,nn) = bwork(k,ikx,4,nn)+val2 ! penalty + bwork(k,ikx,5,nn) = bwork(k,ikx,5,nn)+valqc ! nonlin qc penalty end if end do diff --git a/src/gsi/stpcalc.f90 b/src/gsi/stpcalc.f90 index 8b67e35742..80fac64d61 100644 --- a/src/gsi/stpcalc.f90 +++ b/src/gsi/stpcalc.f90 @@ -280,7 +280,7 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & real(r_kind),dimension(4)::sges real(r_kind),dimension(ioutpen):: outpen,outstp logical :: cxterm,change_dels,ifound - logical :: print_verbose + logical :: print_verbose,pjcalc !************************************************************************************ @@ -296,6 +296,8 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & outpen = zero nsteptot=0 istp_use=0 + kprt=3 + pjcalc=.false. pj=zero_quad ! Begin calculating contributions to penalty and stepsize for various terms @@ -386,12 +388,13 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & pstart=zero_quad pbc=zero_quad + if(iter == 0 .and. kprt >= 2)pjcalc=.true. ! penalty, b and c for background terms pstart(1,1) = qdot_prod_sub(xhatsave,yhatsave) - pj(1,1)=pstart(1,1) + if(pjcalc)pj(1,1)=pstart(1,1) ! two terms in next line should be the same, but roundoff makes average more accurate. @@ -406,7 +409,7 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & if (ljcdfi .and. nobs_bins>1) then call stpjcdfi(dval,sval,pstart(1,2),pstart(2,2),pstart(3,2)) - pj(2,1)=pstart(1,2) + if(pjcalc)pj(2,1)=pstart(1,2) end if ! Penalty, b, c for dry pressure @@ -416,13 +419,15 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & else call stpjcpdry(dval,sval,pstart(1,3),pstart(2,3),pstart(3,3),nobs_bins) end if - pj(3,1)=pstart(1,3) + if(pjcalc)pj(3,1)=pstart(1,3) end if ! iterate over number of stepsize iterations (istp_iter - currently set to a maximum of 5) dels = one_tenth_quad stepsize: do ii=1,istp_iter + pjcalc=.false. + if(iter == 0 .and. kprt >= 2 .and. ii == 1)pjcalc=.true. iis=ii ! Delta stepsize change_dels=.true. @@ -457,11 +462,76 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & ! penalties for moisture constraint if(.not. ltlint)then + if(.not.ljc4tlevs) then + call stplimq(dval(ibin_anl),sval(ibin_anl),sges,pbc(1,4),pbc(1,5),nstep,ntguessig) + if(pjcalc)then + pj(4,1)=pbc(1,4)+pbc(ipenloc,4) + pj(5,1)=pbc(1,5)+pbc(ipenloc,5) + end if + else + do ibin=1,nobs_bins + if (nobs_bins /= nfldsig) then + it=ntguessig + else + it=ibin + end if + call stplimq(dval(ibin),sval(ibin),sges,pbcqmin(1,ibin),pbcqmax(1,ibin),nstep,it) + end do + pbc(:,4)=zero_quad + pbc(:,5)=zero_quad + do ibin=1,nobs_bins + do j=1,nstep + pbc(j,4) = pbc(j,4)+pbcqmin(j,ibin) + pbc(j,5) = pbc(j,5)+pbcqmax(j,ibin) + end do + end do + if(pjcalc)then + do ibin=1,nobs_bins + pj(4,ibin)=pj(4,ibin)+pbcqmin(1,ibin)+pbcqmin(ipenloc,ibin) + pj(5,ibin)=pj(5,ibin)+pbcqmax(1,ibin)+pbcqmax(ipenloc,ibin) + end do + end if + end if +! penalties for gust constraint + if(getindex(cvars2d,'gust')>0) & + call stplimg(dval(1),sval(1),sges,pbc(1,6),nstep) + if(pjcalc)pj(6,1)=pbc(1,6)+pbc(ipenloc,6) + +! penalties for vis constraint + if(getindex(cvars2d,'vis')>0) & + call stplimv(dval(1),sval(1),sges,pbc(1,7),nstep) + if(pjcalc)pj(7,1)=pbc(1,7)+pbc(ipenloc,7) + +! penalties for pblh constraint + if(getindex(cvars2d,'pblh')>0) & + call stplimp(dval(1),sval(1),sges,pbc(1,8),nstep) + if(pjcalc)pj(8,1)=pbc(1,8)+pbc(ipenloc,8) + +! penalties for wspd10m constraint + if(getindex(cvars2d,'wspd10m')>0) & + call stplimw10m(dval(1),sval(1),sges,pbc(1,9),nstep) + if(pjcalc)pj(9,1)=pbc(1,9)+pbc(ipenloc,9) + +! penalties for howv constraint + if(getindex(cvars2d,'howv')>0) & + call stplimhowv(dval(1),sval(1),sges,pbc(1,10),nstep) + if(pjcalc)pj(10,1)=pbc(1,10)+pbc(ipenloc,10) + +! penalties for lcbas constraint + if(getindex(cvars2d,'lcbas')>0) & + call stpliml(dval(1),sval(1),sges,pbc(1,11),nstep) + if(pjcalc)pj(11,1)=pbc(1,11)+pbc(ipenloc,11) + +! penalties for cldch constraint + if(getindex(cvars2d,'cldch')>0) & + call stplimcldch(dval(1),sval(1),sges,pbc(1,12),nstep) + if(pjcalc)pj(12,1)=pbc(1,12)+pbc(ipenloc,12) + if (ljclimqc) then - if (getindex(cvars3d,'ql')>0) then + if (getindex(cvars3d,'ql')>0) then if(.not.ljc4tlevs) then call stplimqc(dval(ibin_anl),sval(ibin_anl),sges,pbc(1,13),nstep,ntguessig,'ql') - if(ii == 1) pj(13,1)=pbc(1,13)+pbc(ipenloc,13) + if(pjcalc) pj(13,1)=pbc(1,13)+pbc(ipenloc,13) else do ibin=1,nobs_bins if (nobs_bins /= nfldsig) then @@ -477,17 +547,17 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & pbc(j,13) = pbc(j,13)+pbcql(j,ibin) end do end do - if(ii == 1)then + if(pjcalc)then do ibin=1,nobs_bins pj(13,ibin)=pj(13,ibin)+pbcql(1,ibin)+pbcql(ipenloc,ibin) end do end if end if - end if - if (getindex(cvars3d,'qi')>0) then + end if + if (getindex(cvars3d,'qi')>0) then if(.not.ljc4tlevs) then call stplimqc(dval(ibin_anl),sval(ibin_anl),sges,pbc(1,14),nstep,ntguessig,'qi') - if(ii == 1) pj(14,1)=pbc(1,14)+pbc(ipenloc,14) + if(pjcalc) pj(14,1)=pbc(1,14)+pbc(ipenloc,14) else do ibin=1,nobs_bins if (nobs_bins /= nfldsig) then @@ -503,17 +573,17 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & pbc(j,14) = pbc(j,14)+pbcqi(j,ibin) end do end do - if(ii == 1)then + if(pjcalc)then do ibin=1,nobs_bins pj(14,ibin)=pj(14,ibin)+pbcqi(1,ibin)+pbcqi(ipenloc,ibin) end do end if end if - end if - if (getindex(cvars3d,'qr')>0) then + end if + if (getindex(cvars3d,'qr')>0) then if(.not.ljc4tlevs) then call stplimqc(dval(ibin_anl),sval(ibin_anl),sges,pbc(1,15),nstep,ntguessig,'qr') - if(ii == 1) pj(15,1)=pbc(1,15)+pbc(ipenloc,15) + if(pjcalc) pj(15,1)=pbc(1,15)+pbc(ipenloc,15) else do ibin=1,nobs_bins if (nobs_bins /= nfldsig) then @@ -529,17 +599,17 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & pbc(j,15) = pbc(j,15)+pbcqr(j,ibin) end do end do - if(ii == 1)then + if(pjcalc)then do ibin=1,nobs_bins pj(15,ibin)=pj(15,ibin)+pbcqr(1,ibin)+pbcqr(ipenloc,ibin) end do end if end if - end if - if (getindex(cvars3d,'qs')>0) then + end if + if (getindex(cvars3d,'qs')>0) then if(.not.ljc4tlevs) then call stplimqc(dval(ibin_anl),sval(ibin_anl),sges,pbc(1,16),nstep,ntguessig,'qs') - if(ii == 1) pj(16,1)=pbc(1,16)+pbc(ipenloc,16) + if(pjcalc) pj(16,1)=pbc(1,16)+pbc(ipenloc,16) else do ibin=1,nobs_bins if (nobs_bins /= nfldsig) then @@ -555,17 +625,17 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & pbc(j,16) = pbc(j,16)+pbcqs(j,ibin) end do end do - if(ii == 1)then + if(pjcalc)then do ibin=1,nobs_bins pj(16,ibin)=pj(16,ibin)+pbcqs(1,ibin)+pbcqs(ipenloc,ibin) end do end if end if - end if - if (getindex(cvars3d,'qg')>0) then + end if + if (getindex(cvars3d,'qg')>0) then if(.not.ljc4tlevs) then call stplimqc(dval(ibin_anl),sval(ibin_anl),sges,pbc(1,17),nstep,ntguessig,'qg') - if(ii == 1) pj(17,1)=pbc(1,17)+pbc(ipenloc,17) + if(pjcalc) pj(17,1)=pbc(1,17)+pbc(ipenloc,17) else do ibin=1,nobs_bins if (nobs_bins /= nfldsig) then @@ -581,78 +651,14 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & pbc(j,17) = pbc(j,17)+pbcqg(j,ibin) end do end do - if(ii == 1)then + if(pjcalc)then do ibin=1,nobs_bins pj(17,ibin)=pj(17,ibin)+pbcqg(1,ibin)+pbcqg(ipenloc,ibin) end do end if end if - end if + end if end if ! ljclimqc - if(.not.ljc4tlevs) then - call stplimq(dval(ibin_anl),sval(ibin_anl),sges,pbc(1,4),pbc(1,5),nstep,ntguessig) - if(ii == 1)then - pj(4,1)=pbc(1,4)+pbc(ipenloc,4) - pj(5,1)=pbc(1,5)+pbc(ipenloc,5) - end if - else - do ibin=1,nobs_bins - if (nobs_bins /= nfldsig) then - it=ntguessig - else - it=ibin - end if - call stplimq(dval(ibin),sval(ibin),sges,pbcqmin(1,ibin),pbcqmax(1,ibin),nstep,it) - end do - pbc(:,4)=zero_quad - pbc(:,5)=zero_quad - do ibin=1,nobs_bins - do j=1,nstep - pbc(j,4) = pbc(j,4)+pbcqmin(j,ibin) - pbc(j,5) = pbc(j,5)+pbcqmax(j,ibin) - end do - end do - if(ii == 1)then - do ibin=1,nobs_bins - pj(4,ibin)=pj(4,ibin)+pbcqmin(1,ibin)+pbcqmin(ipenloc,ibin) - pj(5,ibin)=pj(5,ibin)+pbcqmax(1,ibin)+pbcqmax(ipenloc,ibin) - end do - end if - end if -! penalties for gust constraint - if(getindex(cvars2d,'gust')>0) & - call stplimg(dval(1),sval(1),sges,pbc(1,6),nstep) - if(ii == 1)pj(6,1)=pbc(1,6)+pbc(ipenloc,6) - -! penalties for vis constraint - if(getindex(cvars2d,'vis')>0) & - call stplimv(dval(1),sval(1),sges,pbc(1,7),nstep) - if(ii == 1)pj(7,1)=pbc(1,7)+pbc(ipenloc,7) - -! penalties for pblh constraint - if(getindex(cvars2d,'pblh')>0) & - call stplimp(dval(1),sval(1),sges,pbc(1,8),nstep) - if(ii == 1)pj(8,1)=pbc(1,8)+pbc(ipenloc,8) - -! penalties for wspd10m constraint - if(getindex(cvars2d,'wspd10m')>0) & - call stplimw10m(dval(1),sval(1),sges,pbc(1,9),nstep) - if(ii == 1)pj(9,1)=pbc(1,9)+pbc(ipenloc,9) - -! penalties for howv constraint - if(getindex(cvars2d,'howv')>0) & - call stplimhowv(dval(1),sval(1),sges,pbc(1,10),nstep) - if(ii == 1)pj(10,1)=pbc(1,10)+pbc(ipenloc,10) - -! penalties for lcbas constraint - if(getindex(cvars2d,'lcbas')>0) & - call stpliml(dval(1),sval(1),sges,pbc(1,11),nstep) - if(ii == 1)pj(11,1)=pbc(1,11)+pbc(ipenloc,11) - -! penalties for cldch constraint - if(getindex(cvars2d,'cldch')>0) & - call stplimcldch(dval(1),sval(1),sges,pbc(1,12),nstep) - if(ii == 1)pj(12,1)=pbc(1,12)+pbc(ipenloc,12) end if ! penalties for Jo @@ -667,7 +673,7 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & end do end do enddo - if(ii == 1)then + if(pjcalc)then do ibin=1,size(pbcjoi,3) do j=1,size(pbcjoi,2) pj(n0+j,ibin)=pj(n0+j,ibin)+pbcjoi(ipenloc,j,ibin)+pbcjoi(1,j,ibin) @@ -866,7 +872,6 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & exit stepsize end if end do stepsize - kprt=3 if(kprt >= 2 .and. iter == 0)then call mpl_allreduce(ipen,nobs_bins,pj) if(mype == 0)call prnt_j(pj,n0,ipen,kprt) diff --git a/src/gsi/stpjcmod.f90 b/src/gsi/stpjcmod.f90 index 289cc7672a..2c811912a0 100644 --- a/src/gsi/stpjcmod.f90 +++ b/src/gsi/stpjcmod.f90 @@ -77,7 +77,7 @@ subroutine stplimq(rval,sval,sges,outmin,outmax,nstep,itbin) ! !$$$ use gridmod, only: lat1,lon1,nsig,istart,wgtfactlats - use jfunc, only: factqmin,factqmax + use jfunc, only: factqmin,factqmax,superfact use guess_grids, only: ges_qsat use mpimod, only: mype implicit none @@ -125,8 +125,8 @@ subroutine stplimq(rval,sval,sges,outmin,outmax,nstep,itbin) /(ges_qsat(i,j,k,itbin)*ges_qsat(i,j,k,itbin)) else if(qx > ges_qsat(i,j,k,itbin))then - outmax(kk)=outmax(kk)+(factqmax*wgtfactlats(ii))*(qx-ges_qsat(i,j,k,itbin))* & - (qx-ges_qsat(i,j,k,itbin))/(ges_qsat(i,j,k,itbin)*ges_qsat(i,j,k,itbin)) + outmax(kk)=outmax(kk)+(factqmax*wgtfactlats(ii))*(qx-superfact*ges_qsat(i,j,k,itbin))**2 & + /ges_qsat(i,j,k,itbin)**2 end if end if end do @@ -143,9 +143,9 @@ subroutine stplimq(rval,sval,sges,outmin,outmax,nstep,itbin) if(q < zero)then outmin(1)=outmin(1)+(factqmin*wgtfactlats(ii))*q*q/(ges_qsat(i,j,k,itbin)*ges_qsat(i,j,k,itbin)) else - if(q > ges_qsat(i,j,k,itbin))then - outmax(1)=outmax(1)+(factqmax*wgtfactlats(ii))*(q-ges_qsat(i,j,k,itbin))*& - (q-ges_qsat(i,j,k,itbin))/(ges_qsat(i,j,k,itbin)*ges_qsat(i,j,k,itbin)) + if(q > superfact*ges_qsat(i,j,k,itbin))then + outmax(1)=outmax(1)+(factqmax*wgtfactlats(ii))*(q-superfact*ges_qsat(i,j,k,itbin))**2 & + /ges_qsat(i,j,k,itbin)**2 end if end if end do @@ -242,7 +242,7 @@ subroutine stplimqc(rval,sval,sges,out,nstep,itbin,cldtype) endif if (mype==0) write(6,*)'stplimqc: factqc = ', factqc if (mype==0) write(6,*)'stplimqc: ier ier1 = ', ier, ier1 - if ( factqc==0 ) return + if ( factqc <= zero) return if ( ier/=0 .or. ier1/=0 ) return ! Loop over interior of subdomain diff --git a/src/gsi/stpw.f90 b/src/gsi/stpw.f90 index ee1569d731..423019680d 100644 --- a/src/gsi/stpw.f90 +++ b/src/gsi/stpw.f90 @@ -105,7 +105,7 @@ subroutine stpw(whead,rval,sval,out,sges,nstep) real(r_kind) valu,facu,valv,facv,w1,w2,w3,w4,w5,w6,w7,w8 real(r_kind) cg_t,t_pg,var_jb real(r_kind) uu,vv - real(r_kind),dimension(max(1,nstep))::pen + real(r_kind),dimension(max(1,nstep))::pen,penu,penv real(r_kind),pointer,dimension(:):: ru,rv,su,sv type(wNode), pointer :: wptr @@ -157,41 +157,39 @@ subroutine stpw(whead,rval,sval,out,sges,nstep) do kk=1,nstep uu=facu+sges(kk)*valu vv=facv+sges(kk)*valv - pen(kk)= (uu*uu+vv*vv)*wptr%err2 + penu(kk)= (uu*uu)*wptr%err2 + penv(kk)= (vv*vv)*wptr%err2 end do else - pen(1)= (wptr%ures*wptr%ures+wptr%vres*wptr%vres)*wptr%err2 + penu(1)= (wptr%ures*wptr%ures)*wptr%err2 + penv(1)= (wptr%vres*wptr%vres)*wptr%err2 end if ! Modify penalty term if nonlinear QC - if (vqc .and. nlnqc_iter .and. wptr%pg > tiny_r_kind .and. & + t_pg=zero + cg_t=zero + var_jb=zero + ibb=0 + ikk=0 + if (vqc .and. nlnqc_iter .and. wptr%pg > tiny_r_kind .and. & wptr%b > tiny_r_kind) then t_pg=wptr%pg*varqc_iter cg_t=cg_term/wptr%b - else - t_pg=zero - cg_t=zero - endif - + else if(njqc .and. wptr%jb > tiny_r_kind .and. wptr%jb <10.0_r_kind) then ! for Dr. Jim purser' non liear quality control - if(njqc .and. wptr%jb > tiny_r_kind .and. wptr%jb <10.0_r_kind) then var_jb =wptr%jb - else - var_jb=zero - endif + else if(nvqc .and. wptr%ib >0) then ! mix model VQC - if(nvqc .and. wptr%ib >0) then ibb=wptr%ib ikk=wptr%ik - else - ibb=0 - ikk=0 endif - - call vqc_stp(pen,nstep,t_pg,cg_t,var_jb,ibb,ikk) - + call vqc_stp(penu,nstep,t_pg,cg_t,var_jb,ibb,ikk) + call vqc_stp(penv,nstep,t_pg,cg_t,var_jb,ibb,ikk) + do kk=1,nstep + pen(kk)=penu(kk)+penv(kk) + end do out(1) = out(1)+pen(1)*wptr%raterr2 do kk=2,nstep out(kk) = out(kk)+(pen(kk)-pen(1))*wptr%raterr2 diff --git a/src/gsi/update_guess.f90 b/src/gsi/update_guess.f90 index 041ad50e1d..a66495bc27 100644 --- a/src/gsi/update_guess.f90 +++ b/src/gsi/update_guess.f90 @@ -113,7 +113,7 @@ subroutine update_guess(sval,sbias) use mpimod, only: mype use constants, only: zero,one,fv,max_varname_length,qmin,qcmin,tgmin,& r100,one_tenth,tiny_r_kind - use jfunc, only: iout_iter,bcoption,tsensible,clip_supersaturation + use jfunc, only: iout_iter,bcoption,tsensible,clip_supersaturation,superfact use gridmod, only: lat2,lon2,nsig,& regional,twodvar_regional,regional_ozone,& l_reg_update_hydro_delz @@ -269,7 +269,7 @@ subroutine update_guess(sval,sbias) call gsi_bundlegetpointer (gsi_metguess_bundle(it),guess(ic),ptr3dges,istatus) if (trim(guess(ic))=='q') then call upd_positive_fldr3_(ptr3dges,ptr3dinc, qmin) - if(clip_supersaturation) ptr3dges(:,:,:) = min(ptr3dges(:,:,:),ges_qsat(:,:,:,it)) + if(clip_supersaturation) ptr3dges(:,:,:) = min(ptr3dges(:,:,:),superfact*ges_qsat(:,:,:,it)) cycle endif if (trim(guess(ic))=='oz') then From 52c5ae6802c5bc933b6b65b0a3512e044b3d167e Mon Sep 17 00:00:00 2001 From: "john.derber" Date: Fri, 5 Nov 2021 02:16:05 +0000 Subject: [PATCH 3/5] fix setupw --- fix | 2 +- src/gsi/setupw.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/fix b/fix index 13ad7b5abb..dc537d7106 160000 --- a/fix +++ b/fix @@ -1 +1 @@ -Subproject commit 13ad7b5abb20567aacc52b8a5d2d487628e5b1c5 +Subproject commit dc537d7106c940c62df0588834119eae2d80525c diff --git a/src/gsi/setupw.f90 b/src/gsi/setupw.f90 index 11a5727ce6..ba798ffb4f 100644 --- a/src/gsi/setupw.f90 +++ b/src/gsi/setupw.f90 @@ -1244,7 +1244,7 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav call vqc_setup(vals,ratio_errors,error,cvar,& cg_t,ibb,ikk,var_jb,rat_err2v,wgt,valqcv) rwgt = rwgt+0.5_r_kind*wgt/wgtlim - valqc=half*(valqcu+valqcv) + valqc=valqcu+valqcv ! Accumulate statistics for obs belonging to this task if (muse(i)) then From 333ec2b43f2ccdc61bf1551faec962689357cc7b Mon Sep 17 00:00:00 2001 From: "john.derber" Date: Thu, 13 Apr 2023 15:18:26 +0000 Subject: [PATCH 4/5] Changes to write_incr for ticket # 558 --- src/gsi/write_incr.f90 | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/src/gsi/write_incr.f90 b/src/gsi/write_incr.f90 index e312199998..ad0b6e3da7 100644 --- a/src/gsi/write_incr.f90 +++ b/src/gsi/write_incr.f90 @@ -76,6 +76,7 @@ subroutine write_fv3_inc_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) use general_sub2grid_mod, only: sub2grid_info use gsi_bundlemod, only: gsi_bundle, gsi_bundlegetpointer + use gsi_bundlemod, only: assignment(=) use control_vectors, only: control_vector use constants, only: one, rad2deg, r1000 @@ -158,7 +159,6 @@ subroutine write_fv3_inc_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) ! set up state space based off of xhatsave ! Convert from control space directly to physical ! space for comparison with obs. - call allocate_preds(sbiasinc) do iii=1,nobs_bins call allocate_state(svalinc(iii)) end do @@ -168,7 +168,10 @@ subroutine write_fv3_inc_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) do iii=1,ntlevs_ens call allocate_state(evalinc(iii)) end do + + call allocate_preds(sbiasinc) call control2state(xhatsave,mvalinc,sbiasinc) + call deallocate_preds(sbiasinc) if (l4dvar) then if (l_hyb_ens) then @@ -193,6 +196,12 @@ subroutine write_fv3_inc_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) end do end if end if + do iii=1,ntlevs_ens + call deallocate_state(evalinc(iii)) + end do + do iii=1,nsubwin + call deallocate_state(mvalinc(iii)) + end do ! Check hydrometeors in control variables iql = getindex(svars3d,'ql') @@ -202,7 +211,7 @@ subroutine write_fv3_inc_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) iqg = getindex(svars3d,'qg') istatus=0 - call gsi_bundlegetpointer(gfs_bundle,'q', sub_qanl, iret); istatus=istatus+iret + call gsi_bundlegetpointer(svalinc(ibin2),'q', sub_qanl, iret); istatus=istatus+iret if (iql>0) call gsi_bundlegetpointer(svalinc(ibin2),'ql', sub_ql, iret); istatus=istatus+iret if (iqi>0) call gsi_bundlegetpointer(svalinc(ibin2),'qi', sub_qi, iret); istatus=istatus+iret if (iqr>0) call gsi_bundlegetpointer(svalinc(ibin2),'qr', sub_qr, iret); istatus=istatus+iret @@ -527,6 +536,10 @@ subroutine write_fv3_inc_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) endif ! ! cleanup and exit call nccheck_incr(nf90_close(ncid_out)) + deallocate(out3d) + do iii=1,nobs_bins + call deallocate_state(svalinc(iii)) + end do if ( mype == mype_out ) then write(6,*) "FV3 netCDF increment written, file= "//trim(filename)//".nc" end if From 768a703843ad75c01f6624c78e14efceec4b59fb Mon Sep 17 00:00:00 2001 From: "john.derber" Date: Fri, 14 Apr 2023 07:14:59 +0000 Subject: [PATCH 5/5] Fix issue found by Ting in review. --- src/gsi/write_incr.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gsi/write_incr.f90 b/src/gsi/write_incr.f90 index ad0b6e3da7..b22be8de30 100644 --- a/src/gsi/write_incr.f90 +++ b/src/gsi/write_incr.f90 @@ -211,7 +211,7 @@ subroutine write_fv3_inc_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) iqg = getindex(svars3d,'qg') istatus=0 - call gsi_bundlegetpointer(svalinc(ibin2),'q', sub_qanl, iret); istatus=istatus+iret + call gsi_bundlegetpointer(gfs_bundle,'q', sub_qanl, iret); istatus=istatus+iret if (iql>0) call gsi_bundlegetpointer(svalinc(ibin2),'ql', sub_ql, iret); istatus=istatus+iret if (iqi>0) call gsi_bundlegetpointer(svalinc(ibin2),'qi', sub_qi, iret); istatus=istatus+iret if (iqr>0) call gsi_bundlegetpointer(svalinc(ibin2),'qr', sub_qr, iret); istatus=istatus+iret