Bugfix for pseudo_RH option (#602)

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,

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
 

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))