removed misleading use_qsatensmean option for pseudo_RH, plus

minor bugfixes for running in debug mode.

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
@@ -214,7 +212,8 @@ subroutine read_control()
 allocate(grdin(npts,ncdim,nbackgrounds,nanals_per_iotask))
 ! if only updating the sfc fields, qsat will not be calculated in readgriddata
 ! only allocate if needed.
-allocate(qsat(npts,nlevs,nbackgrounds,nanals_per_iotask))
+q_ind = getindex(cvars3d, 'q')
+if (q_ind > 0)  allocate(qsat(npts,nlevs,nbackgrounds,nanals_per_iotask))
 if (paranc) then
    if (nproc == 0) t1 = mpi_wtime()
    call readgriddata_pnc(cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,nbackgrounds, &
@@ -228,23 +227,6 @@ subroutine read_control()
    end if
    !print *,'min/max qsat',nanal,'=',minval(qsat),maxval(qsat)
    q_ind = getindex(cvars3d, 'q')
-   if (pseudo_RH .and. 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'
@@ -256,34 +238,18 @@ subroutine read_control()
    !   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
 
-deallocate(qsat)
-
 end subroutine read_control
 
 subroutine write_control(no_inflate_flag)
@@ -300,6 +266,17 @@ 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))
@@ -346,34 +323,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)
@@ -428,7 +377,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/enkf_obsmod.f90

@@ -262,7 +262,8 @@ subroutine readobs()
 allocate(corrlengthsq(nobstot),lnsigl(nobstot),obtimel(nobstot))
 lnsigl=1.e10
 do nob=1,nobstot
-   oblnp(nob) = -log(obpress(nob)) ! distance measured in log(p) units
+   ! max to prevent crash got 0 pressure obs.
+   oblnp(nob) = -log(max(obpress(nob),0.01)) ! distance measured in log(p) units
    if (obloclon(nob) < zero) obloclon(nob) = obloclon(nob) + 360._r_single
    radlon=deg2rad*obloclon(nob)
    radlat=deg2rad*obloclat(nob)

src/enkf/gridio_gfs.f90

@@ -782,8 +782,8 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,
 
   !==> get U,V,temp,q,ps on gaussian grid.
   ! u is first nlevs, v is second, t is third, then tracers.
+  clip=tiny_r_kind
   if (use_gfs_nemsio) then
-     clip=tiny_r_kind
      do k=1,nlevs
         call nemsio_readrecv(gfile,'ugrd','mid layer',k,nems_wrk,iret=iret)
         if (iret/=0) then

src/enkf/observer_gfs.f90

@@ -95,8 +95,9 @@ subroutine setup_linhx(rlat, rlon, time, ix, delx, ixp, delxp, iy, dely,  &
   delx = max(zero,min(delx,one))
 
   iyp = 1
-  do while (iyp <= nlons .and. lonsgrd(ix*nlons + iyp) <= rlon)
+  do while (lonsgrd(ix*nlons + iyp) <= rlon)
     iyp = iyp + 1
+    if (iyp > nlons) exit
   enddo
   iy = iyp - 1
   if(iy < 1)     iy = iy + nlons

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))
@@ -747,6 +737,7 @@ subroutine read_namelist()
       statesfcfileprefixes(nstatefields+1)="bfg_"//datestring//"_fhr"//charfhr_state(nstatefields+1)//"_"
    end if
    nstatefields = nstatefields+1
+   if (nstatefields == 7) exit
 end do
 
 do nb=1,nbackgrounds