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 e4631f4e2d..8afb012fcf 100644 --- a/src/enkf/gridio_gfs.f90 +++ b/src/enkf/gridio_gfs.f90 @@ -926,6 +926,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, endif ! use_full_hydro enddo else if (use_gfs_ncio) then + clip=tiny_r_kind call read_vardata(dset, 'ugrd', ug3d,errcode=iret) if (iret /= 0) then print *,'error reading ugrd' @@ -1203,36 +1204,36 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, call copytogrdin(ug,grdin(:,levels(n3d) + soilt4_ind,nb,ne)) endif if (slc1_ind > 0) then - call read_vardata(dset_sfc, 'slc1', values_2d, errcode=iret) + call read_vardata(dset_sfc, 'soill1', values_2d, errcode=iret) if (iret /= 0) then - print *,'error reading slc1' + print *,'error reading soill1' call stop2(22) endif ug = reshape(values_2d,(/nlons*nlats/)) call copytogrdin(ug,grdin(:,levels(n3d) + slc1_ind,nb,ne)) endif if (slc2_ind > 0) then - call read_vardata(dset_sfc, 'slc2', values_2d, errcode=iret) + call read_vardata(dset_sfc, 'soill2', values_2d, errcode=iret) if (iret /= 0) then - print *,'error reading slc2' + print *,'error reading soill2' call stop2(22) endif ug = reshape(values_2d,(/nlons*nlats/)) call copytogrdin(ug,grdin(:,levels(n3d) + slc2_ind,nb,ne)) endif if (slc3_ind > 0) then - call read_vardata(dset_sfc, 'slc3', values_2d, errcode=iret) + call read_vardata(dset_sfc, 'soill3', values_2d, errcode=iret) if (iret /= 0) then - print *,'error reading slc3' + print *,'error reading soill3' call stop2(22) endif ug = reshape(values_2d,(/nlons*nlats/)) call copytogrdin(ug,grdin(:,levels(n3d) + slc3_ind,nb,ne)) endif if (slc4_ind > 0) then - call read_vardata(dset_sfc, 'slc4', values_2d, errcode=iret) + call read_vardata(dset_sfc, 'soill4', values_2d, errcode=iret) if (iret /= 0) then - print *,'error reading slc2' + print *,'error reading soill4' call stop2(22) endif ug = reshape(values_2d,(/nlons*nlats/)) @@ -2122,6 +2123,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 @@ -2143,10 +2145,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. @@ -3584,7 +3588,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, & @@ -3615,7 +3619,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 @@ -3623,6 +3626,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 @@ -3978,20 +3982,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", & @@ -4036,7 +4041,7 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, ! note: same logic/threshold used in global_cycle to produce ! mask on model grid. - call read_vardata(dsfg, 'slc1', values_2d, errcode=iret) + call read_vardata(dsfg, 'soill1', values_2d, errcode=iret) mask = 0 do j=1,nlats diff --git a/src/enkf/inflation.f90 b/src/enkf/inflation.f90 index 225967028c..51d1dd1106 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) diff --git a/src/enkf/readconvobs.f90 b/src/enkf/readconvobs.f90 index a5383069a1..65db770b6d 100644 --- a/src/enkf/readconvobs.f90 +++ b/src/enkf/readconvobs.f90 @@ -24,7 +24,6 @@ module readconvobs ! reflectivity and radial velocity assimilation. POC: xuguang.wang@ou.edu ! 2017-12-13 shlyaeva - added netcdf diag read/write capability ! 2019-03-21 CAPS(C. Tong) - added direct reflectivity DA capability -! 2022-03-23 draper - added option to not scale qobs by forecast qsat. ! ! attributes: ! language: f95 diff --git a/src/enkf/statevec.f90 b/src/enkf/statevec.f90 index 5ad70346aa..44ad5df9b4 100644 --- a/src/enkf/statevec.f90 +++ b/src/enkf/statevec.f90 @@ -136,7 +136,7 @@ subroutine init_statevec() do i = 1, ns2d if (getindex(vars2d_supported, svars2d(i))<0) then if (nproc .eq. 0) then - print *,'Error: 2D variable ', svars2d(i), ' is not supported in current version.' + print *,'Error: state 2D variable ', svars2d(i), ' is not supported in current version.' print *,'Supported variables: ', vars2d_supported endif call stop2(502) @@ -145,7 +145,7 @@ subroutine init_statevec() do i = 1, ns3d if (getindex(vars3d_supported, svars3d(i))<0) then if (nproc .eq. 0) then - print *,'Error: 3D variable ', svars3d(i), ' is not supported in current version.' + print *,'Error: state 3D variable ', svars3d(i), ' is not supported in current version.' print *,'Supported variables: ', vars3d_supported endif call stop2(502) diff --git a/src/gsi/read_prepbufr.f90 b/src/gsi/read_prepbufr.f90 index 87d5aa4bd8..526dc9ee78 100644 --- a/src/gsi/read_prepbufr.f90 +++ b/src/gsi/read_prepbufr.f90 @@ -1669,13 +1669,15 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& pmq(k)=nint(qcmark(8,k)) end do -! 181, 183, 187, and 188 are the screen-level obs over land - global_2m_land = ( (kx==181 .or. kx==183 .or. kx==188 .or. kx==188 ) .and. hofx_2m_sfcfile ) +! 187, 181, and 183 are the screen-level obs over land +! note: don't need the hofx_2m_sfcfile if set usage in convinfo, and qm updated in the input file + global_2m_land = ( (kx==187 .or. kx==181 .or. kx==183) .and. hofx_2m_sfcfile ) ! If temperature ob, extract information regarding virtual ! versus sensible temperature if(tob) then ! use tvirtual if tsensible flag not set, and not in either 2Dregional or global_2m DA mode + ! for now, keeping 2m obs as sensible, for global system. if ( (.not. tsensible) .and. .not. (twodvar_regional .or. global_2m_land) ) then do k=1,levs @@ -1979,18 +1981,17 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& pqm(k)=2 ! otherwise, type 183 will be discarded. qm=2 tqm(k)=2 - if (kx==187) obserr(3,k)=2.2 - if (kx==181) obserr(3,k)=1.5 - if (kx==183) obserr(3,k)=2.6 + if (kx==187) obserr(3,k)=2.0_r_double + if (kx==181) obserr(3,k)=2.0_r_double + if (kx==183) obserr(3,k)=2.0_r_double endif if (qob .and. qm == 9 ) then qm = 2 ! qob err specified as fraction of qsat, multiplied by 10. - if (kx==187) obserr(2,k)=1.0 - if (kx==181) obserr(2,k)=1.0 - if (kx==183) obserr(2,k)=1.0 + if (kx==187) obserr(2,k)=1.0_r_double + if (kx==181) obserr(2,k)=1.0_r_double + if (kx==183) obserr(2,k)=1.0_r_double endif - endif ! Set usage variable usage = zero diff --git a/src/gsi/setupq.f90 b/src/gsi/setupq.f90 index ad6d727ce9..993d57046a 100644 --- a/src/gsi/setupq.f90 +++ b/src/gsi/setupq.f90 @@ -376,7 +376,7 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav do k=1,nobs ikx=nint(data(ikxx,k)) itype=ictype(ikx) - landsfctype =( itype==181 .or. itype==183 .or. itype==187 .or. itype==188 ) + landsfctype =( itype==181 .or. itype==183 .or. itype==187 ) do l=k+1,nobs if (twodvar_regional .or. (hofx_2m_sfcfile .and. landsfctype) ) then duplogic=data(ilat,k) == data(ilat,l) .and. & diff --git a/src/gsi/setupt.f90 b/src/gsi/setupt.f90 index 8d1c308d7f..2f1f57f583 100644 --- a/src/gsi/setupt.f90 +++ b/src/gsi/setupt.f90 @@ -448,7 +448,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav do k=1,nobs ikx=nint(data(ikxx,k)) itype=ictype(ikx) - landsfctype =( itype==181 .or. itype==183 .or. itype==187 .or. itype==188 ) + landsfctype =( itype==181 .or. itype==183 .or. itype==187 ) do l=k+1,nobs if (twodvar_regional .or. (hofx_2m_sfcfile .and. landsfctype) ) then duplogic=data(ilat,k) == data(ilat,l) .and. & @@ -483,9 +483,6 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! Run a buddy-check ! Note: buddy check crashes for hofx_2m_sfcfile option. -! Ccurrent params have buddy radius of 108 km, max diff of 8 K. -! The gross error check removes O-F > 7., so this is probably removing -! most obs that fail the buddy check already if (twodvar_regional .and. buddycheck_t) call buddy_check_t(is,data,luse,mype,nele,nobs,muse,buddyuse) ! If requested, save select data for output to diagnostic file