diff --git a/src/enkf/gridinfo_gfs.f90 b/src/enkf/gridinfo_gfs.f90 index 317ca2221c..de71153c69 100644 --- a/src/enkf/gridinfo_gfs.f90 +++ b/src/enkf/gridinfo_gfs.f90 @@ -67,6 +67,7 @@ module gridinfo 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' /) ! supported variable names in anavinfo contains diff --git a/src/enkf/gridio_gfs.f90 b/src/enkf/gridio_gfs.f90 index 44432549fe..98779a0e27 100644 --- a/src/enkf/gridio_gfs.f90 +++ b/src/enkf/gridio_gfs.f90 @@ -2026,6 +2026,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,n character(nemsio_charkind) :: field character(len=nf90_max_name) :: time_units logical :: hasfield + character(len=max_varname_length), dimension(n3d) :: no_vars3d real(r_kind) kap,kapr,kap1,clip real(r_single) compress_err @@ -2047,10 +2048,12 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,n call set_ncio_file_flags(vars3d, n3d, vars2d, n2d, write_sfc_file, write_atm_file) - if (write_sfc_file .and. nproc==0 ) then + if (write_sfc_file ) then ! adding the sfc increments requires adjusting several other variables. This is done is a separate ! program. - write(6,*)'gridio/writegriddata: not coded to write sfc analysis, use separate add_incr program instead' + if (nproc == 0) write(6,*)'gridio/writegriddata: not coded to write sfc analysis, will write increment for sfc fields' + no_vars3d='' + call writeincrement(nanal1,nanal2,no_vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_flag) endif nocompress = .true. @@ -3488,7 +3491,7 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, integer :: ql_ind, qi_ind, qr_ind, qs_ind, qg_ind ! netcdf things - integer(i_kind) :: dimids3(3), ncstart(3), nccount(3) + integer(i_kind) :: dimids3(3), ncstart(3), nccount(3), dimids2(2) integer(i_kind) :: ncid_out, lon_dimid, lat_dimid, lev_dimid, ilev_dimid integer(i_kind) :: lonvarid, latvarid, levvarid, pfullvarid, ilevvarid, & hyaivarid, hybivarid, uvarid, vvarid, delpvarid, delzvarid, & @@ -3519,7 +3522,6 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, call set_ncio_file_flags(vars3d, n3d, vars2d, n2d, write_sfc_file, write_atm_file) - if ( write_atm_file) then use_full_hydro = .false. clip = tiny_r_kind read(datestring,*) iadateout @@ -3527,6 +3529,7 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, ncstart = (/1, 1, 1/) nccount = (/nlons, nlats, nlevs/) + if ( write_atm_file) then ne = 0 ensmemloop: do nanal=nanal1,nanal2 ne = ne + 1 @@ -3882,20 +3885,21 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, ! create dimensions based on analysis resolution, not guess call nccheck_incr(nf90_def_dim(ncid_out, "longitude", nlons, lon_dimid)) call nccheck_incr(nf90_def_dim(ncid_out, "latitude", nlats, lat_dimid)) + dimids2 = (/ lon_dimid, lat_dimid /) ! create variables call nccheck_incr(nf90_def_var(ncid_out, "longitude", nf90_real, (/lon_dimid/), lonvarid)) call nccheck_incr(nf90_def_var(ncid_out, "latitude", nf90_real, (/lat_dimid/), latvarid)) - call nccheck_incr(nf90_def_var(ncid_out, "tmp2m_inc", nf90_real, dimids3(1:2), tmp2mvarid)) - call nccheck_incr(nf90_def_var(ncid_out, "spfh2m_inc", nf90_real, dimids3(1:2), spfh2mvarid)) - call nccheck_incr(nf90_def_var(ncid_out, "soilt1_inc", nf90_real, dimids3(1:2), soilt1varid)) - call nccheck_incr(nf90_def_var(ncid_out, "soilt2_inc", nf90_real, dimids3(1:2), soilt2varid)) - call nccheck_incr(nf90_def_var(ncid_out, "soilt3_inc", nf90_real, dimids3(1:2), soilt3varid)) - call nccheck_incr(nf90_def_var(ncid_out, "soilt4_inc", nf90_real, dimids3(1:2), soilt4varid)) - call nccheck_incr(nf90_def_var(ncid_out, "slc1_inc", nf90_real, dimids3(1:2), slc1varid)) - call nccheck_incr(nf90_def_var(ncid_out, "slc2_inc", nf90_real, dimids3(1:2), slc2varid)) - call nccheck_incr(nf90_def_var(ncid_out, "slc3_inc", nf90_real, dimids3(1:2), slc3varid)) - call nccheck_incr(nf90_def_var(ncid_out, "slc4_inc", nf90_real, dimids3(1:2), slc4varid)) - call nccheck_incr(nf90_def_var(ncid_out, "soilsnow_mask", nf90_int, dimids3(1:2), maskvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "tmp2m_inc", nf90_real, dimids2, tmp2mvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "spfh2m_inc", nf90_real, dimids2, spfh2mvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt1_inc", nf90_real, dimids2, soilt1varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt2_inc", nf90_real, dimids2, soilt2varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt3_inc", nf90_real, dimids2, soilt3varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt4_inc", nf90_real, dimids2, soilt4varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc1_inc", nf90_real, dimids2, slc1varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc2_inc", nf90_real, dimids2, slc2varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc3_inc", nf90_real, dimids2, slc3varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc4_inc", nf90_real, dimids2, slc4varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilsnow_mask", nf90_int, dimids2, maskvarid)) ! place global attributes to serial calc_increment output call nccheck_incr(nf90_put_att(ncid_out, nf90_global, "source", "GSI EnKF")) call nccheck_incr(nf90_put_att(ncid_out, nf90_global, "comment", & diff --git a/src/enkf/inflation.f90 b/src/enkf/inflation.f90 index 225967028c..a3a9935f6d 100644 --- a/src/enkf/inflation.f90 +++ b/src/enkf/inflation.f90 @@ -77,7 +77,8 @@ module inflation use constants, only: one, zero, rad2deg, deg2rad use covlocal, only: latval, taper use controlvec, only: ncdim, cvars3d, cvars2d, nc3d, nc2d, clevels -use gridinfo, only: latsgrd, logp, npts, nlevs_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 loadbal, only: indxproc, numptsperproc, npts_max, anal_chunk, anal_chunk_prior use smooth_mod, only: smooth @@ -101,9 +102,10 @@ 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 +integer(i_kind) i,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 +real(r_single), allocatable, dimension(:,:) :: tmp_chunk2,covinfglobal,store_presmooth real(r_single) r fnanalsml = one/(real(nanals-1,r_single)) @@ -231,7 +233,33 @@ subroutine inflate_ens() do nn=1,ncdim call mpi_allreduce(mpi_in_place,covinfglobal(1,nn),npts,mpi_real4,mpi_sum,mpi_comm_world,ierr) enddo + ! do not apply smoothing to soil temp. or soil moisture (not globally defined) + + ind = 0 + do i = 1,8 + this_ind = getindex(cvars2d, vars2d_landonly(i)) + if (this_ind>0) then + ind=ind+1 + soil_index(ind)=this_ind + endif + enddo + + if (ind>0) then + allocate(store_presmooth(npts,ind)) + do i = 1, ind + store_presmooth(:,i) = covinfglobal(:,clevels(nc3d)+soil_index(i)) + enddo + endif + call smooth(covinfglobal) + + if (ind>0) then + do i = 1, ind + covinfglobal(:,clevels(nc3d) + soil_index(i)) = store_presmooth(:,i) + enddo + deallocate(store_presmooth) + endif + where (covinfglobal < covinflatemin) covinfglobal = covinflatemin where (covinfglobal > covinflatemax) covinfglobal = covinflatemax do i=1,numptsperproc(nproc+1)