add ability to taper analysis perts near top of model (#729)

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
@@ -65,6 +66,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, 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
@@ -133,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
@@ -164,19 +166,35 @@ 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) == zero .and. ak(k) < taperanalperts_aktop) then
+ taper_vert(nlevs-k+1) = zero
+ 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
 
 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
@@ -463,7 +481,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 +491,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 +508,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

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,7 +58,7 @@ 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
@@ -105,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)
@@ -168,11 +169,13 @@ 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))
  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) == zero .and. ak(k) < taperanalperts_aktop) then
+ taper_vert(nlevs-k+1) = zero
+ 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

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
@@ -16,6 +17,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, 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
@@ -124,6 +126,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 +169,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 +179,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 +194,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

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
@@ -63,6 +64,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 +213,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 +852,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

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
@@ -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
@@ -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 == 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 

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(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,&
  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