diff --git a/src/enkf/gridinfo_fv3reg.f90 b/src/enkf/gridinfo_fv3reg.f90 index 9a16f0ca03..337c9ba682 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 @@ -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 diff --git a/src/enkf/gridinfo_gfs.f90 b/src/enkf/gridinfo_gfs.f90 index de71153c69..efbd7a2959 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,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 diff --git a/src/enkf/gridinfo_nmmb.f90 b/src/enkf/gridinfo_nmmb.f90 index df6c22ee1e..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 @@ -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 diff --git a/src/enkf/gridinfo_wrf.f90 b/src/enkf/gridinfo_wrf.f90 index e67b827b41..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 @@ -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 diff --git a/src/enkf/inflation.f90 b/src/enkf/inflation.f90 index 51d1dd1106..c80cc99c10 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 @@ -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 diff --git a/src/enkf/params.f90 b/src/enkf/params.f90 index 36b0c9c207..f2a52d9a1a 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(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