Merge branch 'NOAA-EMC:develop' into develop

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