From 0a6c84abec83511cfbea660a4531eff9a824b8e4 Mon Sep 17 00:00:00 2001
From: Jeffrey S Whitaker <jwhitake@hercules-login-2.hpc.msstate.edu>
Date: Sun, 24 Mar 2024 18:45:39 -0500
Subject: [PATCH 1/9] add ability to taper analysis perts near top of model

---
 src/enkf/gridinfo_gfs.f90 | 32 ++++++++++++++++++++++++++------
 src/enkf/inflation.f90    | 15 +++++++++++----
 src/enkf/params.f90       |  7 ++++++-
 3 files changed, 43 insertions(+), 11 deletions(-)

diff --git a/src/enkf/gridinfo_gfs.f90 b/src/enkf/gridinfo_gfs.f90
index de71153c69..bb1bd994d8 100644
--- a/src/enkf/gridinfo_gfs.f90
+++ b/src/enkf/gridinfo_gfs.f90
@@ -45,7 +45,8 @@ module gridinfo
 
 use mpisetup, only: nproc, mpi_integer, mpi_real4
 use mpimod, only: mpi_comm_world
-use params, only: datapath,nlevs,nlons,nlats,use_gfs_nemsio,use_gfs_ncio,fgfileprefixes
+use params, only: datapath,nlevs,nlons,nlats,use_gfs_nemsio,use_gfs_ncio,fgfileprefixes,&
+                  taperanalperts,taperanalperts_aktop,taperanalperts_akbot
 use kinds, only: r_kind, i_kind, r_double, r_single
 use constants, only: one,zero,pi,cp,rd,grav,rearth,max_varname_length
 use specmod, only: sptezv_s, sptez_s, init_spec_vars, isinitialized, asin_gaulats, &
@@ -57,14 +58,14 @@ module gridinfo
 public :: getgridinfo, gridinfo_cleanup
 integer(i_kind),public :: nlevs_pres, idvc
 real(r_single),public :: ptop
-real(r_single),public, allocatable, dimension(:) :: lonsgrd, latsgrd
+real(r_single),public, allocatable, dimension(:) :: lonsgrd, latsgrd, taper_vert
 ! arrays passed to kdtree2 routines must be single
 real(r_single),public, allocatable, dimension(:,:) :: gridloc
 real(r_single),public, allocatable, dimension(:,:) :: logp
 integer,public :: npts
 integer,public :: ntrunc
 ! supported variable names in anavinfo
