From accb07e291e0fa707e0ee0adf17c15a83d8e06d8 Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA <33430543+ClaraDraper-NOAA@users.noreply.github.com> Date: Mon, 7 Aug 2023 09:50:18 -0600 Subject: [PATCH] Bugfix for pseudo_RH option (#602) --- regression/regression_namelists.sh | 4 +- src/enkf/controlvec.f90 | 101 +++++++---------------------- src/enkf/params.f90 | 12 +--- 3 files changed, 25 insertions(+), 92 deletions(-) diff --git a/regression/regression_namelists.sh b/regression/regression_namelists.sh index 824c6f071..3668a6f8c 100755 --- a/regression/regression_namelists.sh +++ b/regression/regression_namelists.sh @@ -2155,7 +2155,7 @@ export gsi_namelist=" &nam_enkf datestring=${global_adate},datapath='${DATA}/', analpertwtnh=${analpertwt},analpertwtsh=${analpertwt},analpertwttr=${analpertwt}, - covinflatemax=1.e2,covinflatemin=1,pseudo_rh=.true.,iassim_order=0, + covinflatemax=1.e2,covinflatemin=1,pseudo_rh=.false.,iassim_order=0, corrlengthnh=${corrlength},corrlengthsh=${corrlength},corrlengthtr=${corrlength}, lnsigcutoffnh=${lnsigcutoff},lnsigcutoffsh=${lnsigcutoff},lnsigcutofftr=${lnsigcutoff}, lnsigcutoffpsnh=${lnsigcutoff},lnsigcutoffpssh=${lnsigcutoff},lnsigcutoffpstr=${lnsigcutoff}, @@ -2169,7 +2169,7 @@ export gsi_namelist=" use_gfs_nemsio=${use_gfs_nemsio},use_gfs_ncio=${use_gfs_ncio},imp_physics=$imp_physics,lupp=$lupp, univaroz=.false.,adp_anglebc=.true.,angord=4,use_edges=.false.,emiss_bc=.true., letkf_flag=${letkf_flag},nobsl_max=${nobsl_max},denkf=${denkf},getkf=${getkf}., - nhr_anal=${IAUFHRS_ENKF},nhr_state=${IAUFHRS_ENKF},use_qsatensmean=.true., + nhr_anal=${IAUFHRS_ENKF},nhr_state=${IAUFHRS_ENKF}, lobsdiag_forenkf=$lobsdiag_forenkf, write_spread_diag=$write_spread_diag, modelspace_vloc=$modelspace_vloc, diff --git a/src/enkf/controlvec.f90 b/src/enkf/controlvec.f90 index 4aa2613c6..40dff4dba 100644 --- a/src/enkf/controlvec.f90 +++ b/src/enkf/controlvec.f90 @@ -51,7 +51,7 @@ module controlvec use gridinfo, only: getgridinfo, gridinfo_cleanup, & npts, vars3d_supported, vars2d_supported use params, only: nlevs, nbackgrounds, fgfileprefixes, reducedgrid, & - nanals, pseudo_rh, use_qsatensmean, nlons, nlats,& + nanals, pseudo_rh, nlons, nlats,& nanals_per_iotask, ntasks_io, nanal1, nanal2, & fgsfcfileprefixes, paranc, write_fv3_incr, write_ensmean use kinds, only: r_kind, i_kind, r_double, r_single @@ -64,7 +64,6 @@ module controlvec public :: read_control, write_control, controlvec_cleanup, init_controlvec real(r_single), public, allocatable, dimension(:,:,:,:) :: grdin real(r_double), public, allocatable, dimension(:,:,:,:) :: qsat -real(r_double), public, allocatable, dimension(:,:,:) :: qsatmean integer(i_kind), public :: nc2d, nc3d, ncdim character(len=max_varname_length), allocatable, dimension(:), public :: cvars3d @@ -160,7 +159,7 @@ subroutine init_controlvec() do i = 1, nc2d if (getindex(vars2d_supported, cvars2d(i))<0) then if (nproc .eq. 0) then - print *,'Error: 2D variable ', cvars2d(i), ' is not supported in current version.' + print *,'Error: control 2D variable ', cvars2d(i), ' is not supported in current version.' print *,'Supported variables: ', vars2d_supported endif call stop2(502) @@ -169,7 +168,7 @@ subroutine init_controlvec() do i = 1, nc3d if (getindex(vars3d_supported, cvars3d(i))<0) then if (nproc .eq. 0) then - print *,'Error: 3D variable ', cvars3d(i), ' is not supported in current version.' + print *,'Error: control 3D variable ', cvars3d(i), ' is not supported in current version.' print *,'Supported variables: ', vars3d_supported endif call stop2(502) @@ -192,7 +191,6 @@ subroutine read_control() ! read ensemble members on IO tasks implicit none real(r_double) :: t1,t2 -real(r_double), allocatable, dimension(:) :: qsat_tmp integer(i_kind) :: nb,nlev,ne integer(i_kind) :: q_ind integer(i_kind) :: ierr @@ -229,56 +227,18 @@ subroutine read_control() end if !print *,'min/max qsat',nanal,'=',minval(qsat),maxval(qsat) q_ind = getindex(cvars3d, 'q') - if (use_qsatensmean .and. q_ind>0 ) then - allocate(qsatmean(npts,nlevs,nbackgrounds)) - allocate(qsat_tmp(npts)) - ! compute ensemble mean qsat - qsatmean = 0_r_double - do ne=1,nanals_per_iotask - do nb=1,nbackgrounds - do nlev=1,nlevs - call mpi_allreduce(qsat(:,nlev,nb,ne),qsat_tmp,npts,mpi_real8,mpi_sum,mpi_comm_io,ierr) - qsatmean(:,nlev,nb) = qsatmean(:,nlev,nb) + qsat_tmp - enddo - enddo - enddo - deallocate(qsat_tmp) - qsatmean = qsatmean/real(nanals) - !print *,'min/max qsat ensmean',nanal,'=',minval(qsat),maxval(qsat) - endif if (nproc == 0) then t2 = mpi_wtime() print *,'time in readgridata on root',t2-t1,'secs' end if - !do ne=1,nanals_per_iotask - ! nanal = ne + (nproc-1)*nanals_per_iotask - ! print *,'min/max ps ens mem',nanal,'=',& - ! minval(grdin(:,ncdim,nbackgrounds/2+1,ne)),maxval(grdin(:,ncdim,nbackgrounds/2+1,ne)) - ! print *,'min/max qsat',nanal,'=',& - ! minval(qsat(:,:,nbackgrounds/2+1,ne)),maxval(qsat(:,:,nbackgrounds/2+1,ne)) - !enddo - !if (use_qsatensmean) then - ! print *,'min/max qsatmean proc',nproc,'=',& - ! minval(qsatmean(:,:,nbackgrounds/2+1)),maxval(qsatmean(:,:,nbackgrounds/2+1)) - !endif if (pseudo_rh .and. q_ind > 0) then - if (use_qsatensmean) then - do ne=1,nanals_per_iotask - do nb=1,nbackgrounds - ! create normalized humidity analysis variable. - grdin(:,(q_ind-1)*nlevs+1:q_ind*nlevs,nb,ne) = & - grdin(:,(q_ind-1)*nlevs+1:q_ind*nlevs,nb,ne)/qsatmean(:,:,nb) - enddo - enddo - else - do ne=1,nanals_per_iotask - do nb=1,nbackgrounds - ! create normalized humidity analysis variable. - grdin(:,(q_ind-1)*nlevs+1:q_ind*nlevs,nb,ne) = & - grdin(:,(q_ind-1)*nlevs+1:q_ind*nlevs,nb,ne)/qsat(:,:,nb,ne) - enddo - enddo - endif + do ne=1,nanals_per_iotask + do nb=1,nbackgrounds + ! create normalized humidity analysis variable. + grdin(:,(q_ind-1)*nlevs+1:q_ind*nlevs,nb,ne) = & + grdin(:,(q_ind-1)*nlevs+1:q_ind*nlevs,nb,ne)/qsat(:,:,nb,ne) + enddo + enddo end if endif @@ -299,6 +259,18 @@ subroutine write_control(no_inflate_flag) if (nproc <= ntasks_io-1) then + ! scale q by ensemble qsat, prior to averaging + q_ind = getindex(cvars3d, 'q') + if (pseudo_rh .and. q_ind > 0) then + do ne=1,nanals_per_iotask + do nb=1,nbackgrounds + grdin(:,(q_ind-1)*nlevs+1:q_ind*nlevs,nb,ne) = & + grdin(:,(q_ind-1)*nlevs+1:q_ind*nlevs,nb,ne)*qsat(:,:,nb,ne) + enddo + enddo + endif + + allocate(grdin_mean_tmp(npts,ncdim)) if (nproc == 0) then allocate(grdin_mean(npts,ncdim,nbackgrounds,1)) @@ -345,34 +317,6 @@ subroutine write_control(no_inflate_flag) 100 format('ens. mean anal. increment min/max ',a,2x,g19.12,2x,g19.12) deallocate(grdin_mean_tmp) - q_ind = getindex(cvars3d, 'q') - if (pseudo_rh .and. q_ind > 0) then - if (use_qsatensmean) then - do ne=1,nanals_per_iotask - do nb=1,nbackgrounds - ! re-scale normalized spfh with sat. sphf of ensmean first guess - grdin(:,(q_ind-1)*nlevs+1:q_ind*nlevs,nb,ne) = & - grdin(:,(q_ind-1)*nlevs+1:q_ind*nlevs,nb,ne)*qsatmean(:,:,nb) - enddo - enddo - else - do ne=1,nanals_per_iotask - do nb=1,nbackgrounds - ! re-scale normalized spfh with sat. sphf of first guess - grdin(:,(q_ind-1)*nlevs+1:q_ind*nlevs,nb,ne) = & - grdin(:,(q_ind-1)*nlevs+1:q_ind*nlevs,nb,ne)*qsat(:,:,nb,ne) - enddo - enddo - endif - if (nproc == 0 .and. write_ensmean) then - ! write_ensmean implies use_qsatensmean - do nb=1,nbackgrounds - ! re-scale normalized spfh with sat. sphf of ensmean first guess - grdin_mean(:,(q_ind-1)*nlevs+1:q_ind*nlevs,nb,1) = & - grdin_mean(:,(q_ind-1)*nlevs+1:q_ind*nlevs,nb,1)*qsatmean(:,:,nb) - enddo - endif - end if if (.not. paranc) then if (write_fv3_incr) then call writeincrement(nanal1(nproc),nanal2(nproc),cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin,no_inflate_flag) @@ -427,7 +371,6 @@ subroutine controlvec_cleanup() if (allocated(index_pres)) deallocate(index_pres) if (allocated(grdin)) deallocate(grdin) if (allocated(qsat)) deallocate(qsat) -if (allocated(qsatmean)) deallocate(qsatmean) call gridinfo_cleanup() end subroutine controlvec_cleanup diff --git a/src/enkf/params.f90 b/src/enkf/params.f90 index b21d88abd..62701c24a 100644 --- a/src/enkf/params.f90 +++ b/src/enkf/params.f90 @@ -226,12 +226,6 @@ module params ! EFSOI calculation applications logical,public :: efsoi_flag = .false. -! if true, use ensemble mean qsat in definition of -! normalized humidity analysis variable (instead of -! qsat for each member, which is the default behavior -! when pseudo_rh=.true. If pseudo_rh=.false, use_qsatensmean -! is ignored. -logical,public :: use_qsatensmean = .false. logical,public :: write_spread_diag = .false. ! if true, use jacobian from GSI stored in diag file to compute ! ensemble perturbations in observation space. @@ -261,7 +255,7 @@ module params namelist /nam_enkf/datestring,datapath,iassim_order,nvars,& covinflatemax,covinflatemin,deterministic,sortinc,& mincorrlength_fact,corrlengthnh,corrlengthtr,corrlengthsh,& - varqc,huber,nlons,nlats,smoothparm,use_qsatensmean,& + varqc,huber,nlons,nlats,smoothparm,& readin_localization, zhuberleft,zhuberright,& obtimelnh,obtimeltr,obtimelsh,reducedgrid,& lnsigcutoffnh,lnsigcutofftr,lnsigcutoffsh,& @@ -681,10 +675,6 @@ subroutine read_namelist() letkf_flag) then print *,'warning: no time localization in LETKF!' endif - if ((write_ensmean .and. pseudo_rh) .and. .not. use_qsatensmean) then - print *,'write_ensmean=T requires use_qsatensmean=T when pseudo_rh=T' - call stop2(19) - endif print *, trim(adjustl(datapath))