-character(len=max_varname_length),public, dimension(13) :: vars3d_supported = (/'u   ', 'v   ', 'tv  ', 'q   ', 'oz  ', 'cw  ', 'tsen', 'prse', &
+character(len=max_varname_length),public, dimension(15) :: vars3d_supported = (/'u   ', 'v   ', 'tv  ', 'dprs', 'delz', 'q   ', 'oz  ', 'cw  ', 'tsen', 'prse', &
                                                                                 'ql  ', 'qi  ', 'qr  ', 'qs  ', 'qg  '/) 
 character(len=max_varname_length),public, dimension(13)  :: vars2d_supported = (/'ps ', 'pst', 'sst', 't2m', 'q2m', 'st1', 'st2', 'st3', 'st4', 'sl1', 'sl2', 'sl3', 'sl4' /)
 character(len=max_varname_length),public, dimension(8)  :: vars2d_landonly = (/'st1', 'st2', 'st3', 'st4', 'sl1', 'sl2', 'sl3', 'sl4' /)
@@ -173,6 +174,8 @@ subroutine getgridinfo(fileprefix, reducedgrid)
    allocate(presslmn(nlons*nlats,nlevs))
    allocate(pressimn(nlons*nlats,nlevs+1))
    allocate(spressmn(nlons*nlats))
+   allocate(taper_vert(nlevs))
+   taper_vert=one
    if (use_gfs_nemsio) then
       call nemsio_readrecv(gfile,'pres','sfc',1,nems_wrk,iret=iret)
       if (iret/=0) then
@@ -221,7 +224,6 @@ subroutine getgridinfo(fileprefix, reducedgrid)
       enddo
       call nemsio_close(gfile, iret=iret)
       ptop = ak(nlevs+1)
-      deallocate(ak,bk)
    else if (use_gfs_ncio) then
       call read_vardata(dset, 'pressfc', values_2d,errcode=iret)
       if (iret /= 0) then
@@ -238,7 +240,7 @@ subroutine getgridinfo(fileprefix, reducedgrid)
          pressimn(:,k) = 0.01_r_kind*ak(nlevs-k+2)+bk(nlevs-k+2)*spressmn(:)
       enddo
       ptop = 0.01_r_kind*ak(1)
-      deallocate(ak,bk,values_2d)
+      deallocate(values_2d)
    else
 ! get pressure from ensemble mean,
 ! distribute to all processors.
@@ -278,7 +280,6 @@ subroutine getgridinfo(fileprefix, reducedgrid)
       enddo
       call sigio_axdata(sigdata,iret)
       ptop = ak(nlevs+1)
-      deallocate(ak,bk)
    endif
    if (reducedgrid) then
       call reducedgrid_init(nlons,nlats,asin_gaulats)
@@ -334,11 +335,28 @@ subroutine getgridinfo(fileprefix, reducedgrid)
       logp(:,nlevs_pres) = -log(spressmn(:))
    endif
    deallocate(spressmn,presslmn,pressimn)
+  ! vertical taper function for ens perts
+  if (taperanalperts) then
+     do k=1,nlevs
+        if (k < nlevs/2 .and. (ak(k) <= taperanalperts_akbot .and. ak(k) >= taperanalperts_aktop)) then
+           taper_vert(nlevs-k+1)= log(ak(k) - taperanalperts_aktop)/log(taperanalperts_akbot - taperanalperts_aktop)
+        else if (bk(k) .eq. 0. .and. ak(k) < taperanalperts_aktop) then
+           taper_vert(nlevs-k+1) = 0.
+        endif
+     enddo
+     print *,'vertical taper for anal perts:'
+     do k=1,nlevs
+        print *,k,ak(nlevs-k+1),bk(nlevs-k+1),taper_vert(k)
+     enddo
+  endif
+  if (allocated(ak)) deallocate(ak)
+  if (allocated(bk)) deallocate(bk)
 end if
 call mpi_bcast(npts,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
 if (nproc .ne. 0) then
    ! allocate arrays on other (non-root) tasks
    allocate(latsgrd(npts),lonsgrd(npts))
+   allocate(taper_vert(nlevs))
    allocate(logp(npts,nlevs_pres)) ! log(ens mean first guess press) on mid-layers
    allocate(gridloc(3,npts))
    ! initialize reducedgrid_mod on other tasks.
@@ -352,6 +370,7 @@ subroutine getgridinfo(fileprefix, reducedgrid)
 enddo
 call mpi_bcast(lonsgrd,npts,mpi_real4,0,MPI_COMM_WORLD,ierr)
 call mpi_bcast(latsgrd,npts,mpi_real4,0,MPI_COMM_WORLD,ierr)
+call mpi_bcast(taper_vert,nlevs,mpi_real4,0,MPI_COMM_WORLD,ierr)
 call mpi_bcast(ptop,1,mpi_real4,0,MPI_COMM_WORLD,ierr)
 !==> precompute cartesian coords of analysis grid points.
 do nn=1,npts
@@ -365,6 +384,7 @@ end subroutine getgridinfo
 subroutine gridinfo_cleanup()
 if (allocated(lonsgrd)) deallocate(lonsgrd)
 if (allocated(latsgrd)) deallocate(latsgrd)
+if (allocated(taper_vert)) deallocate(taper_vert)
 if (allocated(logp)) deallocate(logp)
 if (allocated(gridloc)) deallocate(gridloc)
 end subroutine gridinfo_cleanup
diff --git a/src/enkf/inflation.f90 b/src/enkf/inflation.f90
index 51d1dd1106..baab12427a 100644
--- a/src/enkf/inflation.f90
+++ b/src/enkf/inflation.f90
@@ -71,14 +71,14 @@ module inflation
                   analpertwtnh_rtpp,analpertwtsh_rtpp,analpertwttr_rtpp,&
                   latbound, delat, datapath, covinflatemax, save_inflation, &
                   covinflatemin, nlons, nlats, smoothparm, nbackgrounds,&
-                  covinflatenh,covinflatesh,covinflatetr,lnsigcovinfcutoff
+                  covinflatenh,covinflatesh,covinflatetr,lnsigcovinfcutoff,taperanalperts
 use kinds, only: r_single, i_kind
 use mpeu_util, only: getindex
 use constants, only: one, zero, rad2deg, deg2rad
 use covlocal, only: latval, taper
-use controlvec, only: ncdim, cvars3d, cvars2d, nc3d, nc2d, clevels
+use controlvec, only: ncdim, cvars3d, cvars2d, nc3d, nc2d, clevels, index_pres
 ! note: vars2d_landonly currently only defined for gridio_gfs, but smoothing only coded for gfs.
-use gridinfo, only: latsgrd, logp, npts, nlevs_pres, vars2d_landonly
+use gridinfo, only: latsgrd, logp, npts, nlevs_pres, vars2d_landonly, taper_vert
 use loadbal, only: indxproc, numptsperproc, npts_max, anal_chunk, anal_chunk_prior
 use smooth_mod, only: smooth
 
@@ -102,7 +102,7 @@ subroutine inflate_ens()
 real(r_single),dimension(ndiag) :: sumcoslat,suma,suma2,sumi,sumf,sumitot,sumatot, &
      sumcoslattot,suma2tot,sumftot
 real(r_single) fnanalsml,coslat
-integer(i_kind) i,nn,iunit,ierr,nb,nnlvl,ps_ind, this_ind, ind
+integer(i_kind) i,k,nlev,nn,iunit,ierr,nb,nnlvl,ps_ind, this_ind, ind
 integer(i_kind), dimension(8) :: soil_index
 character(len=500) filename
 real(r_single), allocatable, dimension(:,:) :: tmp_chunk2,covinfglobal,store_presmooth
@@ -302,11 +302,18 @@ subroutine inflate_ens()
 
 ! apply inflation.
 do nn=1,ncdim
+ nlev = index_pres(nn) ! vertical index for i'th control variable
+ if (nlev .eq. nlevs+1) nlev=-1 ! 2d field
  do i=1,numptsperproc(nproc+1)
 
    ! inflate posterior perturbations.
    anal_chunk(:,i,nn,nb) = tmp_chunk2(i,nn)*anal_chunk(:,i,nn,nb)
 
+   ! optionally 'deflate' perturbations to reduce spread near top of model
+   if (taperanalperts .and. nlev > 0) then
+      anal_chunk(:,i,nn,nb) = taper_vert(nlev)*anal_chunk(:,i,nn,nb)
+   endif
+
    ! area mean surface pressure posterior spread, inflation.
    ! (this diagnostic only makes sense for grids that are regular in longitude)
    if (ps_ind > 0 .and. nn == clevels(nc3d) + ps_ind) then 
diff --git a/src/enkf/params.f90 b/src/enkf/params.f90
index 36b0c9c207..6f76530c49 100644
--- a/src/enkf/params.f90
+++ b/src/enkf/params.f90
@@ -254,6 +254,11 @@ module params
 ! write ensemble mean analysis (or analysis increment)
 logical,public :: write_ensmean = .false.
 
+! taper analysis ens perturbations at top of model (gfs only)
+logical, public :: taperanalperts = .false.
+real, public :: taperanalperts_akbot = 500.
+real, public :: taperanalperts_aktop = -1.
+
 namelist /nam_enkf/datestring,datapath,iassim_order,nvars,&
                    covinflatemax,covinflatemin,deterministic,sortinc,&
                    mincorrlength_fact,corrlengthnh,corrlengthtr,corrlengthsh,&
@@ -286,7 +291,7 @@ module params
                    fv3_native, paranc, nccompress, write_fv3_incr,incvars_to_zero,write_ensmean, &
                    corrlengthrdrnh,corrlengthrdrsh,corrlengthrdrtr,&
                    lnsigcutoffrdrnh,lnsigcutoffrdrsh,lnsigcutoffrdrtr,&
-                   l_use_enkf_directZDA
+                   l_use_enkf_directZDA,taperanalperts,taperanalperts_akbot,taperanalperts_aktop
 namelist /nam_wrf/arw,nmm,nmm_restart
 namelist /nam_fv3/fv3fixpath,nx_res,ny_res,ntiles,l_pres_add_saved,l_fv3reg_filecombined, &
                   fv3_io_layout_nx,fv3_io_layout_ny

From 1bfb7567999d4a5c1ac61f74bbc89f8eefdab91a Mon Sep 17 00:00:00 2001
From: Jeffrey S Whitaker <jwhitake@hercules-login-2.hpc.msstate.edu>
Date: Sun, 24 Mar 2024 18:55:40 -0500
Subject: [PATCH 2/9] update

---
 src/enkf/gridinfo_gfs.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/enkf/gridinfo_gfs.f90 b/src/enkf/gridinfo_gfs.f90
index bb1bd994d8..1ca6a43af0 100644
--- a/src/enkf/gridinfo_gfs.f90
+++ b/src/enkf/gridinfo_gfs.f90
@@ -65,7 +65,7 @@ module gridinfo
 integer,public :: npts
 integer,public :: ntrunc
 ! supported variable names in anavinfo
-character(len=max_varname_length),public, dimension(15) :: vars3d_supported = (/'u   ', 'v   ', 'tv  ', 'dprs', 'delz', 'q   ', 'oz  ', 'cw  ', 'tsen', 'prse', &
+character(len=max_varname_length),public, dimension(13) :: vars3d_supported = (/'u   ', 'v   ', 'tv  ', 'q   ', 'oz  ', 'cw  ', 'tsen', 'prse', &
                                                                                 'ql  ', 'qi  ', 'qr  ', 'qs  ', 'qg  '/) 
 character(len=max_varname_length),public, dimension(13)  :: vars2d_supported = (/'ps ', 'pst', 'sst', 't2m', 'q2m', 'st1', 'st2', 'st3', 'st4', 'sl1', 'sl2', 'sl3', 'sl4' /)
 character(len=max_varname_length),public, dimension(8)  :: vars2d_landonly = (/'st1', 'st2', 'st3', 'st4', 'sl1', 'sl2', 'sl3', 'sl4' /)

From 5c6d350cdf0ce4c145609d32c6a560bff7365290 Mon Sep 17 00:00:00 2001
From: Jeffrey S Whitaker <jwhitake@hercules-login-2.hpc.msstate.edu>
Date: Sun, 24 Mar 2024 19:59:51 -0500
Subject: [PATCH 3/9] add taper_vert to regional models (but set to 1.0)

---
 src/enkf/gridinfo_fv3reg.f90 | 8 ++++++--
 src/enkf/gridinfo_nmmb.f90   | 6 ++++++
 src/enkf/gridinfo_wrf.f90    | 4 ++++
 3 files changed, 16 insertions(+), 2 deletions(-)

diff --git a/src/enkf/gridinfo_fv3reg.f90 b/src/enkf/gridinfo_fv3reg.f90
index 9a16f0ca03..ef81e4a9ac 100644
--- a/src/enkf/gridinfo_fv3reg.f90
+++ b/src/enkf/gridinfo_fv3reg.f90
@@ -65,6 +65,7 @@ module gridinfo
 public :: ak,bk,eta1_ll,eta2_ll
 real(r_single),public :: ptop
 real(r_single),public, allocatable, dimension(:) :: lonsgrd, latsgrd
+real(r_single),public, allocatalbe, dimension(:) :: taper_vert
 ! arrays passed to kdtree2 routines must be single
 real(r_single),public, allocatable, dimension(:,:) :: gridloc
 real(r_single),public, allocatable, dimension(:,:) :: logp
@@ -283,7 +284,8 @@ subroutine getgridinfo(fileprefix, reducedgrid)
 if(nproc == 0)  then
   nn = 0
   npts=nx_res*ny_res
-  allocate(latsgrd(npts),lonsgrd(npts))
+  allocate(latsgrd(npts),lonsgrd(npts),taper_vert(nlevs))
+  taper_vert=one
 endif
 if(any (mpi_id_group == nproc))  then ! if paranc=.fales., this equal to "nproc== 00 
   call mpi_comm_rank(mpi_comm_userread,iope, ierror)
@@ -463,7 +465,7 @@ subroutine getgridinfo(fileprefix, reducedgrid)
 allocate(gridloc(3,npts))
 if (nproc .ne. 0) then
    ! allocate arrays on other (non-root) tasks
-   allocate(latsgrd(npts),lonsgrd(npts))
+   allocate(latsgrd(npts),lonsgrd(npts),taper_vert(nlevs))
    allocate(logp(npts,nlevs_pres)) ! log(ens mean first guess press) on mid-layers
    allocate(eta1_ll(nlevsp1),eta2_ll(nlevsp1))
 endif
@@ -473,6 +475,7 @@ subroutine getgridinfo(fileprefix, reducedgrid)
 enddo
 call mpi_bcast(lonsgrd,npts,mpi_real4,0,MPI_COMM_WORLD,ierr)
 call mpi_bcast(latsgrd,npts,mpi_real4,0,MPI_COMM_WORLD,ierr)
+call mpi_bcast(taper_vert,nlevs,mpi_real4,0,MPI_COMM_WORLD,ierr)
 call mpi_bcast(eta1_ll,nlevsp1,mpi_real4,0,MPI_COMM_WORLD,ierr)
 call mpi_bcast(eta2_ll,nlevsp1,mpi_real4,0,MPI_COMM_WORLD,ierr)
 call mpi_bcast(ptop,1,mpi_real4,0,MPI_COMM_WORLD,ierr)
@@ -489,6 +492,7 @@ end subroutine getgridinfo
 subroutine gridinfo_cleanup()
 if (allocated(lonsgrd)) deallocate(lonsgrd)
 if (allocated(latsgrd)) deallocate(latsgrd)
+if (allocated(taper_vert)) deallocate(taper_vert)
 if (allocated(logp)) deallocate(logp)
 if (allocated(gridloc)) deallocate(gridloc)
 end subroutine gridinfo_cleanup
diff --git a/src/enkf/gridinfo_nmmb.f90 b/src/enkf/gridinfo_nmmb.f90
index df6c22ee1e..05a4b5a932 100644
--- a/src/enkf/gridinfo_nmmb.f90
+++ b/src/enkf/gridinfo_nmmb.f90
@@ -16,6 +16,7 @@ module gridinfo
 integer(i_kind),public :: nlevs_pres
 real(r_single),public :: ptop
 real(r_single),public, allocatable, dimension(:) :: lonsgrd, latsgrd
+real(r_single),public, allocatalbe, dimension(:) :: taper_vert
 ! arrays passed to kdtree2 routines must be single
 real(r_single),public, allocatable, dimension(:,:) :: gridloc
 real(r_single),public, allocatable, dimension(:,:) :: logp
@@ -124,6 +125,8 @@ subroutine getgridinfo(fileprefix, reducedgrid)
    allocate(latsgrd(npts),lonsgrd(npts))
    allocate(logp(npts,nlevs_pres)) ! log(ens mean first guess press) on mid-layers
    allocate(gridloc(3,npts))
+   allocate(taper_vert(nlevs))
+   taper_vert=one
    lonsgrd = lons; latsgrd = lats
    print *,'min/max lonsgrd',minval(lonsgrd),maxval(lonsgrd)
    print *,'min/max latsgrd',minval(latsgrd),maxval(latsgrd)
@@ -165,6 +168,7 @@ subroutine getgridinfo(fileprefix, reducedgrid)
 if (nproc .ne. 0) then
    ! allocate arrays on other (non-root) tasks
    allocate(latsgrd(npts),lonsgrd(npts))
+   allocate(taper_vert(nlevs))
    allocate(logp(npts,nlevs_pres)) ! log(ens mean first guess press) on mid-layers
    allocate(gridloc(3,npts))
 endif
@@ -174,6 +178,7 @@ subroutine getgridinfo(fileprefix, reducedgrid)
 enddo
 call mpi_bcast(lonsgrd,npts,mpi_real4,0,MPI_COMM_WORLD,ierr)
 call mpi_bcast(latsgrd,npts,mpi_real4,0,MPI_COMM_WORLD,ierr)
+call mpi_bcast(taper_vert,nlevs,mpi_real4,0,MPI_COMM_WORLD,ierr)
 call mpi_bcast(ptop,1,mpi_real4,0,MPI_COMM_WORLD,ierr)
 
 !==> precompute cartesian coords of analysis grid points.
@@ -188,6 +193,7 @@ end subroutine getgridinfo
 subroutine gridinfo_cleanup()
 if (allocated(lonsgrd)) deallocate(lonsgrd)
 if (allocated(latsgrd)) deallocate(latsgrd)
+if (allocated(taper_vert)) deallocate(taper_vert)
 if (allocated(logp)) deallocate(logp)
 if (allocated(gridloc)) deallocate(gridloc)
 end subroutine gridinfo_cleanup
diff --git a/src/enkf/gridinfo_wrf.f90 b/src/enkf/gridinfo_wrf.f90
index e67b827b41..ef734d4a6a 100644
--- a/src/enkf/gridinfo_wrf.f90
+++ b/src/enkf/gridinfo_wrf.f90
@@ -63,6 +63,7 @@ module gridinfo
   real(r_single),      dimension(:,:), allocatable, public     :: gridloc
   real(r_single),      dimension(:),   allocatable, public     :: lonsgrd
   real(r_single),      dimension(:),   allocatable, public     :: latsgrd
+  real(r_single),      dimension(:),   allocatable, public     :: taper_vert
   real(r_single),                                   public     :: ptop
   integer(i_long),                                  public     :: npts
   integer(i_kind),                                  public     :: nlevs_pres
@@ -211,7 +212,9 @@ subroutine getgridinfo_arw(fileprefix)
     ! Allocate memory for global arrays
     if(.not. allocated(lonsgrd)) allocate(lonsgrd(npts))
     if(.not. allocated(latsgrd)) allocate(latsgrd(npts))
+    if(.not. allocated(taper_vert)) allocate(taper_vert(nlevs))
     if(.not. allocated(logp))    allocate(logp(npts,nlevs_pres))
+    taper_vert = one
 
     !======================================================================
     ! Begin: Ingest all grid variables required for EnKF routines and
@@ -848,6 +851,7 @@ end subroutine dot2cross
   subroutine gridinfo_cleanup()
     if (allocated(lonsgrd))       deallocate(lonsgrd)
     if (allocated(latsgrd))       deallocate(latsgrd)
+    if (allocated(taper_vert))    deallocate(taper_vert)
     if (allocated(logp))          deallocate(logp)
     if (allocated(gridloc))       deallocate(gridloc)
   end subroutine gridinfo_cleanup

From c2f8bbf0fbfed0e49a4f56c1be78e7fa13a479d5 Mon Sep 17 00:00:00 2001
From: Jeffrey S Whitaker <jwhitake@hercules-login-2.hpc.msstate.edu>
Date: Sun, 24 Mar 2024 20:25:22 -0500
Subject: [PATCH 4/9] fix typo

---
 src/enkf/gridinfo_fv3reg.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/enkf/gridinfo_fv3reg.f90 b/src/enkf/gridinfo_fv3reg.f90
index ef81e4a9ac..ae671f22f5 100644
--- a/src/enkf/gridinfo_fv3reg.f90
+++ b/src/enkf/gridinfo_fv3reg.f90
@@ -65,7 +65,7 @@ module gridinfo
 public :: ak,bk,eta1_ll,eta2_ll
 real(r_single),public :: ptop
 real(r_single),public, allocatable, dimension(:) :: lonsgrd, latsgrd
-real(r_single),public, allocatalbe, dimension(:) :: taper_vert
+real(r_single),public, allocatable, dimension(:) :: taper_vert
 ! arrays passed to kdtree2 routines must be single
 real(r_single),public, allocatable, dimension(:,:) :: gridloc
 real(r_single),public, allocatable, dimension(:,:) :: logp

From 9f7c049e10814c9e5f1a4d5ccb3ffc124c83c9d1 Mon Sep 17 00:00:00 2001
From: Jeffrey S Whitaker <jwhitake@hercules-login-2.hpc.msstate.edu>
Date: Sun, 24 Mar 2024 20:35:31 -0500
Subject: [PATCH 5/9] fix typo

---
 src/enkf/gridinfo_nmmb.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/enkf/gridinfo_nmmb.f90 b/src/enkf/gridinfo_nmmb.f90
index 05a4b5a932..724430a047 100644
--- a/src/enkf/gridinfo_nmmb.f90
+++ b/src/enkf/gridinfo_nmmb.f90
@@ -16,7 +16,7 @@ module gridinfo
 integer(i_kind),public :: nlevs_pres
 real(r_single),public :: ptop
 real(r_single),public, allocatable, dimension(:) :: lonsgrd, latsgrd
-real(r_single),public, allocatalbe, dimension(:) :: taper_vert
+real(r_single),public, allocatable, dimension(:) :: taper_vert
 ! arrays passed to kdtree2 routines must be single
 real(r_single),public, allocatable, dimension(:,:) :: gridloc
 real(r_single),public, allocatable, dimension(:,:) :: logp

From 432220d40eaaca7fd77117406e43c61089a934ee Mon Sep 17 00:00:00 2001
From: Jeffrey S Whitaker <jwhitake@hercules-login-2.hpc.msstate.edu>
Date: Sun, 24 Mar 2024 20:45:28 -0500
Subject: [PATCH 6/9] fix name clash

---
 src/enkf/gridinfo_nmmb.f90 | 3 ++-
 src/enkf/gridinfo_wrf.f90  | 5 +++--
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/enkf/gridinfo_nmmb.f90 b/src/enkf/gridinfo_nmmb.f90
index 724430a047..33b487354e 100644
--- a/src/enkf/gridinfo_nmmb.f90
+++ b/src/enkf/gridinfo_nmmb.f90
@@ -1,6 +1,7 @@
 module gridinfo
 
-use mpisetup
+use mpisetup, only: nproc, mpi_integer, mpi_real4
+use mpimod, only: mpi_comm_world
 use params, only: datapath,nlevs,datestring,&
                   nmmb,regional,nlons,nlats,nbackgrounds,fgfileprefixes
 use kinds, only: r_kind, i_kind, r_double, r_single
diff --git a/src/enkf/gridinfo_wrf.f90 b/src/enkf/gridinfo_wrf.f90
index ef734d4a6a..4ad80aaa60 100644
--- a/src/enkf/gridinfo_wrf.f90
+++ b/src/enkf/gridinfo_wrf.f90
@@ -32,12 +32,13 @@ module gridinfo
 
   ! Define associated modules
 
-  use constants, only: rearth_equator, omega, pi, deg2rad, zero, rad2deg,    &
+  use constants, only: rearth_equator, omega, pi, deg2rad, zero, one, rad2deg, &
                        rearth,max_varname_length
   use kinds,     only: i_kind, r_kind, r_single, i_long, r_double
   use params,    only: datapath, nlevs, nlons, nlats,           &
                        arw, nmm
-  use mpisetup
+  use mpisetup, only: nproc, mpi_integer, mpi_real4,mpi_status
+  use mpimod, only: mpi_comm_world
   use netcdf_io
 
   implicit none

From a0b482bfc059d97d9ffa23e541c528634c726235 Mon Sep 17 00:00:00 2001
From: Jeffrey S Whitaker <jwhitake@hercules-login-2.hpc.msstate.edu>
Date: Mon, 25 Mar 2024 13:08:43 -0500
Subject: [PATCH 7/9] define tapervert using ak,bk if taperanalperts=T

---
 src/enkf/gridinfo_fv3reg.f90 | 28 ++++++++++++++++++++++------
 1 file changed, 22 insertions(+), 6 deletions(-)

diff --git a/src/enkf/gridinfo_fv3reg.f90 b/src/enkf/gridinfo_fv3reg.f90
index ae671f22f5..3c3a288468 100644
--- a/src/enkf/gridinfo_fv3reg.f90
+++ b/src/enkf/gridinfo_fv3reg.f90
@@ -47,7 +47,8 @@ module gridinfo
 use mpimod, only: mpi_comm_world
 use params, only: datapath,nlevs,nlons,nlats,use_gfs_nemsio, fgfileprefixes, &
                   fv3fixpath, nx_res,ny_res, ntiles,l_fv3reg_filecombined,paranc, &
-                  fv3_io_layout_nx,fv3_io_layout_ny
+                  fv3_io_layout_nx,fv3_io_layout_ny,taperanalperts,taperanalperts_akbot, &
+                  taperanalperts_aktop
           
 use kinds, only: r_kind, i_kind, r_double, r_single
 use constants, only: one,zero,pi,cp,rd,grav,rearth,max_varname_length
@@ -165,12 +166,28 @@ subroutine getgridinfo(fileprefix, reducedgrid)
        eta2_ll(i)=bk(i)
    enddo
 
-
-
-
    ptop = eta1_ll(nlevsp1)
    call nc_check( nf90_close(file_id),&
    myname_,'close '//trim(filename) )
+
+   ! vertical taper function for ens perts
+   allocate(taper_vert(nlevs))
+   if (taperanalperts) then
+      do k=1,nlevs
+         if (k < nlevs/2 .and. (ak(k) <= taperanalperts_akbot .and. ak(k) >= taperanalperts_aktop)) then
+            taper_vert(nlevs-k+1)= log(ak(k) - taperanalperts_aktop)/log(taperanalperts_akbot - taperanalperts_aktop)
+         else if (bk(k) .eq. 0. .and. ak(k) < taperanalperts_aktop) then
+            taper_vert(nlevs-k+1) = 0.
+         endif
+      enddo
+      print *,'vertical taper for anal perts:'
+      do k=1,nlevs
+         print *,k,ak(nlevs-k+1),bk(nlevs-k+1),taper_vert(k)
+      enddo
+   else
+      taper_vert = one
+   endif
+
    deallocate(ak,bk)
 endif ! root task
 
@@ -284,8 +301,7 @@ subroutine getgridinfo(fileprefix, reducedgrid)
 if(nproc == 0)  then
   nn = 0
   npts=nx_res*ny_res
-  allocate(latsgrd(npts),lonsgrd(npts),taper_vert(nlevs))
-  taper_vert=one
+  allocate(latsgrd(npts),lonsgrd(npts))
 endif
 if(any (mpi_id_group == nproc))  then ! if paranc=.fales., this equal to "nproc== 00 
   call mpi_comm_rank(mpi_comm_userread,iope, ierror)

From 3116c80495f8258d8ef119ed16652408c21f7664 Mon Sep 17 00:00:00 2001
From: Jeffrey S Whitaker <jwhitake@hercules-login-3.hpc.msstate.edu>
Date: Thu, 28 Mar 2024 12:08:16 -0500
Subject: [PATCH 8/9] address @RussTreadon-NOAA's comments

---
 src/enkf/gridinfo_fv3reg.f90 | 8 ++++----
 src/enkf/gridinfo_gfs.f90    | 8 ++++----
 src/enkf/inflation.f90       | 6 +++---
 src/enkf/params.f90          | 4 ++--
 4 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/src/enkf/gridinfo_fv3reg.f90 b/src/enkf/gridinfo_fv3reg.f90
index 3c3a288468..337c9ba682 100644
--- a/src/enkf/gridinfo_fv3reg.f90
+++ b/src/enkf/gridinfo_fv3reg.f90
@@ -135,7 +135,7 @@ subroutine getgridinfo(fileprefix, reducedgrid)
 
 !when paranc=.false, fv3_io_layout_nx=fv3_io_layout_ny=1
 ! read data on root task
-if (nproc .eq. 0) then
+if (nproc == 0) then
 
    !  read ak,bk from ensmean fv_core.res.nc
    !  read nx,ny and nz from fv_core.res.nc
@@ -176,8 +176,8 @@ subroutine getgridinfo(fileprefix, reducedgrid)
       do k=1,nlevs
          if (k < nlevs/2 .and. (ak(k) <= taperanalperts_akbot .and. ak(k) >= taperanalperts_aktop)) then
             taper_vert(nlevs-k+1)= log(ak(k) - taperanalperts_aktop)/log(taperanalperts_akbot - taperanalperts_aktop)
-         else if (bk(k) .eq. 0. .and. ak(k) < taperanalperts_aktop) then
-            taper_vert(nlevs-k+1) = 0.
+         else if (bk(k) == zero .and. ak(k) < taperanalperts_aktop) then
+            taper_vert(nlevs-k+1) = zero
          endif
       enddo
       print *,'vertical taper for anal perts:'
@@ -194,7 +194,7 @@ subroutine getgridinfo(fileprefix, reducedgrid)
 allocate(nxlocgroup(fv3_io_layout_nx,fv3_io_layout_ny))
 allocate(nylocgroup(fv3_io_layout_nx,fv3_io_layout_ny))
 
-if(nproc.eq.0) then 
+if(nproc == 0) then 
   ii=0
   do j=1,fv3_io_layout_ny
     do i=1,fv3_io_layout_nx
diff --git a/src/enkf/gridinfo_gfs.f90 b/src/enkf/gridinfo_gfs.f90
index 1ca6a43af0..efbd7a2959 100644
--- a/src/enkf/gridinfo_gfs.f90
+++ b/src/enkf/gridinfo_gfs.f90
@@ -106,7 +106,7 @@ subroutine getgridinfo(fileprefix, reducedgrid)
 kapr = cp/rd
 kap1 = kap + one
 nlevs_pres=nlevs+1
-if (nproc .eq. 0) then
+if (nproc == 0) then
 filename = trim(adjustl(datapath))//trim(adjustl(fileprefix))//"ensmean"
 if (use_gfs_nemsio) then
      call nemsio_init(iret=iret)
@@ -169,7 +169,7 @@ subroutine getgridinfo(fileprefix, reducedgrid)
 ! initialize spectral module on all tasks.
 if (.not. isinitialized) call init_spec_vars(nlons,nlats,ntrunc,4)
 
-if (nproc .eq. 0) then
+if (nproc == 0) then
    ! get pressure, lat/lon information from ensemble mean file.
    allocate(presslmn(nlons*nlats,nlevs))
    allocate(pressimn(nlons*nlats,nlevs+1))
@@ -340,8 +340,8 @@ subroutine getgridinfo(fileprefix, reducedgrid)
      do k=1,nlevs
         if (k < nlevs/2 .and. (ak(k) <= taperanalperts_akbot .and. ak(k) >= taperanalperts_aktop)) then
            taper_vert(nlevs-k+1)= log(ak(k) - taperanalperts_aktop)/log(taperanalperts_akbot - taperanalperts_aktop)
-        else if (bk(k) .eq. 0. .and. ak(k) < taperanalperts_aktop) then
-           taper_vert(nlevs-k+1) = 0.
+        else if (bk(k) == zero .and. ak(k) < taperanalperts_aktop) then
+           taper_vert(nlevs-k+1) = zero
         endif
      enddo
      print *,'vertical taper for anal perts:'
diff --git a/src/enkf/inflation.f90 b/src/enkf/inflation.f90
index baab12427a..c80cc99c10 100644
--- a/src/enkf/inflation.f90
+++ b/src/enkf/inflation.f90
@@ -113,7 +113,7 @@ subroutine inflate_ens()
 if (analpertwtnh_rtpp > 1.e-5_r_single .and. &
     analpertwtnh_rtpp > 1.e-5_r_single .and. &
     analpertwttr_rtpp > 1.e-5_r_single) then
-if (nproc .eq. 0) print *,'performing RTPP inflation...'
+if (nproc == 0) print *,'performing RTPP inflation...'
 nbloop: do nb=1,nbackgrounds ! loop over time levels in background
 ! First perform RTPP ensemble inflation,
 ! as first described in:
@@ -139,7 +139,7 @@ subroutine inflate_ens()
     abs(analpertwttr) < 1.e-5_r_single .and. &
     abs(analpertwtsh) < 1.e-5_r_single) return
 
-if (nproc .eq. 0) print *,'performing RTPS inflation...'
+if (nproc == 0) print *,'performing RTPS inflation...'
 
 ! now perform RTPS inflation
 nbloop2: do nb=1,nbackgrounds ! loop over time levels in background
@@ -303,7 +303,7 @@ subroutine inflate_ens()
 ! apply inflation.
 do nn=1,ncdim
  nlev = index_pres(nn) ! vertical index for i'th control variable
- if (nlev .eq. nlevs+1) nlev=-1 ! 2d field
+ if (nlev == nlevs+1) nlev=-1 ! 2d field
  do i=1,numptsperproc(nproc+1)
 
    ! inflate posterior perturbations.
diff --git a/src/enkf/params.f90 b/src/enkf/params.f90
index 6f76530c49..7f7dc14d6c 100644
--- a/src/enkf/params.f90
+++ b/src/enkf/params.f90
@@ -256,8 +256,8 @@ module params
 
 ! taper analysis ens perturbations at top of model (gfs only)
 logical, public :: taperanalperts = .false.
-real, public :: taperanalperts_akbot = 500.
-real, public :: taperanalperts_aktop = -1.
+real(r_kind), public :: taperanalperts_akbot = 500.
+real(r_kind), public :: taperanalperts_aktop = -1.
 
 namelist /nam_enkf/datestring,datapath,iassim_order,nvars,&
                    covinflatemax,covinflatemin,deterministic,sortinc,&

From 9df1215a4aff6998e8cd93eef9307c686b42f9c3 Mon Sep 17 00:00:00 2001
From: Jeffrey S Whitaker <jwhitake@hercules-login-3.hpc.msstate.edu>
Date: Thu, 28 Mar 2024 12:32:43 -0500
Subject: [PATCH 9/9] specify kind of default namelist parameters

---
 src/enkf/params.f90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/enkf/params.f90 b/src/enkf/params.f90
index 7f7dc14d6c..f2a52d9a1a 100644
--- a/src/enkf/params.f90
+++ b/src/enkf/params.f90
@@ -256,8 +256,8 @@ module params
 
 ! taper analysis ens perturbations at top of model (gfs only)
 logical, public :: taperanalperts = .false.
-real(r_kind), public :: taperanalperts_akbot = 500.
-real(r_kind), public :: taperanalperts_aktop = -1.
+real(r_kind), public :: taperanalperts_akbot = 500.0_r_kind
+real(r_kind), public :: taperanalperts_aktop = -1.0_r_kind
 
 namelist /nam_enkf/datestring,datapath,iassim_order,nvars,&
                    covinflatemax,covinflatemin,deterministic,sortinc,&