From e55a937bea4ba780b03e14c7098421b2ab1fe109 Mon Sep 17 00:00:00 2001 From: RussTreadon-NOAA <26926959+RussTreadon-NOAA@users.noreply.github.com> Date: Sat, 11 Feb 2023 14:04:08 -0500 Subject: [PATCH 1/2] Regression test updates: global_4dvar bug fix, oom fix, enhance error checking (#532) --- regression/global_4dvar.sh | 12 ++++----- regression/regression_driver.sh | 1 - regression/regression_param.sh | 4 +-- regression/regression_test.sh | 17 +++++++++--- regression/regression_test_enkf.sh | 43 ++++++++++++++++++++---------- 5 files changed, 50 insertions(+), 27 deletions(-) diff --git a/regression/global_4dvar.sh b/regression/global_4dvar.sh index ec2a8396e0..cac0d28e6b 100755 --- a/regression/global_4dvar.sh +++ b/regression/global_4dvar.sh @@ -98,7 +98,7 @@ SINGLEOB="" # bftab_sst= bufr table for sst ONLY needed for sst retrieval (retrieval=.true.) # aeroinfo = text file with information about assimilation of aerosol data -anavinfo=$fixgsi/global_anavinfo.l${LEVS}.txt +anavinfo=$fixgsi/global_anavinfo_qlqi.l${LEVS}.txt berror=$fixgsi/Big_Endian/global_berror.l${LEVS}y${NLAT}.f77 locinfo=$fixgsi/global_hybens_info.l${LEVS}.txt satinfo=$fixgsi/global_satinfo.txt @@ -311,9 +311,9 @@ $gsi_namelist EOF cp gsiparm.anl gsiparm.anl.obsvr -echo "run gsi now" +echo "run gsi observer" eval "$APRUN $tmpdir/gsi.x < gsiparm.anl > stdout.obsvr 2>&1" -rc=$? +ra=$? # Run gsi identity model 4dvar under Parallel Operating Environment (poe) on NCEP IBM rm -f siganl sfcanl.gsi satbias_out fort.2* @@ -331,8 +331,8 @@ cat < gsiparm.anl $gsi_namelist EOF -echo "run gsi now" +echo "run gsi 4dvar" eval "$APRUN $tmpdir/gsi.x < gsiparm.anl > stdout 2>&1" -rc=$? - +rb=$? +rc=$((ra+rb)) exit $rc diff --git a/regression/regression_driver.sh b/regression/regression_driver.sh index 621ccbf485..e1d3b18dc7 100755 --- a/regression/regression_driver.sh +++ b/regression/regression_driver.sh @@ -47,7 +47,6 @@ for jn in `seq ${RSTART} ${REND}`; do $scripts/regression_wait.sh ${job[$jn]} ${rcname} $check_resource rc=$? if [ $rc -ne 0 ]; then - rm -f ${rcname} exit 1 fi done diff --git a/regression/regression_param.sh b/regression/regression_param.sh index 85f2949bfe..a2808ddfc0 100755 --- a/regression/regression_param.sh +++ b/regression/regression_param.sh @@ -150,8 +150,8 @@ case $regtest in topts[1]="0:15:00" ; popts[1]="20/1/" ; ropts[1]="/1" topts[2]="0:15:00" ; popts[2]="20/2/" ; ropts[2]="/1" elif [[ "$machine" = "Orion" ]]; then - topts[1]="0:15:00" ; popts[1]="4/4/" ; ropts[1]="/1" - topts[2]="0:15:00" ; popts[2]="6/6/" ; ropts[2]="/1" + topts[1]="0:15:00" ; popts[1]="20/1/" ; ropts[1]="/1" + topts[2]="0:15:00" ; popts[2]="20/2/" ; ropts[2]="/1" elif [[ "$machine" = "Jet" ]]; then topts[1]="0:15:00" ; popts[1]="4/4/" ; ropts[1]="/1" topts[2]="0:15:00" ; popts[2]="6/6/" ; ropts[2]="/1" diff --git a/regression/regression_test.sh b/regression/regression_test.sh index deb34ff244..0bcb9f4d90 100755 --- a/regression/regression_test.sh +++ b/regression/regression_test.sh @@ -576,11 +576,20 @@ mkdir -p $vfydir $ncp $output $vfydir/ +# Final check for any failed tests +count=$(grep -i "fail" $output |wc -l) +if [ $count -gt 0 ]; then + (( failed_test = $failed_test + $count )) +fi + +# Remove job log files is no failures detected cd $scripts -rm -f ${exp1}.out -rm -f ${exp2}.out -rm -f ${exp3}.out -rm -f ${exp2_scale}.out +if [ $count -eq 0 ]; then + rm -f ${exp1}.out + rm -f ${exp2}.out + rm -f ${exp3}.out + rm -f ${exp2_scale}.out +fi if [[ "$clean" = ".true." ]]; then rm -rf $savdir diff --git a/regression/regression_test_enkf.sh b/regression/regression_test_enkf.sh index 213ee726da..f52a5d451f 100755 --- a/regression/regression_test_enkf.sh +++ b/regression/regression_test_enkf.sh @@ -34,7 +34,7 @@ maxtime=1200 maxmem=${maxmem:-3400000} # set in regression_param maxmem=$((${memnode:-64}*1024*1024)) -# Copy stdout and sanl files +# Copy stdout and incr files # from $savdir to $tmpdir list="$exp1 $exp2 $exp3" for exp in $list; do @@ -43,7 +43,7 @@ for exp in $list; do imem=1 while [[ $imem -le $nmem ]]; do member="_mem"`printf %03i $imem` - $ncp $savdir/$exp/sanl_${global_adate}_fhr06$member $tmpdir/sanl$member.$exp + $ncp $savdir/$exp/incr_${global_adate}_fhr06$member $tmpdir/incr$member.$exp (( imem = $imem + 1 )) done done @@ -282,10 +282,13 @@ nmem=10 imem=1 while [[ $imem -le $nmem ]]; do member="_mem"`printf %03i $imem` - if ! cmp -s sanl$member.${exp1} sanl$member.${exp2} -then - echo 'sanl'$member'.'${exp1}' sanl'$member'.'${exp2}' are NOT identical' -fi + ncdump incr$member.${exp1} > incr$member.${exp1}.out + ncdump incr$member.${exp2} > incr$member.${exp2}.out + if [ ! diff incr$member.${exp1}.out incr$member.${exp2}.out ]; then + echo 'incr'$member'.'${exp1}' incr'$member'.'${exp2}' are NOT identical' + else + rm -f incr$member.${exp1}.out incr$member.${exp2}.out + fi (( imem = $imem + 1 )) done echo @@ -379,11 +382,14 @@ else imem=1 while [[ $imem -le $nmem ]]; do member="_mem"`printf %03i $imem` - if ! cmp -s sanl$member.${exp1} sanl$member.${exp3} - then - echo 'sanl'$member'.'${exp1}' sanl'$member'.'${exp3}' are NOT identical' + ncdump incr$member.${exp1} > incr$member.${exp1}.out + ncdump incr$member.${exp3} > incr$member.${exp3}.out + if [ ! diff incr$member.${exp1}.out incr$member.${exp3}.out ]; then + echo 'incr'$member'.'${exp1}' incr'$member'.'${exp3}' are NOT identical' + else + rm -f incr$member.${exp1}.out incr$member.${exp3}.out fi - (( imem = $imem + 1 )) + (( imem = $imem + 1 )) done echo } >> $output @@ -411,11 +417,20 @@ mkdir -p $vfydir $ncp $output $vfydir/ +# Final check for any failed tests +count=$(grep -i "fail" $output |wc -l) +if [ $count -gt 0 ]; then + (( failed_test = $failed_test + $count )) +fi + +# Remove job log files is no failures detected cd $scripts -rm -f ${exp1}.out -rm -f ${exp2}.out -rm -f ${exp3}.out -rm -f ${exp2_scale}.out +if [ $count -eq 0 ]; then + rm -f ${exp1}.out + rm -f ${exp2}.out + rm -f ${exp3}.out + rm -f ${exp2_scale}.out +fi if [[ "$clean" = ".true." ]]; then rm -rf $savdir From 8857995a2467a2f0c608a10d08ca24fdf8589fcc Mon Sep 17 00:00:00 2001 From: Yongming-OUMAP <66268887+Wangy1111@users.noreply.github.com> Date: Mon, 13 Feb 2023 17:12:15 -0600 Subject: [PATCH 2/2] GitHub Issue NOAA-EMC/GSI#468 Enhancements to SDL and VDL for simultaneous multiscale EnVar and parallel ensemble IO for EnVar for FV3-LAM (#504) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The following capabilities developed by OU MAP lab are included (1) Further development for simultaneous multiscale EnVar for both global and regional DA (1a) spatial scale-dependent localization (SDL; contributed by Ting Lei and Daryl Kleist/EMC) is implemented in EnVar as described in Huang et al 2021, MWR for the global NWP application. (1b) variable-dependent localization (VDL) method by Wang and Wang 2022, JAMES is implemented in EnVar. (2)Development of parallel ensemble IO for EnVar for FV3-LAM  Implement an approach to simultaneously read in all ensemble members for EnVar. Specifically, parallel ensemble IO for both conventional and radar EnVar for FV3-LAM is implemented by reading in all ensemble members simultaneously. (3) Direct assimilation of radar reflectivity for EnVar for RRFS The direct radar reflectivity assimilation approach by Wang and Wang 2017, MWR is implemented and tested for FV3-LAM. Fixes #468 --- src/enkf/gridinfo_fv3reg.f90 | 3 +- src/gsi/apply_scaledepwgts.f90 | 173 ++++++ src/gsi/cplr_get_fv3_regional_ensperts.f90 | 519 +++++++++++++++++- src/gsi/general_specmod.f90 | 31 ++ src/gsi/get_gefs_ensperts_dualres.f90 | 32 ++ src/gsi/gsi_files.cmake | 1 + src/gsi/gsi_rfv3io_mod.f90 | 593 ++++++++++++++++++++- src/gsi/gsimod.F90 | 44 +- src/gsi/hybrid_ensemble_isotropic.F90 | 213 ++++++-- src/gsi/hybrid_ensemble_parameters.f90 | 21 + src/gsi/obsmod.F90 | 2 +- src/gsi/read_dbz_nc.f90 | 8 +- src/gsi/setupdbz.f90 | 48 +- src/gsi/update_guess.f90 | 6 +- 14 files changed, 1576 insertions(+), 118 deletions(-) create mode 100644 src/gsi/apply_scaledepwgts.f90 diff --git a/src/enkf/gridinfo_fv3reg.f90 b/src/enkf/gridinfo_fv3reg.f90 index 53e5f5b3de..ef5b242901 100644 --- a/src/enkf/gridinfo_fv3reg.f90 +++ b/src/enkf/gridinfo_fv3reg.f90 @@ -43,7 +43,8 @@ module gridinfo ! !$$$ -use mpisetup, only: nproc, mpi_integer, mpi_real4, mpi_comm_world,mpi_status +use mpisetup, only: nproc, mpi_integer, mpi_real4,mpi_status +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 diff --git a/src/gsi/apply_scaledepwgts.f90 b/src/gsi/apply_scaledepwgts.f90 new file mode 100644 index 0000000000..62c455e011 --- /dev/null +++ b/src/gsi/apply_scaledepwgts.f90 @@ -0,0 +1,173 @@ +!$$$ program documentation block +! +! program history: +! +! 2018-03-28 T. Lei and D. Kleist - consoliated and added codes +! for the scale dependent scale localization scheme +! +!$$$ end documentation block + +function fwgtofwvlen (rvlft,rvrgt,rcons,rlen,rinput) +!$$$ subprogram documentation block +! +! subprogram: fwgtofwvlen +! +! abstract: Calculation of spectral filter functions +! +!$$$ end documentation block + + use kinds, only: r_kind,i_kind,r_single + implicit none + + real(r_kind),intent(in) :: rvlft,rvrgt,rcons,rlen,rinput + real(r_kind) :: fwgtofwvlen + real(r_kind) :: rlen1,rtem1,rconshalf + + rlen1=rlen/10.0_r_kind ! rlen corresponds to a (-5,5) region + rconshalf=0.5_r_kind*rcons + if(rinput > rvlft .and. rinput < rvrgt) then + fwgtofwvlen=rcons + else + rtem1=min(abs(rinput-rvlft),abs(rinput-rvrgt)) + fwgtofwvlen=rconshalf*(1.0_r_kind+tanh(5.0_r_kind-rtem1/rlen1)) + endif + +end function fwgtofwvlen +! . . . . +subroutine init_mult_spc_wgts(jcap_in) +!$$$ subprogram documentation block +! +! subprogram: init_mult_spc_wgts +! +!$$$ end documentation block + + use kinds, only: r_kind,i_kind,r_single + use hybrid_ensemble_parameters,only: s_ens_hv,sp_loc,grd_ens,grd_loc,sp_ens + use hybrid_ensemble_parameters,only: n_ens,p_sploc2ens,grd_sploc + use hybrid_ensemble_parameters,only: use_localization_grid + use gridmod,only: use_sp_eqspace + use general_specmod, only: general_init_spec_vars + use constants, only: zero,half,one,two,three,rearth,pi + use constants, only: rad2deg + use mpimod, only: mype + use general_sub2grid_mod, only: general_sub2grid_create_info + use egrid2agrid_mod,only: g_create_egrid2agrid + use general_sub2grid_mod, only: sub2grid_info + use gsi_io, only: verbose + use hybrid_ensemble_parameters, only: nsclgrp + use hybrid_ensemble_parameters, only: spc_multwgt,spcwgt_params,i_ensloccov4scl + implicit none + + integer(i_kind),intent(in ) :: jcap_in + real(r_kind),allocatable :: totwvlength(:) + + integer(i_kind) i,ii,j,k,l,n,kk,nsigend + integer(i_kind) ig + real(r_kind) rwv0,rtem1,rtem2 + real (r_kind):: fwgtofwvlen + integer(i_kind) :: l_sum_spc_weights + + ! Spectral scale decomposition is differernt between SDL-cross and SDL-nocross + if( i_ensloccov4scl == 1 )then + l_sum_spc_weights = 1 + else + l_sum_spc_weights = 0 + end if + + allocate(totwvlength(jcap_in)) + + rwv0=2*pi*rearth*0.001_r_kind + do i=1,jcap_in + totwvlength(i)= rwv0/real(i) + enddo + do i=1,jcap_in + rtem1=0 + do ig=1,nsclgrp + if(ig /= 2) then + spc_multwgt(i,ig)=fwgtofwvlen(spcwgt_params(1,ig),spcwgt_params(2,ig),& + spcwgt_params(3,ig),spcwgt_params(4,ig),totwvlength(i)) + if(l_sum_spc_weights == 0 ) then + rtem1=rtem1+spc_multwgt(i,ig) + else + rtem1=rtem1+spc_multwgt(i,ig)*spc_multwgt(i,ig) + endif + endif + enddo + rtem2 =1.0_r_kind - rtem1 + if(abs(rtem2) >= zero) then + + if(l_sum_spc_weights == 0 ) then + spc_multwgt(i,2)=rtem2 + else + spc_multwgt(i,2)=sqrt(rtem2) + endif + endif + enddo + spc_multwgt=max(spc_multwgt,0.0_r_kind) + + deallocate(totwvlength) + return +end subroutine init_mult_spc_wgts + +subroutine apply_scaledepwgts(grd_in,sp_in,wbundle,spwgts,wbundle2) +! +! Program history log: +! 2017-03-30 J. Kay, X. Wang - copied from Kleist's apply_scaledepwgts and +! add the calculation of scale-dependent weighting for mixed resolution ensemble +! POC: xuguang.wang@ou.edu +! + use constants, only: one + use control_vectors, only: nrf_var,cvars2d,cvars3d,control_vector + use kinds, only: r_kind,i_kind + use kinds, only: r_single + use mpimod, only: mype,nvar_id,levs_id + use hybrid_ensemble_parameters, only: oz_univ_static + use general_specmod, only: general_spec_multwgt + use gsi_bundlemod, only: gsi_bundle + use general_sub2grid_mod, only: general_sub2grid,general_grid2sub + use general_specmod, only: spec_vars + use general_sub2grid_mod, only: sub2grid_info + use mpimod, only: mpi_comm_world,mype,npe,ierror + use file_utility, only : get_lun + implicit none + +! Declare passed variables + type(gsi_bundle),intent(in) :: wbundle + type(gsi_bundle),intent(inout) :: wbundle2 + type(spec_vars),intent (in):: sp_in + type(sub2grid_info),intent(in)::grd_in + real(r_kind),dimension(0:sp_in%jcap),intent(in):: spwgts + +! Declare local variables + integer(i_kind) ii,kk + integer(i_kind) i,j,lunit + + real(r_kind),dimension(grd_in%lat2,grd_in%lon2):: slndt,sicet,sst + real(r_kind),dimension(grd_in%nlat*grd_in%nlon*grd_in%nlevs_alloc) :: hwork + real(r_kind),dimension(grd_in%nlat,grd_in%nlon,grd_in%nlevs_alloc) :: work + real(r_kind),dimension(sp_in%nc):: spc1 + character*64 :: fname1 + character*5:: varname1 + +! Beta1 first +! Get from subdomains to + call general_sub2grid(grd_in,wbundle%values,hwork) + work=reshape(hwork,(/grd_in%nlat,grd_in%nlon,grd_in%nlevs_alloc/)) + + do kk=1,grd_in%nlevs_alloc +! Transform from physical space to spectral space + call general_g2s0(grd_in,sp_in,spc1,work(:,:,kk)) + +! Apply spectral weights + call general_spec_multwgt(sp_in,spc1,spwgts) +! Transform back to physical space + call general_s2g0(grd_in,sp_in,spc1,work(:,:,kk)) + + end do + +! Transfer work back to subdomains + hwork=reshape(work,(/grd_in%nlat*grd_in%nlon*grd_in%nlevs_alloc/)) + call general_grid2sub(grd_in,hwork,wbundle2%values) + + return +end subroutine apply_scaledepwgts diff --git a/src/gsi/cplr_get_fv3_regional_ensperts.f90 b/src/gsi/cplr_get_fv3_regional_ensperts.f90 index f99ca39790..6e94b29c6c 100644 --- a/src/gsi/cplr_get_fv3_regional_ensperts.f90 +++ b/src/gsi/cplr_get_fv3_regional_ensperts.f90 @@ -8,14 +8,18 @@ module get_fv3_regional_ensperts_mod procedure, pass(this) :: get_fv3_regional_ensperts => get_fv3_regional_ensperts_run procedure, pass(this) :: ens_spread_dualres_regional => ens_spread_dualres_regional_fv3_regional procedure, pass(this) :: general_read_fv3_regional + procedure, pass(this) :: general_read_fv3_regional_parallel_over_ens + procedure, pass(this) :: parallel_read_fv3_step2 + procedure, nopass :: fill_regional_2d end type get_fv3_regional_ensperts_class - type(sub2grid_info):: grd_fv3lam_ens_dynvar_io_nouv,grd_fv3lam_ens_tracer_io_nouv,grd_fv3lam_ens_uv + type(sub2grid_info):: grd_fv3lam_ens_dynvar_io_nouv,grd_fv3lam_ens_tracer_io_nouv,grd_fv3lam_ens_uv,grd_fv3lam_ens_phyvar_io_nouv character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_ens_io_dynmetvars3d_nouv ! copy of cvars3d excluding uv 3-d fields character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_ens_io_tracermetvars3d_nouv ! copy of cvars3d excluding uv 3-d fields character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_ens_io_dynmetvars2d_nouv character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_ens_io_tracermetvars2d_nouv + character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_ens_io_phymetvars3d_nouv contains subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) @@ -46,8 +50,8 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) use kinds, only: r_kind,i_kind,r_single use constants, only: zero,one,half,zero_single,rd_over_cp,one_tenth - use mpimod, only: mpi_comm_world,ierror,mype - use hybrid_ensemble_parameters, only: n_ens,grd_ens + use mpimod, only: mpi_comm_world,ierror,mype,npe + use hybrid_ensemble_parameters, only: n_ens,grd_ens,parallelization_over_ensmembers use hybrid_ensemble_parameters, only: l_both_fv3sar_gfs_ens, n_ens_gfs,n_ens_fv3sar use hybrid_ensemble_parameters, only: ntlevs_ens,ensemble_path use control_vectors, only: cvars2d,cvars3d,nc2d,nc3d @@ -69,6 +73,8 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) use gsi_rfv3io_mod, only: fv3lam_io_dynmetvars2d_nouv,fv3lam_io_tracermetvars2d_nouv use netcdf , only: nf90_open, nf90_close,nf90_nowrite,nf90_inquire,nf90_format_netcdf4 use netcdf_mod , only: nc_check + use gsi_rfv3io_mod, only: fv3lam_io_phymetvars3d_nouv + use obsmod, only: if_model_dbz implicit none @@ -79,7 +85,11 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig):: u,v,tv,oz,rh real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2):: ps - real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)::w,ql,qi,qr,qg,qs,qnr + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)::w,ql,qi,qr,qg,qs,qnr,dbz + real(r_kind),dimension(:,:,:),allocatable :: gg_u,gg_v,gg_tv,gg_rh + real(r_kind),dimension(:,:,:),allocatable :: gg_w,gg_dbz,gg_qr,gg_qs, & + gg_qi,gg_qg,gg_oz,gg_cwmr + real(r_kind),dimension(:,:),allocatable :: gg_ps real(r_single),pointer,dimension(:,:,:):: w3 =>NULL() real(r_single),pointer,dimension(:,:):: w2 =>NULL() @@ -96,9 +106,9 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) integer(i_kind):: i,j,k,n,mm1,istatus integer(i_kind):: ndynvario2d,ntracerio2d - integer(r_kind):: ndynvario3d,ntracerio3d + integer(r_kind):: ndynvario3d,ntracerio3d,nphyvario3d integer(i_kind):: inner_vars,numfields - integer(i_kind):: ilev,ic2,ic3 + integer(i_kind):: ilev,ic2,ic3,iope integer(i_kind):: m integer(i_kind)::loc_id,ncfmt @@ -125,8 +135,10 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) !clt setup varnames for IO ndynvario2d=0 ntracerio2d=0 + nphyvario3d=0 ndynvario3d=size(fv3lam_io_dynmetvars3d_nouv) ntracerio3d=size(fv3lam_io_tracermetvars3d_nouv) + nphyvario3d=size(fv3lam_io_phymetvars3d_nouv) if (allocated(fv3lam_io_dynmetvars2d_nouv))then ndynvario2d=size(fv3lam_io_dynmetvars2d_nouv) endif @@ -137,6 +149,10 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) allocate(fv3lam_ens_io_tracermetvars3d_nouv(ndynvario3d)) fv3lam_ens_io_dynmetvars3d_nouv=fv3lam_io_dynmetvars3d_nouv fv3lam_ens_io_tracermetvars3d_nouv=fv3lam_io_tracermetvars3d_nouv + if ( nphyvario3d > 0 )then + allocate(fv3lam_ens_io_phymetvars3d_nouv(nphyvario3d)) + fv3lam_ens_io_phymetvars3d_nouv=fv3lam_io_phymetvars3d_nouv + end if if (ndynvario2d > 0 ) then allocate(fv3lam_ens_io_dynmetvars2d_nouv(ndynvario2d)) fv3lam_ens_io_dynmetvars2d_nouv=fv3lam_io_dynmetvars2d_nouv @@ -168,6 +184,24 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) call general_sub2grid_create_info(grd_fv3lam_ens_dynvar_io_nouv,inner_vars,grd_ens%nlat,& grd_ens%nlon,grd_ens%nsig,numfields,regional,names=names,lnames=lnames) + if( nphyvario3d > 0 )then + inner_vars=1 + numfields=inner_vars*(nphyvario3d*grd_ens%nsig) + deallocate(lnames,names) + allocate(lnames(1,numfields),names(1,numfields)) + ilev=1 + do i=1,nphyvario3d + do k=1,grd_ens%nsig + lnames(1,ilev)=k + names(1,ilev)=fv3lam_ens_io_phymetvars3d_nouv(i) + ilev=ilev+1 + enddo + enddo + + call general_sub2grid_create_info(grd_fv3lam_ens_phyvar_io_nouv,inner_vars,grd_ens%nlat,& + grd_ens%nlon,grd_ens%nsig,numfields,regional,names=names,lnames=lnames) + end if + inner_vars=1 numfields=inner_vars*(ntracerio3d*grd_ens%nsig+ntracerio2d) deallocate(lnames,names) @@ -254,16 +288,73 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) mm1=mype+1 kap1=rd_over_cp+one kapr=one/rd_over_cp + + if( parallelization_over_ensmembers ) then + if(n_ens_fv3sar>npe) then + parallelization_over_ensmembers=.false. +130 format('Disabling parallelization_over_ensmembers because number of ensemble members (',I0,') is greater than number of MPI ranks (',I0,').') + if(mype==0) then + write(6,130) n_ens_fv3sar,npe + endif + endif + endif ! parallelization_over_ensmembers + + if(parallelization_over_ensmembers .and. mype==0) then + write(6,'(I0,A)') mype,': will read ensemble data in parallel (parallelization_over_ensmembers=.true.)' + endif + + if( parallelization_over_ensmembers )then + do n=1,n_ens_fv3sar + write(ensfilenam_str,22) trim(adjustl(ensemble_path)),ens_fhrlevs(m),n +22 format(a,'fv3SAR',i2.2,'_ens_mem',i3.3) + iope=(n-1)*npe/n_ens_fv3sar + ! DEFINE INPUT FILE NAME + fv3_filename%grid_spec=trim(ensfilenam_str)//'-fv3_grid_spec' + fv3_filename%ak_bk=trim(ensfilenam_str)//'-fv3_akbk' + fv3_filename%dynvars=trim(ensfilenam_str)//'-fv3_dynvars' + fv3_filename%tracers=trim(ensfilenam_str)//"-fv3_tracer" + fv3_filename%phyvars=trim(ensfilenam_str)//'-fv3_phyvars' + fv3_filename%sfcdata=trim(ensfilenam_str)//"-fv3_sfcdata" + fv3_filename%couplerres=trim(ensfilenam_str)//"-coupler.res" + + + if( mype==iope) then + allocate(gg_u(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) + allocate(gg_v(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) + allocate(gg_tv(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) + allocate(gg_rh(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) + allocate(gg_oz(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) + allocate(gg_ps(grd_ens%nlat,grd_ens%nlon)) + if ( .not. if_model_dbz ) then + call this%general_read_fv3_regional_parallel_over_ens(iope,fv3_filename,gg_ps,gg_u,gg_v,gg_tv,gg_rh,gg_oz) + else + allocate(gg_w(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) + allocate(gg_dbz(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) + allocate(gg_qr(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) + allocate(gg_qs(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) + allocate(gg_qi(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) + allocate(gg_qg(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) + allocate(gg_cwmr(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) + call this%general_read_fv3_regional_parallel_over_ens(iope,fv3_filename,gg_ps,gg_u,gg_v,gg_tv,gg_rh,gg_oz, & + g_ql=gg_cwmr,g_qi=gg_qi,g_qr=gg_qr,g_qs=gg_qs,g_qg=gg_qg,g_w=gg_w,g_dbz=gg_dbz) + end if + end if + end do + if(mype==0) then + write(6,'(I0,A)') mype,': reading ensemble data in parallel is done (parallelization_over_ensmembers=.true.)' + endif + end if + call MPI_Barrier(mpi_comm_world,ierror) ! ! LOOP OVER ENSEMBLE MEMBERS do n_fv3sar=1,n_ens_fv3sar n=n_ens_gfs+n_fv3sar write(ensfilenam_str,22) trim(adjustl(ensemble_path)),ens_fhrlevs(m),n_fv3sar -22 format(a,'fv3SAR',i2.2,'_ens_mem',i3.3) ! DEFINE INPUT FILE NAME fv3_filename%grid_spec=trim(ensfilenam_str)//'-fv3_grid_spec' !exmaple thinktobe fv3_filename%ak_bk=trim(ensfilenam_str)//'-fv3_akbk' fv3_filename%dynvars=trim(ensfilenam_str)//'-fv3_dynvars' + fv3_filename%phyvars=trim(ensfilenam_str)//'-fv3_phyvars' fv3_filename%tracers=trim(ensfilenam_str)//"-fv3_tracer" fv3_filename%sfcdata=trim(ensfilenam_str)//"-fv3_sfcdata" fv3_filename%couplerres=trim(ensfilenam_str)//"-coupler.res" @@ -299,15 +390,52 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) endif ! ! READ ENEMBLE MEMBERS DATA - if (mype == 0) write(6,'(a,a)') & - 'CALL READ_FV3_REGIONAL_ENSPERTS FOR ENS DATA with the filename str : ',trim(ensfilenam_str) - if (.not.l_use_dbz_directDA ) then ! Read additional hydrometers and w for dirZDA - call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz) - else - call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz, & - g_ql=ql,g_qi=qi,g_qr=qr,g_qs=qs,g_qg=qg,g_qnr=qnr,g_w=w) + if( .not. parallelization_over_ensmembers )then + if (mype == 0) write(6,'(a,a)') & + 'CALL READ_FV3_REGIONAL_ENSPERTS FOR ENS DATA with the filename str : ',trim(ensfilenam_str) + if (.not. (l_use_dbz_directDA .or. if_model_dbz) ) then ! Read additional hydrometers and w for dirZDA + call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz) + else + if( l_use_dbz_directDA ) then + call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz, & + g_ql=ql,g_qi=qi,g_qr=qr,g_qs=qs,g_qg=qg,g_qnr=qnr,g_w=w) + else if( if_model_dbz )then + call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz, & + g_ql=ql,g_qi=qi,g_qr=qr,g_qs=qs,g_qg=qg,g_qnr=qnr,g_w=w,g_dbz=dbz) + end if + end if end if + if( parallelization_over_ensmembers )then + iope=(n_fv3sar-1)*npe/n_ens_fv3sar + if(mype==iope) then + write(0,'(I0,A,I0,A)') mype,': scatter member ',n_fv3sar,' to other ranks...' + if( if_model_dbz )then + call this%parallel_read_fv3_step2(mype,iope,& + g_ps=ps,g_u=u,g_v=v,g_tv=tv,g_rh=rh,g_ql=ql,& + g_oz=oz,g_w=w,g_qr=qr,g_qs=qs,g_qi=qi,g_qg=qg,g_dbz=dbz,& + gg_ps=gg_ps,gg_tv=gg_tv,gg_u=gg_u,gg_v=gg_v,& + gg_rh=gg_rh,gg_w=gg_w,gg_dbz=gg_dbz,gg_qr=gg_qr,& + gg_qs=gg_qs,gg_qi=gg_qi,gg_qg=gg_qg,gg_ql=gg_cwmr) + else + call this%parallel_read_fv3_step2(mype,iope,& + g_ps=ps,g_u=u,g_v=v,g_tv=tv,g_rh=rh,g_ql=ql,g_oz=oz, & + gg_ps=gg_ps,gg_tv=gg_tv,gg_u=gg_u,gg_v=gg_v,gg_rh=gg_rh) + end if + else + if( if_model_dbz )then + call this%parallel_read_fv3_step2(mype,iope,& + g_ps=ps,g_u=u,g_v=v,g_tv=tv,g_rh=rh,g_ql=ql,& + g_oz=oz,g_w=w,g_qr=qr,g_qs=qs,g_qi=qi,g_qg=qg,g_dbz=dbz) + else + call this%parallel_read_fv3_step2(mype,iope,& + g_ps=ps,g_u=u,g_v=v,g_tv=tv,g_rh=rh,g_ql=ql,g_oz=oz) + endif + endif + + call MPI_Barrier(mpi_comm_world,ierror) + end if + ! SAVE ENSEMBLE MEMBER DATA IN COLUMN VECTOR do ic3=1,nc3d @@ -463,6 +591,16 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) end do end do end do + + case('dbz','DBZ') + do k=1,grd_ens%nsig + do i=1,grd_ens%lon2 + do j=1,grd_ens%lat2 + w3(j,i,k) = dbz(j,i,k) + x3(j,i,k)=x3(j,i,k)+dbz(j,i,k) + end do + end do + end do end select @@ -572,7 +710,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) end subroutine get_fv3_regional_ensperts_run subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g_rh,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,g_qnr,g_w) + g_ql,g_qi,g_qr,g_qs,g_qg,g_qnr,g_w,g_dbz) !$$$ subprogram documentation block ! first compied from general_read_arw_regional . . . . ! subprogram: general_read_fv3_regional read fv3sar model ensemble members @@ -623,7 +761,7 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g use hybrid_ensemble_parameters, only: grd_ens use directDA_radaruse_mod, only: l_use_cvpqx, cvpqx_pval, cld_nt_updt use directDA_radaruse_mod, only: l_cvpnr, cvpnr_pval - + use obsmod, only:if_model_dbz implicit none @@ -632,7 +770,7 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g class(get_fv3_regional_ensperts_class), intent(inout) :: this type (type_fv3regfilenameg) , intent (in) :: fv3_filenameginput real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),intent(out)::g_u,g_v,g_tv,g_rh,g_oz - real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),optional,intent(out)::g_ql,g_qi,g_qr + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),optional,intent(out)::g_ql,g_qi,g_qr,g_dbz real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),optional,intent(out)::g_qs,g_qg,g_qnr,g_w real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2),intent(out):: g_ps @@ -656,11 +794,13 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g character(len=24),parameter :: myname_ = 'general_read_fv3_regional' type(gsi_bundle) :: gsibundle_fv3lam_ens_dynvar_nouv type(gsi_bundle) :: gsibundle_fv3lam_ens_tracer_nouv + type(gsi_bundle) :: gsibundle_fv3lam_ens_phyvar_nouv type(gsi_grid):: grid_ens character(len=:),allocatable :: grid_spec !='fv3_grid_spec' character(len=:),allocatable :: ak_bk !='fv3_akbk' character(len=:),allocatable :: dynvars !='fv3_dynvars' + character(len=:),allocatable :: phyvars !='fv3_phyvars' character(len=:),allocatable :: tracers !='fv3_tracer' character(len=:),allocatable :: sfcdata !='fv3_sfcdata' character(len=:),allocatable :: couplerres!='coupler.res' @@ -676,6 +816,7 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g grid_spec=fv3_filenameginput%grid_spec ak_bk=fv3_filenameginput%ak_bk dynvars=fv3_filenameginput%dynvars + phyvars=fv3_filenameginput%phyvars tracers=fv3_filenameginput%tracers sfcdata=fv3_filenameginput%sfcdata couplerres=fv3_filenameginput%couplerres @@ -701,7 +842,10 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g endif - + if(allocated(fv3lam_ens_io_phymetvars3d_nouv))then + call gsi_bundlecreate(gsibundle_fv3lam_ens_phyvar_nouv,grid_ens,'gsibundle_fv3lam_ens_phyvar_nouv',istatus, & + names3d=fv3lam_ens_io_phymetvars3d_nouv) + end if if(fv3sar_ensemble_opt == 0 ) then @@ -714,6 +858,10 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g fv3_filenameginput%dynvars,fv3_filenameginput) call gsi_fv3ncdf_read(grd_fv3lam_ens_tracer_io_nouv,gsibundle_fv3lam_ens_tracer_nouv,& fv3_filenameginput%tracers,fv3_filenameginput) + if( if_model_dbz ) then + call gsi_fv3ncdf_read(grd_fv3lam_ens_phyvar_io_nouv,gsibundle_fv3lam_ens_phyvar_nouv,& + fv3_filenameginput%phyvars,fv3_filenameginput) + end if else call gsi_fv3ncdf_read_v1(grd_fv3lam_ens_dynvar_io_nouv,gsibundle_fv3lam_ens_dynvar_nouv,& fv3_filenameginput%dynvars,fv3_filenameginput) @@ -724,14 +872,17 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_dynvar_nouv, 'tsen' ,g_tsen ,istatus );ier=ier+istatus call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'q' ,g_q ,istatus );ier=ier+istatus call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'oz' ,g_oz ,istatus );ier=ier+istatus - if (l_use_dbz_directDA) then + if (l_use_dbz_directDA .or. if_model_dbz) then call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'ql' ,g_ql ,istatus );ier=ier+istatus call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qi' ,g_qi ,istatus );ier=ier+istatus call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qr' ,g_qr ,istatus );ier=ier+istatus call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qs' ,g_qs ,istatus );ier=ier+istatus call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qg' ,g_qg ,istatus );ier=ier+istatus + if (l_use_dbz_directDA) & call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qnr',g_qnr ,istatus );ier=ier+istatus call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_dynvar_nouv, 'w' , g_w ,istatus );ier=ier+istatus + if( if_model_dbz )& + call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_phyvar_nouv, 'dbz' , g_dbz ,istatus );ier=ier+istatus end if @@ -834,11 +985,339 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g enddo call gsi_bundledestroy(gsibundle_fv3lam_ens_dynvar_nouv) call gsi_bundledestroy(gsibundle_fv3lam_ens_tracer_nouv) - + call gsi_bundledestroy(gsibundle_fv3lam_ens_phyvar_nouv) return end subroutine general_read_fv3_regional + subroutine general_read_fv3_regional_parallel_over_ens(this,iope,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g_rh,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,g_qnr,g_w,g_dbz) + !$$$ subprogram documentation block + ! first compied from general_read_arw_regional . . . . + ! subprogram: general_read_fv3_regional read fv3sar model ensemble members + ! prgmmr: Ting org: emc/ncep date: 2018 + ! + ! abstract: read ensemble members from the fv3sar model in "restart" or "cold start" netcdf format + ! for use with hybrid ensemble option. + ! + ! program history log: + ! 2018- Ting - intial versions + ! 2022-04-01 Y. Wang and X. Wang - read all fields for each member for + ! parallel ensemble IO capability + ! poc: xuguang.wang@ou.edu + ! + ! input argument list: + ! + ! output argument list: + ! + ! attributes: + ! language: f90 + ! machine: ibm RS/6000 SP + ! + !$$$ end documentation block + + use netcdf, only: nf90_nowrite + use netcdf, only: nf90_open,nf90_close + use netcdf, only: nf90_inq_dimid,nf90_inquire_dimension + use netcdf, only: nf90_inq_varid,nf90_inquire_variable,nf90_get_var + use kinds, only: r_kind,r_single,i_kind + use gridmod, only: eta1_ll,eta2_ll + use constants, only: zero,one,fv,zero_single,one_tenth,h300 + use hybrid_ensemble_parameters, only: grd_ens,q_hyb_ens + use hybrid_ensemble_parameters, only: fv3sar_ensemble_opt + + use mpimod, only: mpi_comm_world,mpi_rtype,mype + use netcdf_mod, only: nc_check + use gsi_rfv3io_mod,only: type_fv3regfilenameg + use gsi_rfv3io_mod,only:n2d + use constants, only: half,zero + use gsi_rfv3io_mod, only: gsi_fv3ncdf_read + use gsi_rfv3io_mod, only: gsi_fv3ncdf_read_v1 + use gsi_rfv3io_mod, only: gsi_fv3ncdf_readuv + use gsi_rfv3io_mod, only: gsi_fv3ncdf_readuv_v1 + use gsi_rfv3io_mod, only: gsi_fv3ncdf2d_read_v1 + use directDA_radaruse_mod, only: l_use_dbz_directDA + use gsi_bundlemod, only: gsi_gridcreate + use gsi_bundlemod, only: gsi_grid + use gsi_bundlemod, only: gsi_bundlecreate,gsi_bundledestroy + use gsi_bundlemod, only: gsi_bundlegetvar + use obsmod, only: if_model_dbz + use directDA_radaruse_mod, only: l_use_cvpqx, cvpqx_pval, cld_nt_updt + use directDA_radaruse_mod, only: l_cvpnr, cvpnr_pval + use gsi_rfv3io_mod, only: gsi_fv3ncdf_read_ens_parallel_over_ens,gsi_fv3ncdf_readuv_ens_parallel_over_ens + + + + implicit none +! +! Declare passed variables + class(get_fv3_regional_ensperts_class), intent(inout) :: this + integer(i_kind), intent (in) :: iope + type (type_fv3regfilenameg) , intent (in) :: fv3_filenameginput + real(r_kind),dimension(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig),intent(out)::g_u,g_v,g_tv,g_rh,g_oz + real(r_kind),dimension(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig),optional,intent(out)::g_ql,g_qi,g_qr,g_dbz + real(r_kind),dimension(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig),optional,intent(out)::g_qs,g_qg,g_qnr,g_w + + real(r_kind),dimension(grd_ens%nlat,grd_ens%nlon),intent(out):: g_ps + real(r_kind),dimension(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig+1) ::g_prsi + real(r_kind),dimension(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig) ::g_prsl ,g_tsen,g_q,g_delp +! +! Declare local parameters + real(r_kind),parameter:: r0_01 = 0.01_r_kind + real(r_kind),parameter:: r10 = 10.0_r_kind + real(r_kind),parameter:: r100 = 100.0_r_kind + ! +! Declare local variables + + integer(i_kind):: i,j,k,kp + + integer(i_kind) iderivative + + + logical ice + + character(len=24),parameter :: myname_ = 'general_read_fv3_regional' + + character(len=:),allocatable :: grid_spec !='fv3_grid_spec' + character(len=:),allocatable :: ak_bk !='fv3_akbk' + character(len=:),allocatable :: dynvars !='fv3_dynvars' + character(len=:),allocatable :: phyvars !='fv3_phyvars' + character(len=:),allocatable :: tracers !='fv3_tracer' + character(len=:),allocatable :: sfcdata !='fv3_sfcdata' + character(len=:),allocatable :: couplerres!='coupler.res' + integer (i_kind) ier,istatus + + + associate( this => this ) ! eliminates warning for unused dummy argument needed for binding + end associate + + if( mype == iope )then + grid_spec=fv3_filenameginput%grid_spec + ak_bk=fv3_filenameginput%ak_bk + dynvars=fv3_filenameginput%dynvars + phyvars=fv3_filenameginput%phyvars + tracers=fv3_filenameginput%tracers + sfcdata=fv3_filenameginput%sfcdata + couplerres=fv3_filenameginput%couplerres + + if(fv3sar_ensemble_opt == 0 ) then + call gsi_fv3ncdf_readuv_ens_parallel_over_ens(g_u,g_v,fv3_filenameginput,iope) + else + write(6,*) "Warning: we can only grab fields from restart files not cold start files for ensemble!" + endif + + if(fv3sar_ensemble_opt == 0) then + if (if_model_dbz) then + call gsi_fv3ncdf_read_ens_parallel_over_ens(fv3_filenameginput%dynvars,fv3_filenameginput,delp=g_delp,tsen=g_tsen,w=g_w,iope=iope) + call gsi_fv3ncdf_read_ens_parallel_over_ens(fv3_filenameginput%tracers,fv3_filenameginput,q=g_q,oz=g_oz,ql=g_ql,qr=g_qr,& + qs=g_qs,qi=g_qi,qg=g_qg,iope=iope) + call gsi_fv3ncdf_read_ens_parallel_over_ens(fv3_filenameginput%phyvars,fv3_filenameginput,dbz=g_dbz,iope=iope) + else + call gsi_fv3ncdf_read_ens_parallel_over_ens(fv3_filenameginput%dynvars,fv3_filenameginput,delp=g_delp,tsen=g_tsen,iope=iope) + call gsi_fv3ncdf_read_ens_parallel_over_ens(fv3_filenameginput%tracers,fv3_filenameginput,q=g_q,oz=g_oz,iope=iope) + end if + else + write(6,*) "Warning: we can only grab fields from restart files not cold start files for ensemble!" + endif + + + if (fv3sar_ensemble_opt == 0) then + g_prsi(:,:,grd_ens%nsig+1)=eta1_ll(grd_ens%nsig+1) !thinkto be done , should use eta1_ll from ensemble grid + do i=grd_ens%nsig,1,-1 + g_prsi(:,:,i)=g_delp(:,:,i)*0.001_r_kind+g_prsi(:,:,i+1) + enddo + g_ps(:,:)=g_prsi(:,:,1) + + endif + + !! tsen2tv !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + do k=1,grd_ens%nsig + do j=1,grd_ens%nlon + do i=1,grd_ens%nlat + g_tv(i,j,k)=g_tsen(i,j,k)*(one+fv*g_q(i,j,k)) + enddo + enddo + enddo + if (.not.q_hyb_ens) then + ice=.true. + iderivative=0 + do k=1,grd_ens%nsig + kp=k+1 + do j=1,grd_ens%nlon + do i=1,grd_ens%nlat + g_prsl(i,j,k)=(g_prsi(i,j,k)+g_prsi(i,j,kp))*half + end do + end do + end do + call genqsat(g_rh,g_tsen(1,1,1),g_prsl(1,1,1),grd_ens%nlat,grd_ens%nlon,grd_ens%nsig,ice,iderivative) + do k=1,grd_ens%nsig + do j=1,grd_ens%nlon + do i=1,grd_ens%nlat + g_rh(i,j,k) = g_q(i,j,k)/g_rh(i,j,k) + end do + end do + end do + else + do k=1,grd_ens%nsig + do j=1,grd_ens%nlon + do i=1,grd_ens%nlat + g_rh(i,j,k) = g_q(i,j,k) + end do + end do + end do + end if + end if ! mype + + return + end subroutine general_read_fv3_regional_parallel_over_ens + + + subroutine parallel_read_fv3_step2(this,mype,iope, & + g_ps,g_u,g_v,g_tv,g_rh,g_ql,g_oz,g_w,g_qr,g_qs,g_qi,& + g_qg,g_dbz, & + gg_ps,gg_tv,gg_u,gg_v,gg_rh,gg_w,gg_dbz,gg_qr,& + gg_qs,gg_qi,gg_qg,gg_ql) + + !$$$ subprogram documentation block + ! . + ! subprogram: parallel_read_fv3_step2 distribute all fields into all tasks + ! prgmmr: Y. Wang and X. Wang org: OU/MAP date: 2022-04-01 + ! + ! abstract: All fields have been read in by general_read_fv3_regional_parallel_over_ens. + ! Different tasks contain the data from different members. + ! This program will divided the full-domain fields into subdomains + ! and assign them to all tasks. poc: xuguang.wang@ou.edu + ! + ! program history log: + ! + ! 2022-04-01 Y. Wang and X. Wang - Changed from the code for WRF-ARW + ! + ! input argument list: + ! + ! output argument list: + ! + ! attributes: + ! language: f90 + ! machine: ibm RS/6000 SP + ! + !$$$ end documentation block + + use hybrid_ensemble_parameters, only: grd_ens + use mpimod, only: mpi_comm_world,ierror,mpi_rtype + use kinds, only: r_kind,r_single,i_kind + use gridmod,only: itotsub + + implicit none + + ! + ! Declare passed variables + class(get_fv3_regional_ensperts_class), intent(inout) :: this + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),intent(out):: & + g_u,g_v,g_tv,g_rh,g_ql,g_oz + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),intent(out),optional::& + g_w,g_qr,g_qs,g_qi,g_qg,g_dbz + integer(i_kind), intent(in) :: mype, iope + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2),intent(out):: g_ps + + ! The gg_ arrays are only sent by the rank doing I/O (mype==iope) + real(r_kind),optional,dimension(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig) :: & + gg_u,gg_v,gg_tv,gg_rh + + real(r_kind),optional,dimension(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig) :: & + gg_w,gg_dbz,gg_qr,gg_qs,gg_qi,gg_qg,gg_ql + real(r_kind),optional,dimension(grd_ens%nlat,grd_ens%nlon):: gg_ps + + ! Declare local variables + real(r_kind),allocatable,dimension(:):: wrk_send_2d + integer(i_kind) :: k + + ! transfer data from root to subdomains on each task + ! scatterv used, since full grids exist only on root task. + allocate(wrk_send_2d(grd_ens%itotsub)) + ! first PS (output from fill_regional_2d is a column vector with a halo) + if(mype==iope) call this%fill_regional_2d(gg_ps,wrk_send_2d) + call mpi_scatterv(wrk_send_2d,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype, & + g_ps,grd_ens%ijn_s(mype+1),mpi_rtype,iope,mpi_comm_world,ierror) + ! then TV,U,V,RH + do k=1,grd_ens%nsig + if (mype==iope) then + call this%fill_regional_2d(gg_tv(:,:,k),wrk_send_2d) + endif + call mpi_scatterv(wrk_send_2d,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype, & + g_tv(1,1,k),grd_ens%ijn_s(mype+1),mpi_rtype,iope,mpi_comm_world,ierror) + if (mype==iope) call this%fill_regional_2d(gg_u(1,1,k),wrk_send_2d) + call mpi_scatterv(wrk_send_2d,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype, & + g_u(1,1,k),grd_ens%ijn_s(mype+1),mpi_rtype,iope,mpi_comm_world,ierror) + if (mype==iope) call this%fill_regional_2d(gg_v(1,1,k),wrk_send_2d) + call mpi_scatterv(wrk_send_2d,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype, & + g_v(1,1,k),grd_ens%ijn_s(mype+1),mpi_rtype,iope,mpi_comm_world,ierror) + if (mype==iope) call this%fill_regional_2d(gg_rh(1,1,k),wrk_send_2d) + call mpi_scatterv(wrk_send_2d,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype, & + g_rh(1,1,k),grd_ens%ijn_s(mype+1),mpi_rtype,iope,mpi_comm_world,ierror) + if( present(g_dbz) )then + if (mype==iope) call this%fill_regional_2d(gg_w(1,1,k),wrk_send_2d) + call mpi_scatterv(wrk_send_2d,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype, & + g_w(1,1,k),grd_ens%ijn_s(mype+1),mpi_rtype,iope,mpi_comm_world,ierror) + if (mype==iope) call this%fill_regional_2d(gg_dbz(1,1,k),wrk_send_2d) + call mpi_scatterv(wrk_send_2d,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype,& + g_dbz(1,1,k),grd_ens%ijn_s(mype+1),mpi_rtype,iope,mpi_comm_world,ierror) + if (mype==iope) call this%fill_regional_2d(gg_qr(1,1,k),wrk_send_2d) + call mpi_scatterv(wrk_send_2d,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype,& + g_qr(1,1,k),grd_ens%ijn_s(mype+1),mpi_rtype,iope,mpi_comm_world,ierror) + if (mype==iope) call this%fill_regional_2d(gg_qs(1,1,k),wrk_send_2d) + call mpi_scatterv(wrk_send_2d,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype,& + g_qs(1,1,k),grd_ens%ijn_s(mype+1),mpi_rtype,iope,mpi_comm_world,ierror) + if (mype==iope) call this%fill_regional_2d(gg_qi(1,1,k),wrk_send_2d) + call mpi_scatterv(wrk_send_2d,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype,& + g_qi(1,1,k),grd_ens%ijn_s(mype+1),mpi_rtype,iope,mpi_comm_world,ierror) + if (mype==iope) call this%fill_regional_2d(gg_qg(1,1,k),wrk_send_2d) + call mpi_scatterv(wrk_send_2d,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype, & + g_qg(1,1,k),grd_ens%ijn_s(mype+1),mpi_rtype,iope,mpi_comm_world,ierror) + if (mype==iope) call this%fill_regional_2d(gg_ql(1,1,k),wrk_send_2d) + call mpi_scatterv(wrk_send_2d,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype,& + g_ql(1,1,k),grd_ens%ijn_s(mype+1),mpi_rtype,iope,mpi_comm_world,ierror) + end if + enddo + deallocate(wrk_send_2d) + end subroutine parallel_read_fv3_step2 + + subroutine fill_regional_2d(fld_in,fld_out) + !$$$ subprogram documentation block + ! . . . . + ! subprogram: fill_regional_2d + ! prgmmr: mizzi org: ncar/mmm date: 2010-08-11 + ! + ! abstract: create a column vector for the subdomain (including halo) + ! from global 2d grid. + ! + ! + ! program history log: + ! 2010-08-11 parrish, initial documentation + ! 2012-03-12 whitaker, remove nx,ny,itotsub from argument list. + ! + ! input argument list: + ! + ! output argument list: + ! + ! attributes: + ! language: f90 + ! machine: ibm RS/6000 SP + ! + !$$$ end documentation block + use kinds, only: r_kind,i_kind + use hybrid_ensemble_parameters, only: grd_ens + implicit none + real(r_kind),dimension(grd_ens%nlat,grd_ens%nlon)::fld_in + real(r_kind),dimension(grd_ens%itotsub)::fld_out + integer(i_kind):: i,j,k + do k=1,grd_ens%itotsub + i=grd_ens%ltosj_s(k) + j=grd_ens%ltosi_s(k) + fld_out(k)=fld_in(j,i) + enddo + return + end subroutine fill_regional_2d + subroutine ens_spread_dualres_regional_fv3_regional(this,mype,en_perts,nelen) !$$$ subprogram documentation block ! . . . . diff --git a/src/gsi/general_specmod.f90 b/src/gsi/general_specmod.f90 index 439e26e431..20feae98de 100644 --- a/src/gsi/general_specmod.f90 +++ b/src/gsi/general_specmod.f90 @@ -64,6 +64,7 @@ module general_specmod ! set subroutines to public public :: general_init_spec_vars public :: general_destroy_spec_vars + public :: general_spec_multwgt ! set passed variables to public public :: spec_vars public :: spec_cut @@ -306,6 +307,36 @@ subroutine general_init_spec_vars(sp,jcap,jcap_test,nlat_a,nlon_a,eqspace) return end subroutine general_init_spec_vars + subroutine general_spec_multwgt(sp,spcwrk,spcwgt) +! 2017-03-30 J. Kay, X. Wang - add general_spec_multwgt for scale-dependent +! weighting of mixed resolution ensemble, +! POC: xuguang.wang@ou.edu +! + implicit none + type(spec_vars),intent(in) :: sp + real(r_kind),dimension(sp%nc),intent(inout) :: spcwrk + real(r_kind),dimension(0:sp%jcap),intent(in) :: spcwgt + + integer(i_kind) ii1,l,m,jmax,ks,n + +!! Code borrowed from spvar in splib + jmax=sp%jcap + + do n=0,jmax + ks=2*n + spcwrk(ks+1)=spcwrk(ks+1)*spcwgt(n) + end do + do n=0,jmax + do l=MAX(1,n-jmax),MIN(n,jmax) + ks=l*(2*jmax+(-1)*(l-1))+2*n + spcwrk(ks+1) = spcwrk(ks+1)*spcwgt(n) + spcwrk(ks+2) = spcwrk(ks+2)*spcwgt(n) + end do + end do + + return + end subroutine general_spec_multwgt + subroutine general_destroy_spec_vars(sp) !$$$ subprogram documentation block ! . . . . diff --git a/src/gsi/get_gefs_ensperts_dualres.f90 b/src/gsi/get_gefs_ensperts_dualres.f90 index fa3d0ecbdd..1cb7586a89 100644 --- a/src/gsi/get_gefs_ensperts_dualres.f90 +++ b/src/gsi/get_gefs_ensperts_dualres.f90 @@ -68,6 +68,7 @@ subroutine get_gefs_ensperts_dualres use gsi_enscouplermod, only: gsi_enscoupler_create_sub2grid_info use gsi_enscouplermod, only: gsi_enscoupler_destroy_sub2grid_info use general_sub2grid_mod, only: sub2grid_info,general_sub2grid_create_info,general_sub2grid_destroy_info + use hybrid_ensemble_parameters, only: nsclgrp,sp_ens,spc_multwgt,global_spectral_filter_sd implicit none real(r_kind),pointer,dimension(:,:) :: ps @@ -78,6 +79,8 @@ subroutine get_gefs_ensperts_dualres real(r_kind),pointer,dimension(:,:):: x2 type(gsi_bundle),allocatable,dimension(:) :: en_read type(gsi_bundle):: en_bar + type(gsi_bundle) :: en_pertstmp1 + type(gsi_bundle) :: en_pertstmp2 ! type(gsi_grid) :: grid_ens real(r_kind) bar_norm,sig_norm,kapr,kap1 ! real(r_kind),allocatable,dimension(:,:):: z,sst2 @@ -91,6 +94,7 @@ subroutine get_gefs_ensperts_dualres ! integer(i_kind) il,jl logical ice,hydrometeor type(sub2grid_info) :: grd_tmp + integer(i_kind) :: ig0,ig ! Create perturbations grid and get variable names from perturbations if(en_perts(1,1,1)%grid%im/=grd_ens%lat2.or. & @@ -130,6 +134,16 @@ subroutine get_gefs_ensperts_dualres if ( istatus /= 0 ) & call die('get_gefs_ensperts_dualres',': trouble creating en_bar bundle, istatus =',istatus) +! Allocate bundle used for temporary usage + if( nsclgrp > 1 .and. global_spectral_filter_sd )then + call gsi_bundlecreate(en_pertstmp1,en_perts(1,1,1)%grid,'aux-ens-read',istatus,names2d=cvars2d,names3d=cvars3d) + call gsi_bundlecreate(en_pertstmp2,en_perts(1,1,1)%grid,'aux-ens-read',istatus,names2d=cvars2d,names3d=cvars3d) + if(istatus/=0) then + write(6,*)' get_gefs_ensperts_dualres: trouble creating en_read like tempbundle' + call stop2(999) + endif + end if + ! Allocate bundle used for reading members allocate(en_read(n_ens)) do n=1,n_ens @@ -318,12 +332,30 @@ subroutine get_gefs_ensperts_dualres end do end do ntlevs_ens_loop !end do over bins + if(nsclgrp > 1 .and. global_spectral_filter_sd) then + do m=1,ntlevs_ens + do n=1,n_ens + en_pertstmp1%values=en_perts(n,1,m)%valuesr4 + do ig=1,nsclgrp + call apply_scaledepwgts(grd_ens,sp_ens,en_pertstmp1,spc_multwgt(:,ig),en_pertstmp2) + en_perts(n,ig,m)%valuesr4=en_pertstmp2%values + enddo + enddo + enddo + endif + do n=n_ens,1,-1 call gsi_bundledestroy(en_read(n),istatus) if ( istatus /= 0 ) & call die('get_gefs_ensperts_dualres',': trouble destroying en_read bundle, istatus = ', istatus) end do deallocate(en_read) + if(nsclgrp > 1 .and. global_spectral_filter_sd) then + call gsi_bundledestroy(en_pertstmp1,istatus) + call gsi_bundledestroy(en_pertstmp2,istatus) + if ( istatus /= 0 ) & + call die('get_gefs_ensperts_dualres',': trouble destroying en_pertstmp2 bundle, istatus = ', istatus) + end if call gsi_enscoupler_destroy_sub2grid_info(grd_tmp) ! mm1=mype+1 diff --git a/src/gsi/gsi_files.cmake b/src/gsi/gsi_files.cmake index 461b49ddf6..1658c83bf4 100644 --- a/src/gsi/gsi_files.cmake +++ b/src/gsi/gsi_files.cmake @@ -87,6 +87,7 @@ anisofilter_glb.f90 antcorr_application.f90 antest_maps0.f90 antest_maps0_glb.f90 +apply_scaledepwgts.f90 atms_spatial_average_mod.f90 balmod.f90 berror.f90 diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90 index eb7a86160f..2edde4723f 100644 --- a/src/gsi/gsi_rfv3io_mod.f90 +++ b/src/gsi/gsi_rfv3io_mod.f90 @@ -68,6 +68,7 @@ module gsi_rfv3io_mod character(len=:),allocatable :: ak_bk !='fv3_akbk' character(len=:),allocatable :: dynvars !='fv3_dynvars' character(len=:),allocatable :: tracers !='fv3_tracer' + character(len=:),allocatable :: phyvars !='fv3_phyvars' character(len=:),allocatable :: sfcdata !='fv3_sfcdata' character(len=:),allocatable :: couplerres!='coupler.res' contains @@ -91,8 +92,9 @@ module gsi_rfv3io_mod type(sub2grid_info) :: grd_fv3lam_tracer_ionouv type(sub2grid_info) :: grd_fv3lam_tracerchem_ionouv type(sub2grid_info) :: grd_fv3lam_tracersmoke_ionouv + type(sub2grid_info) :: grd_fv3lam_phyvar_ionouv type(sub2grid_info) :: grd_fv3lam_uv - integer(i_kind) ,parameter:: ndynvarslist=13, ntracerslist=8 + integer(i_kind) ,parameter:: ndynvarslist=13, ntracerslist=8, nphyvarslist=1 character(len=max_varname_length), dimension(ndynvarslist), parameter :: & vardynvars = [character(len=max_varname_length) :: & @@ -100,12 +102,14 @@ module gsi_rfv3io_mod character(len=max_varname_length), dimension(ntracerslist+naero_cmaq_fv3+7+naero_smoke_fv3), parameter :: & vartracers = [character(len=max_varname_length) :: & 'q','oz','ql','qi','qr','qs','qg','qnr',aeronames_cmaq_fv3,'pm25at','pm25ac','pm25co','pm2_5','amassi','amassj','amassk',aeronames_smoke_fv3] - character(len=max_varname_length), dimension(15+naero_cmaq_fv3+7+naero_smoke_fv3), parameter :: & + character(len=max_varname_length), dimension(nphyvarslist), parameter :: & + varphyvars = [character(len=max_varname_length) :: 'dbz'] + character(len=max_varname_length), dimension(16+naero_cmaq_fv3+7+naero_smoke_fv3), parameter :: & varfv3name = [character(len=max_varname_length) :: & - 'u','v','W','T','delp','sphum','o3mr','liq_wat','ice_wat','rainwat','snowwat','graupel','rain_nc','ps','DZ', & + 'u','v','W','T','delp','sphum','o3mr','liq_wat','ice_wat','rainwat','snowwat','graupel','rain_nc','ref_f3d','ps','DZ', & aeronames_cmaq_fv3,'pm25at','pm25ac','pm25co','pm2_5','amassi','amassj','amassk',aeronames_smoke_fv3], & vgsiname = [character(len=max_varname_length) :: & - 'u','v','w','tsen','delp','q','oz','ql','qi','qr','qs','qg','qnr','ps','delzinc', & + 'u','v','w','tsen','delp','q','oz','ql','qi','qr','qs','qg','qnr','dbz','ps','delzinc', & aeronames_cmaq_fv3,'pm25at','pm25ac','pm25co','pm2_5','amassi','amassj','amassk',aeronames_smoke_fv3] character(len=max_varname_length),dimension(:),allocatable:: name_metvars2d character(len=max_varname_length),dimension(:),allocatable:: name_metvars3d @@ -119,6 +123,8 @@ module gsi_rfv3io_mod public :: gsi_fv3ncdf_read_v1 public :: gsi_fv3ncdf_readuv public :: gsi_fv3ncdf_readuv_v1 + public :: gsi_fv3ncdf_read_ens_parallel_over_ens + public :: gsi_fv3ncdf_readuv_ens_parallel_over_ens public :: read_fv3_files public :: read_fv3_netcdf_guess public :: wrfv3_netcdf @@ -131,6 +137,7 @@ module gsi_rfv3io_mod public :: ijns,ijns2d,displss,displss2d,ijnz,displsz_g public :: fv3lam_io_dynmetvars3d_nouv,fv3lam_io_tracermetvars3d_nouv public :: fv3lam_io_tracerchemvars3d_nouv,fv3lam_io_tracersmokevars3d_nouv + public :: fv3lam_io_phymetvars3d_nouv public :: fv3lam_io_dynmetvars2d_nouv,fv3lam_io_tracermetvars2d_nouv integer(i_kind) mype_u,mype_v,mype_t,mype_q,mype_p,mype_delz,mype_oz,mype_ql @@ -158,6 +165,7 @@ module gsi_rfv3io_mod ! copy of cvars3d excluding uv 3-d fields character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_tracermetvars3d_nouv ! copy of cvars3d excluding uv 3-d fields + character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_phymetvars3d_nouv character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_tracerchemvars3d_nouv character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_tracersmokevars3d_nouv ! copy of cvars3d excluding uv 3-d fields @@ -169,8 +177,10 @@ module gsi_rfv3io_mod !to define names in gsibundle character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_names_gsibundle_tracer_nouv !to define names in gsibundle + character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_names_gsibundle_phyvar_nouv type(gsi_bundle):: gsibundle_fv3lam_dynvar_nouv type(gsi_bundle):: gsibundle_fv3lam_tracer_nouv + type(gsi_bundle):: gsibundle_fv3lam_phyvar_nouv type(gsi_bundle):: gsibundle_fv3lam_tracerchem_nouv type(gsi_bundle):: gsibundle_fv3lam_tracersmoke_nouv @@ -204,6 +214,12 @@ subroutine fv3regfilename_init(this,it) write(filename,"(A11,I2.2)") 'fv3_tracer_',ifilesig(it) this%tracers=trim(filename) endif + if (it == ntguessig) then + this%phyvars='fv3_phyvars' + else + write(filename,"(A12,I2.2)") 'fv3_phyvars_',ifilesig(it) + this%phyvars=trim(filename) + endif if (it == ntguessig) then this%sfcdata='fv3_sfcdata' else @@ -740,6 +756,9 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) ! ! abstract: read guess for FV3 regional model ! program history log: +! 2022-04-01 Y. Wang and X. Wang - add capability to read reflectivity +! for direct radar EnVar DA using reflectivity as state +! variable, poc: xuguang.wang@ou.edu ! attributes: ! language: f90 ! machine: ibm RS/6000 SP @@ -768,6 +787,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) use gsi_metguess_mod, only: gsi_metguess_get use netcdf, only:nf90_open,nf90_close,nf90_inquire,nf90_nowrite, nf90_format_netcdf4 use gsi_chemguess_mod, only: gsi_chemguess_get + use obsmod, only: if_model_dbz implicit none @@ -797,7 +817,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) real(r_kind),dimension(:,:,:),pointer::ges_qg=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_qnr=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_w=>NULL() - + real(r_kind),dimension(:,:,:),pointer::ges_dbz=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_aalj=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_acaj=>NULL() @@ -890,8 +910,8 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) integer(i_kind),dimension(:,:),allocatable:: lnames integer(i_kind),dimension(:,:),allocatable:: uvlnames integer(i_kind):: inner_vars,numfields - integer(i_kind):: ndynvario2d,ntracerio2d,ilev,jdynvar,jtracer - integer(i_kind):: iuv,ndynvario3d,ntracerio3d + integer(i_kind):: ndynvario2d,ntracerio2d,ilev,jdynvar,jtracer,jphyvar + integer(i_kind):: iuv,ndynvario3d,ntracerio3d,nphyvario3d integer(i_kind):: ntracerchemio3d,ntracersmokeio3d integer(i_kind):: loc_id,ncfmt @@ -957,6 +977,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) iuv=0 ndynvario3d=0 ntracerio3d=0 + nphyvario3d=0 do i=1,size(name_metvars3d) vartem=trim(name_metvars3d(i)) if(trim(vartem)=='u'.or.trim(vartem)=='v') then @@ -967,6 +988,8 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) ndynvario3d=ndynvario3d+1 else if (ifindstrloc(vartracers,trim(vartem))> 0) then ntracerio3d=ntracerio3d+1 + else if (ifindstrloc(varphyvars,trim(vartem))> 0) then + nphyvario3d=nphyvario3d+1 else write(6,*)'the metvarname1 ',trim(vartem),' has not been considered yet, stop' call stop2(333) @@ -978,6 +1001,12 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) write(6,*)"the set up for met variable is not as expected, abort" call stop2(222) endif + if ( if_model_dbz ) then + if( nphyvario3d<=0 ) then + write(6,*)"the set up for met variable (phyvar) is not as expected, abort" + call stop2(223) + end if + endif if (fv3sar_bg_opt == 0.and.ifindstrloc(name_metvars3d,'delp') <= 0)then ndynvario3d=ndynvario3d+1 ! for delp endif @@ -987,6 +1016,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) if (l_reg_update_hydro_delz.and.fv3sar_bg_opt==0) ndynvario3d=ndynvario3d+1 ! for delzinc allocate(fv3lam_io_dynmetvars3d_nouv(ndynvario3d)) allocate(fv3lam_io_tracermetvars3d_nouv(ntracerio3d)) + allocate(fv3lam_io_phymetvars3d_nouv(nphyvario3d)) if (laeroana_fv3cmaq) then allocate(fv3lam_io_tracerchemvars3d_nouv(naero_cmaq_fv3+7)) @@ -998,6 +1028,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) jdynvar=0 jtracer=0 + jphyvar=0 do i=1,size(name_metvars3d) vartem=trim(name_metvars3d(i)) if(.not.(trim(vartem)=='u'.or.trim(vartem)=='v'.or.trim(vartem)=='iqr')) then @@ -1012,6 +1043,9 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) else if (ifindstrloc(vartracers,trim(vartem)) > 0) then jtracer=jtracer+1 fv3lam_io_tracermetvars3d_nouv(jtracer)=trim(vartem) + else if (ifindstrloc(varphyvars,trim(vartem)) > 0) then + jphyvar=jphyvar+1 + fv3lam_io_phymetvars3d_nouv(jphyvar)=trim(vartem) else write(6,*)'the metvarname ',vartem,' is not expected, stop' call flush(6) @@ -1027,7 +1061,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) jdynvar=jdynvar+1 fv3lam_io_dynmetvars3d_nouv(jdynvar)="delzinc" endif - if(jdynvar /= ndynvario3d.or.jtracer /= ntracerio3d ) then + if(jdynvar /= ndynvario3d.or.jtracer /= ntracerio3d.or.jphyvar /= nphyvario3d ) then write(6,*)'ndynvario3d is not as expected, stop' call flush(6) call stop2(333) @@ -1035,6 +1069,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) if(mype == 0) then write(6,*) ' fv3lam_io_dynmetvars3d_nouv is ',(trim(fv3lam_io_dynmetvars3d_nouv(i)),i=1,ndynvario3d) write(6,*) ' fv3lam_io_tracermevars3d_nouv is ',(trim(fv3lam_io_tracermetvars3d_nouv(i)),i=1,ntracerio3d) + write(6,*) ' fv3lam_io_phymetvars3d_nouv is ',(trim(fv3lam_io_phymetvars3d_nouv(i)),i=1,nphyvario3d) endif ndynvario2d=0 @@ -1047,7 +1082,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) else if (ifindstrloc(vartracers,trim(vartem)) > 0) then ntracerio2d=ntracerio2d+1 else if(trim(vartem)=='z') then - write(6,*)'the metvarname ',trim(vartem),' will be dealt separately' + if(mype == 0) write(6,*)'the metvarname ',trim(vartem),' will be dealt separately' else if(trim(vartem)=='t2m') then else if(trim(vartem)=='q2m') then else @@ -1165,6 +1200,11 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) ntracerio2d=0 endif + if( if_model_dbz )then + call gsi_bundlecreate(gsibundle_fv3lam_phyvar_nouv,GSI_MetGuess_Bundle(it)%grid,'gsibundle_fv3lam_phyvar_nouv',istatus, & + names3d=fv3lam_io_phymetvars3d_nouv) + end if + if (laeroana_fv3cmaq) then if (allocated(fv3lam_io_tracerchemvars3d_nouv) ) then call gsi_bundlecreate(gsibundle_fv3lam_tracerchem_nouv,GSI_ChemGuess_Bundle(it)%grid,'gsibundle_fv3lam_tracerchem_nouv',istatus, & @@ -1254,6 +1294,22 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) endif + if ( if_model_dbz )then + inner_vars=1 + numfields=inner_vars*(nphyvario3d*grd_a%nsig) + deallocate(lnames,names) + allocate(lnames(1,numfields),names(1,numfields)) + ilev=1 + do i=1,nphyvario3d + do k=1,grd_a%nsig + lnames(1,ilev)=k + names(1,ilev)=trim(fv3lam_io_phymetvars3d_nouv(i)) + ilev=ilev+1 + enddo + enddo + call general_sub2grid_create_info(grd_fv3lam_phyvar_ionouv,inner_vars,grd_a%nlat,& + grd_a%nlon,grd_a%nsig,numfields,regional,names=names,lnames=lnames) + end if inner_vars=2 numfields=grd_a%nsig @@ -1279,15 +1335,19 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'tv' ,ges_tv ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'q' ,ges_q ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'oz' ,ges_oz ,istatus );ier=ier+istatus - if (l_use_dbz_directDA) then + if (l_use_dbz_directDA .or. if_model_dbz) then call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'ql' ,ges_ql ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qi' ,ges_qi ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qr' ,ges_qr ,istatus );ier=ier+istatus - call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'iqr' ,ges_iqr ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qs' ,ges_qs ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qg' ,ges_qg ,istatus );ier=ier+istatus - call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qnr',ges_qnr ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'w' , ges_w ,istatus );ier=ier+istatus + if (l_use_dbz_directDA) then + call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'iqr' ,ges_iqr ,istatus );ier=ier+istatus + call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qnr',ges_qnr ,istatus );ier=ier+istatus + end if + if(if_model_dbz) & + call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'dbz' , ges_dbz ,istatus );ier=ier+istatus end if if (ier/=0) call die(trim(myname),'cannot get pointers for fv3 met-fields, ier =',ier) @@ -1343,6 +1403,10 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) & ,fv3filenamegin(it)%dynvars,fv3filenamegin(it)) call gsi_fv3ncdf_read(grd_fv3lam_tracer_ionouv,gsibundle_fv3lam_tracer_nouv & & ,fv3filenamegin(it)%tracers,fv3filenamegin(it)) + if( if_model_dbz )then + call gsi_fv3ncdf_read(grd_fv3lam_phyvar_ionouv,gsibundle_fv3lam_phyvar_nouv & + & ,fv3filenamegin(it)%phyvars,fv3filenamegin(it)) + end if if (laeroana_fv3cmaq) then call gsi_fv3ncdf_read(grd_fv3lam_tracerchem_ionouv,gsibundle_fv3lam_tracerchem_nouv & & ,fv3filenamegin(it)%tracers,fv3filenamegin(it)) @@ -1439,6 +1503,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) call gsi_copy_bundle(gsibundle_fv3lam_tracersmoke_nouv,GSI_ChemGuess_Bundle(it)) endif + if(if_model_dbz) call gsi_copy_bundle(gsibundle_fv3lam_phyvar_nouv,GSI_MetGuess_Bundle(it)) call GSI_BundleGetPointer ( gsibundle_fv3lam_dynvar_nouv, 'tsen' ,ges_tsen_readin ,istatus );ier=ier+istatus !! tsen2tv !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do k=1,nsig @@ -2151,13 +2216,15 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) character(len=max_varname_length) :: varname,vgsiname character(len=max_varname_length) :: name character(len=max_varname_length) :: filenamein2 - + real(r_kind),allocatable,dimension(:,:):: uu2d_tmp + integer(i_kind) :: countloc_tmp(3),startloc_tmp(3) integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(3),startloc(3) integer(i_kind) ilev,ilevtot,inative - integer(i_kind) kbgn,kend + integer(i_kind) kbgn,kend,len + logical :: phy_smaller_domain integer(i_kind) gfile_loc,iret,var_id - integer(i_kind) nz,nzp1,mm1 + integer(i_kind) nz,nzp1,mm1,nx_phy ! for io_layout > 1 real(r_kind),allocatable,dimension(:,:):: uu2d_layout integer(i_kind) :: nio @@ -2212,6 +2279,22 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) inative=nzp1-ilev startloc=(/1,1,inative/) countloc=(/nxcase,nycase,1/) + ! Variable ref_f3d in phy_data.nc has a smaller domain size than + ! dynvariables and tracers as well as a reversed order in vertical + if ( trim(adjustl(varname)) == 'ref_f3d' )then + iret=nf90_inquire_dimension(gfile_loc,1,name,len) + if(trim(name)=='xaxis_1') nx_phy=len + if( nx_phy == nxcase )then + allocate(uu2d_tmp(nxcase,nycase)) + countloc_tmp=(/nxcase,nycase,1/) + phy_smaller_domain = .false. + else + allocate(uu2d_tmp(nxcase-6,nycase-6)) + countloc_tmp=(/nxcase-6,nycase-6,1/) + phy_smaller_domain = .true. + end if + startloc_tmp=(/1,1,ilev/) + end if if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 @@ -2224,7 +2307,22 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) enddo else iret=nf90_inq_varid(gfile_loc,trim(adjustl(varname)),var_id) - iret=nf90_get_var(gfile_loc,var_id,uu2d,start=startloc,count=countloc) + if ( trim(adjustl(varname)) == 'ref_f3d' )then + uu2d = 0.0_r_kind + iret=nf90_get_var(gfile_loc,var_id,uu2d_tmp,start=startloc_tmp,count=countloc_tmp) + where(uu2d_tmp < 0.0_r_kind) + uu2d_tmp = 0.0_r_kind + endwhere + + if( phy_smaller_domain )then + uu2d(4:nxcase-3,4:nycase-3) = uu2d_tmp + else + uu2d(1:nxcase,1:nycase) = uu2d_tmp + end if + deallocate(uu2d_tmp) + else + iret=nf90_get_var(gfile_loc,var_id,uu2d,start=startloc,count=countloc) + end if endif call fv3_h_to_ll(uu2d,hwork(1,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,grid_reverse_flag) @@ -2642,6 +2740,395 @@ subroutine gsi_fv3ncdf_readuv_v1(grd_uv,ges_u,ges_v,fv3filenamegin) end subroutine gsi_fv3ncdf_readuv_v1 +subroutine gsi_fv3ncdf_read_ens_parallel_over_ens(filenamein,fv3filenamegin, & + delp,tsen,w,q,oz,ql,qr,qs,qi,qg,dbz,iope) +!$$$ subprogram documentation block +! . . . . +! subprogram: gsi_fv3ncdf_read_ens_parallel_over_ens +! program history log: +! 2022-04-01 Y. Wang and X. Wang, changed from gsi_fv3ncdf_read_ens +! for FV3LAM ensemble parallel IO in hybrid EnVar +! poc: xuguang.wang@ou.edu +! +! abstract: read in fields excluding u and v +! program history log: +! +! input argument list: +! filenamein - file name to read from +! iope - pe to read in the field +! +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ end documentation block + + + use kinds, only: r_kind,i_kind + use mpimod, only: mpi_comm_world,mpi_rtype,mype + use mpimod, only: MPI_INFO_NULL + use netcdf, only: nf90_open,nf90_close,nf90_get_var,nf90_noerr + use netcdf, only: nf90_nowrite,nf90_inquire,nf90_inquire_dimension + use netcdf, only: nf90_inquire_variable + use netcdf, only: nf90_inq_varid + use gridmod, only: nsig,nlon,nlat + use mod_fv3_lola, only: fv3_h_to_ll + use gsi_bundlemod, only: gsi_bundle + use general_sub2grid_mod, only: sub2grid_info,general_grid2sub + + implicit none + character(*),intent(in):: filenamein + type (type_fv3regfilenameg),intent(in) ::fv3filenamegin + integer(i_kind) ,intent(in ) :: iope + real(r_kind),allocatable,dimension(:,:):: uu2d, uu2d_tmp + real(r_kind),allocatable,dimension(:):: wrk_send_2d + real(r_kind),dimension(nlat,nlon,nsig):: hwork + real(r_kind),dimension(nlat,nlon,nsig),intent(out),optional:: delp,tsen,w,q,oz,ql,qr,qs,qi,qg,dbz + character(len=max_varname_length) :: varname,vgsiname + character(len=max_varname_length) :: name + character(len=max_varname_length), allocatable,dimension(:) :: varname_files + + integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(3),startloc(3),countloc_tmp(3),startloc_tmp(3) + integer(i_kind) ilev,ilevtot,inative,ivar + integer(i_kind) kbgn,kend + integer(i_kind) gfile_loc,iret,var_id + integer(i_kind) nz,nzp1,mm1,len,nx_phy + logical :: phy_smaller_domain +! for io_layout > 1 + real(r_kind),allocatable,dimension(:,:):: uu2d_layout + integer(i_kind) :: nio + integer(i_kind),allocatable :: gfile_loc_layout(:) + character(len=180) :: filename_layout + + mm1=mype+1 + nloncase=nlon + nlatcase=nlat + nxcase=nx + nycase=ny + kbgn=1 + kend=nsig + + if( mype == iope )then + allocate(uu2d(nxcase,nycase)) + if( present(delp).or.present(tsen).or.present(w) )then ! dynvars + if( present(w) )then + allocate(varname_files(3)) + varname_files = (/'T ','delp','W '/) + else + allocate(varname_files(2)) + varname_files = (/'T ','delp'/) + end if + end if + if( present(q).or.present(ql).or.present(qr) )then ! tracers + if(present(qr))then + allocate(varname_files(7)) + varname_files = (/'sphum ','o3mr ','liq_wat','ice_wat','rainwat','snowwat','graupel'/) + else + allocate(varname_files(2)) + varname_files = (/'sphum',' o3mr'/) + end if + end if + if( present(dbz) )then ! phyvars: dbz + allocate(varname_files(1)) + varname_files = (/'ref_f3d'/) + end if + + if(fv3_io_layout_y > 1) then + allocate(gfile_loc_layout(0:fv3_io_layout_y-1)) + do nio=0,fv3_io_layout_y-1 + write(filename_layout,'(a,a,I4.4)') trim(filenamein),'.',nio + iret=nf90_open(filename_layout,nf90_nowrite,gfile_loc_layout(nio),comm=mpi_comm_world,info=MPI_INFO_NULL) + if(iret/=nf90_noerr) then + write(6,*)' gsi_fv3ncdf_read: problem opening ',trim(filename_layout),gfile_loc_layout(nio),', Status = ',iret + call flush(6) + call stop2(333) + endif + enddo + else + iret=nf90_open(filenamein,nf90_nowrite,gfile_loc) + if(iret/=nf90_noerr) then + write(6,*)' gsi_fv3ncdf_read: problem opening ',trim(filenamein),gfile_loc,', Status = ',iret + call flush(6) + call stop2(333) + endif + endif + do ivar = 1, size(varname_files) + do ilevtot=kbgn,kend + ilev=ilevtot + nz=nsig + nzp1=nz+1 + inative=nzp1-ilev + startloc=(/1,1,inative/) + countloc=(/nxcase,nycase,1/) + varname = trim(varname_files(ivar)) + ! Variable ref_f3d in phy_data.nc has a smaller domain size than + ! dynvariables and tracers as well as a reversed order in vertical + if ( trim(adjustl(varname)) == 'ref_f3d' )then + iret=nf90_inquire_dimension(gfile_loc,1,name,len) + if(trim(name)=='xaxis_1') nx_phy=len + if( nx_phy == nxcase )then + allocate(uu2d_tmp(nxcase,nycase)) + countloc_tmp=(/nxcase,nycase,1/) + phy_smaller_domain = .false. + else + allocate(uu2d_tmp(nxcase-6,nycase-6)) + countloc_tmp=(/nxcase-6,nycase-6,1/) + phy_smaller_domain = .true. + end if + startloc_tmp=(/1,1,ilev/) + end if + + if(fv3_io_layout_y > 1) then + do nio=0,fv3_io_layout_y-1 + countloc=(/nxcase,ny_layout_len(nio),1/) + allocate(uu2d_layout(nxcase,ny_layout_len(nio))) + iret=nf90_inq_varid(gfile_loc_layout(nio),trim(adjustl(varname)),var_id) + iret=nf90_get_var(gfile_loc_layout(nio),var_id,uu2d_layout,start=startloc,count=countloc) + uu2d(:,ny_layout_b(nio):ny_layout_e(nio))=uu2d_layout + deallocate(uu2d_layout) + enddo + else + iret=nf90_inq_varid(gfile_loc,trim(adjustl(varname)),var_id) + if ( trim(adjustl(varname)) == 'ref_f3d' )then + uu2d = 0.0_r_kind + iret=nf90_get_var(gfile_loc,var_id,uu2d_tmp,start=startloc_tmp,count=countloc_tmp) + where(uu2d_tmp < 0.0_r_kind) + uu2d_tmp = 0.0_r_kind + endwhere + if(phy_smaller_domain)then + uu2d(4:nxcase-3,4:nycase-3) = uu2d_tmp + else + uu2d = uu2d_tmp + end if + deallocate(uu2d_tmp) + else + iret=nf90_get_var(gfile_loc,var_id,uu2d,start=startloc,count=countloc) + end if + endif + call fv3_h_to_ll(uu2d,hwork(:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,grid_reverse_flag) + enddo ! ilevtot + if( present(delp).or.present(tsen).or.present(w) )then ! dynvars + if(ivar == 1)then + tsen = hwork + else if(ivar == 2)then + delp = hwork + end if + if( present(w) .and. ivar == 3 )then + w = hwork + end if + end if + if( present(q).or.present(ql).or.present(qr) )then ! tracers + if(ivar == 1)then + q = hwork + else if(ivar == 2)then + oz = hwork + end if + if(present(qr))then + if(ivar == 3)then + ql = hwork + else if(ivar == 4)then + qi = hwork + else if(ivar == 5)then + qr = hwork + else if(ivar == 6)then + qs = hwork + else if(ivar == 7)then + qg = hwork + end if + end if + end if + if( present(dbz) )then ! phyvars: dbz + dbz = hwork + end if + + end do + + if(fv3_io_layout_y > 1) then + do nio=1,fv3_io_layout_y-1 + iret=nf90_close(gfile_loc_layout(nio)) + enddo + deallocate(gfile_loc_layout) + else + iret=nf90_close(gfile_loc) + endif + + deallocate (uu2d,varname_files) + end if + + return +end subroutine gsi_fv3ncdf_read_ens_parallel_over_ens + +subroutine gsi_fv3ncdf_readuv_ens_parallel_over_ens(ges_u,ges_v,fv3filenamegin,iope) +!$$$ subprogram documentation block +! . . . . +! subprogram: gsi_fv3ncdf_readuv_ens_parallel_over_ens +! program history log: +! 2022-04-01 Y. Wang and X. Wang, changed from gsi_fv3ncdf_readuv_ens +! for FV3LAM ensemble parallel IO in hybrid EnVar +! poc: xuguang.wang@ou.edu +! +! abstract: read in a field from a netcdf FV3 file in mype_u,mype_v +! then scatter the field to each PE +! program history log: +! +! input argument list: +! +! output argument list: +! ges_u - output sub domain u field +! ges_v - output sub domain v field +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ end documentation block + use kinds, only: r_kind,i_kind + use mpimod, only: mpi_comm_world,mpi_rtype,mype,mpi_info_null + use gridmod, only: nsig,nlon,nlat + use netcdf, only: nf90_open,nf90_close,nf90_get_var,nf90_noerr + use netcdf, only: nf90_nowrite,nf90_inquire,nf90_inquire_dimension + use netcdf, only: nf90_inquire_variable + use netcdf, only: nf90_inq_varid + use mod_fv3_lola, only: fv3_h_to_ll,fv3uv2earth + use general_sub2grid_mod, only: sub2grid_info,general_grid2sub + + implicit none + real(r_kind) ,intent(out ) :: ges_u(nlat,nlon,nsig) + real(r_kind) ,intent(out ) :: ges_v(nlat,nlon,nsig) + type (type_fv3regfilenameg),intent (in) :: fv3filenamegin + integer(i_kind), intent(in) :: iope + real(r_kind),dimension(2,nlat,nlon,nsig):: hwork + character(:), allocatable:: filenamein + real(r_kind),allocatable,dimension(:,:):: u2d,v2d + real(r_kind),allocatable,dimension(:,:):: uc2d,vc2d + character(len=max_varname_length) :: varname,vgsiname + real(r_kind),allocatable,dimension(:,:,:,:):: worksub + integer(i_kind) u_grd_VarId,v_grd_VarId + integer(i_kind) nlatcase,nloncase + integer(i_kind) nxcase,nycase + integer(i_kind) u_countloc(3),u_startloc(3),v_countloc(3),v_startloc(3) + integer(i_kind) inative,ilev,ilevtot + integer(i_kind) kbgn,kend + + integer(i_kind) gfile_loc,iret + integer(i_kind) nz,nzp1,mm1 + +! for fv3_io_layout_y > 1 + real(r_kind),allocatable,dimension(:,:):: u2d_layout,v2d_layout + integer(i_kind) :: nio + integer(i_kind),allocatable :: gfile_loc_layout(:) + character(len=180) :: filename_layout + + mm1=mype+1 + nloncase=nlon + nlatcase=nlat + nxcase=nx + nycase=ny + kbgn=1 + kend=nsig + if( mype == iope )then + allocate(u2d(nxcase,nycase+1)) + allocate(v2d(nxcase+1,nycase)) + allocate(uc2d(nxcase,nycase)) + allocate(vc2d(nxcase,nycase)) + filenamein=fv3filenamegin%dynvars + + if(fv3_io_layout_y > 1) then + allocate(gfile_loc_layout(0:fv3_io_layout_y-1)) + do nio=0,fv3_io_layout_y-1 + write(filename_layout,'(a,a,I4.4)') trim(filenamein),".",nio + iret=nf90_open(filename_layout,nf90_nowrite,gfile_loc_layout(nio),comm=mpi_comm_world,info=MPI_INFO_NULL) + if(iret/=nf90_noerr) then + write(6,*)'problem opening ',trim(filename_layout),gfile_loc_layout(nio),', Status = ',iret + call flush(6) + call stop2(333) + endif + enddo + else + iret=nf90_open(filenamein,nf90_nowrite,gfile_loc) + if(iret/=nf90_noerr) then + write(6,*)' problem opening ',trim(filenamein),', Status = ',iret + call flush(6) + call stop2(333) + endif + endif + do ilevtot=kbgn,kend + ilev=ilevtot + nz=nsig + nzp1=nz+1 + inative=nzp1-ilev + u_countloc=(/nxcase,nycase+1,1/) + v_countloc=(/nxcase+1,nycase,1/) + u_startloc=(/1,1,inative/) + v_startloc=(/1,1,inative/) + + if(fv3_io_layout_y > 1) then + do nio=0,fv3_io_layout_y-1 + u_countloc=(/nxcase,ny_layout_len(nio)+1,1/) + allocate(u2d_layout(nxcase,ny_layout_len(nio)+1)) + call check( nf90_inq_varid(gfile_loc_layout(nio),'u',u_grd_VarId) ) + iret=nf90_get_var(gfile_loc_layout(nio),u_grd_VarId,u2d_layout,start=u_startloc,count=u_countloc) + u2d(:,ny_layout_b(nio):ny_layout_e(nio))=u2d_layout(:,1:ny_layout_len(nio)) + if(nio==fv3_io_layout_y-1) u2d(:,ny_layout_e(nio)+1)=u2d_layout(:,ny_layout_len(nio)+1) + deallocate(u2d_layout) + + v_countloc=(/nxcase+1,ny_layout_len(nio),1/) + allocate(v2d_layout(nxcase+1,ny_layout_len(nio))) + call check( nf90_inq_varid(gfile_loc_layout(nio),'v',v_grd_VarId) ) + iret=nf90_get_var(gfile_loc_layout(nio),v_grd_VarId,v2d_layout,start=v_startloc,count=v_countloc) + v2d(:,ny_layout_b(nio):ny_layout_e(nio))=v2d_layout + deallocate(v2d_layout) + enddo + else + call check( nf90_inq_varid(gfile_loc,'u',u_grd_VarId) ) + iret=nf90_get_var(gfile_loc,u_grd_VarId,u2d,start=u_startloc,count=u_countloc) + call check( nf90_inq_varid(gfile_loc,'v',v_grd_VarId) ) + iret=nf90_get_var(gfile_loc,v_grd_VarId,v2d,start=v_startloc,count=v_countloc) + endif + + if(.not.grid_reverse_flag) then + call reverse_grid_r_uv (u2d,nxcase,nycase+1,1) + call reverse_grid_r_uv (v2d,nxcase+1,nycase,1) + endif + call fv3uv2earth(u2d(:,:),v2d(:,:),nxcase,nycase,uc2d,vc2d) + + ! NOTE on transfor to earth u/v: + ! The u and v before transferring need to be in E-W/N-S grid, which is + ! defined as reversed grid here because it is revered from map view. + ! + ! Have set the following flag for grid orientation + ! grid_reverse_flag=true: E-W/N-S grid + ! grid_reverse_flag=false: W-E/S-N grid + ! + ! So for preparing the wind transferring, need to reverse the grid + ! from + ! W-E/S-N grid to E-W/N-S grid when grid_reverse_flag=false: + ! + ! if(.not.grid_reverse_flag) call reverse_grid_r_uv + ! + ! and the last input parameter for fv3_h_to_ll is alway true: + ! + ! + call fv3_h_to_ll(uc2d,hwork(1,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,.true.) + call fv3_h_to_ll(vc2d,hwork(2,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,.true.) + enddo ! ilevtot + if(fv3_io_layout_y > 1) then + do nio=0,fv3_io_layout_y-1 + iret=nf90_close(gfile_loc_layout(nio)) + enddo + deallocate(gfile_loc_layout) + else + iret=nf90_close(gfile_loc) + endif + deallocate(u2d,v2d,uc2d,vc2d) + ges_u = hwork(1,:,:,:) + ges_v = hwork(2,:,:,:) + end if ! mype + +end subroutine gsi_fv3ncdf_readuv_ens_parallel_over_ens + + subroutine wrfv3_netcdf(fv3filenamegin) !$$$ subprogram documentation block ! . . . . @@ -2653,7 +3140,8 @@ subroutine wrfv3_netcdf(fv3filenamegin) ! program history log: ! 2019-04-18 CAPS(C. Tong) - import direct reflectivity DA capabilities ! 2019-11-22 CAPS(C. Tong) - modify "add_saved" to properly output analyses -! 2021-01-05 x.zhang/lei - add code for updating delz analysis in regional da +! 2021-01-05 x.zhang/lei - add code for updating delz analysis in regional da +! 2022-04-01 Y. Wang and X. Wang - add code for updating reflectivity ! ! input argument list: ! @@ -2680,6 +3168,7 @@ subroutine wrfv3_netcdf(fv3filenamegin) use directDA_radaruse_mod, only: l_cvpnr, cvpnr_pval use gridmod, only: eta1_ll,eta2_ll use constants, only: one + use obsmod, only: if_model_dbz implicit none @@ -2705,6 +3194,7 @@ subroutine wrfv3_netcdf(fv3filenamegin) real(r_kind),pointer,dimension(:,:,:):: ges_qg =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_qnr =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_w =>NULL() + real(r_kind),pointer,dimension(:,:,:):: ges_dbz =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_delzinc =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_delp =>NULL() real(r_kind),dimension(:,: ),allocatable:: ges_ps_write @@ -2794,14 +3284,17 @@ subroutine wrfv3_netcdf(fv3filenamegin) call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'u' , ges_u ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'v' , ges_v ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'q' ,ges_q ,istatus);ier=ier+istatus - if (l_use_dbz_directDA) then + if (l_use_dbz_directDA .or. if_model_dbz) then call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'ql' ,ges_ql ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qi' ,ges_qi ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qr' ,ges_qr ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qs' ,ges_qs ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qg' ,ges_qg ,istatus);ier=ier+istatus + if (l_use_dbz_directDA) & call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qnr',ges_qnr,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'w' , ges_w ,istatus);ier=ier+istatus + if( if_model_dbz )& + call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'dbz' , ges_dbz ,istatus);ier=ier+istatus end if if(i_use_2mq4b > 0 .and. i_use_2mt4b > 0 ) then call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'q2m',ges_q2m,istatus); ier=ier+istatus @@ -2931,6 +3424,7 @@ subroutine wrfv3_netcdf(fv3filenamegin) call gsi_copy_bundle(GSI_MetGuess_Bundle(it),gsibundle_fv3lam_dynvar_nouv) call gsi_copy_bundle(GSI_MetGuess_Bundle(it),gsibundle_fv3lam_tracer_nouv) + if( if_model_dbz ) call gsi_copy_bundle(GSI_MetGuess_Bundle(it),gsibundle_fv3lam_phyvar_nouv) if (laeroana_fv3cmaq) then call gsi_copy_bundle(GSI_ChemGuess_Bundle(it),gsibundle_fv3lam_tracerchem_nouv) end if @@ -2979,6 +3473,10 @@ subroutine wrfv3_netcdf(fv3filenamegin) add_saved,fv3filenamegin%dynvars,fv3filenamegin) call gsi_fv3ncdf_write(grd_fv3lam_tracer_ionouv,gsibundle_fv3lam_tracer_nouv, & add_saved,fv3filenamegin%tracers,fv3filenamegin) + if( if_model_dbz ) then + call gsi_fv3ncdf_write(grd_fv3lam_phyvar_ionouv,gsibundle_fv3lam_phyvar_nouv,& + add_saved,fv3filenamegin%phyvars,fv3filenamegin) + end if call gsi_fv3ncdf_writeuv(grd_fv3lam_uv,ges_u,ges_v,add_saved,fv3filenamegin) if (laeroana_fv3cmaq) then call gsi_fv3ncdf_write(grd_fv3lam_tracerchem_ionouv,gsibundle_fv3lam_tracerchem_nouv, & @@ -3598,6 +4096,7 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file use netcdf, only: nf90_open,nf90_close use netcdf, only: nf90_write,nf90_inq_varid use netcdf, only: nf90_put_var,nf90_get_var + use netcdf, only: nf90_inquire_dimension use gsi_bundlemod, only: gsi_bundle use general_sub2grid_mod, only: sub2grid_info,general_sub2grid implicit none @@ -3609,16 +4108,19 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file type (type_fv3regfilenameg),intent (in) :: fv3filenamegin real(r_kind),dimension(1,grd_ionouv%nlat,grd_ionouv%nlon,grd_ionouv%kbegin_loc:grd_ionouv%kend_alloc):: hwork character(len=max_varname_length) :: filenamein2 - character(len=max_varname_length) :: varname,vgsiname + character(len=max_varname_length) :: varname,vgsiname,name integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(3),startloc(3) + integer(i_kind) countloc_tmp(3),startloc_tmp(3) integer(i_kind) kbgn,kend integer(i_kind) inative,ilev,ilevtot integer(i_kind) :: VarId,gfile_loc - integer(i_kind) mm1,nzp1 + integer(i_kind) mm1,nzp1,len,nx_phy,iret + logical :: phy_smaller_domain real(r_kind),allocatable,dimension(:,:):: work_a real(r_kind),allocatable,dimension(:,:):: work_b real(r_kind),allocatable,dimension(:,:):: workb2,worka2 + real(r_kind),allocatable,dimension(:,:):: work_b_tmp ! for io_layout > 1 real(r_kind),allocatable,dimension(:,:):: work_b_layout @@ -3673,7 +4175,20 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file work_a=hwork(1,:,:,ilevtot) - + if( trim(varname) == 'ref_f3d' )then + iret=nf90_inquire_dimension(gfile_loc,1,name,len) + if(trim(name)=='xaxis_1') nx_phy=len + if( nx_phy == nxcase )then + allocate(work_b_tmp(nxcase,nycase)) + countloc_tmp=(/nxcase,nycase,1/) + phy_smaller_domain = .false. + else + allocate(work_b_tmp(nxcase-6,nycase-6)) + countloc_tmp=(/nxcase-6,nycase-6,1/) + phy_smaller_domain = .true. + end if + startloc_tmp=(/1,1,ilev/) + end if call check( nf90_inq_varid(gfile_loc,trim(varname),VarId) ) @@ -3703,16 +4218,29 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file deallocate(work_b_layout) enddo else - call check( nf90_get_var(gfile_loc,VarId,work_b,start = startloc, count = countloc) ) + if( trim(varname) == 'ref_f3d' )then + work_b = 0.0_r_kind + call check( nf90_get_var(gfile_loc,VarId,work_b_tmp,start = startloc_tmp, count = countloc_tmp) ) + where(work_b_tmp < 0.0_r_kind) + work_b_tmp = 0.0_r_kind + end where + if(phy_smaller_domain)then + work_b(4:nxcase-3,4:nycase-3) = work_b_tmp + else + work_b(1:nxcase,1:nycase) = work_b_tmp + end if + else + call check( nf90_get_var(gfile_loc,VarId,work_b,start = startloc, count = countloc) ) + end if endif call fv3_h_to_ll(work_b(:,:),worka2,nlon_regional,nlat_regional,nloncase,nlatcase,grid_reverse_flag) !!!!!!!! analysis_inc: work_a !!!!!!!!!!!!!!!! work_a(:,:)=work_a(:,:)-worka2(:,:) call fv3_ll_to_h(work_a(:,:),workb2,nloncase,nlatcase,nlon_regional,nlat_regional,grid_reverse_flag) work_b(:,:)=work_b(:,:)+workb2(:,:) - else + else call fv3_ll_to_h(work_a(:,:),work_b(:,:),nloncase,nlatcase,nlon_regional,nlat_regional,grid_reverse_flag) - endif + endif endif if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 @@ -3723,7 +4251,20 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file deallocate(work_b_layout) enddo else - call check( nf90_put_var(gfile_loc,VarId,work_b, start = startloc, count = countloc) ) + if( trim(varname) == 'ref_f3d' )then + if(phy_smaller_domain)then + work_b_tmp = work_b(4:nxcase-3,4:nycase-3) + else + work_b_tmp = work_b(1:nxcase,1:nycase) + end if + where(work_b_tmp < 0.0_r_kind) + work_b_tmp = 0.0_r_kind + end where + call check( nf90_put_var(gfile_loc,VarId,work_b_tmp, start = startloc_tmp, count = countloc_tmp) ) + deallocate(work_b_tmp) + else + call check( nf90_put_var(gfile_loc,VarId,work_b, start = startloc, count = countloc) ) + end if endif enddo !ilevtotl loop @@ -4443,6 +4984,8 @@ subroutine getfv3lamfilevname(vgsinamein,fv3filenamegref,filenameout,vname) filenameout=fv3filenamegref%dynvars else if(ifindstrloc(vartracers,vgsinamein)> 0 ) then filenameout=fv3filenamegref%tracers + else if(ifindstrloc(varphyvars,vgsinamein)> 0) then + filenameout=fv3filenamegref%phyvars else write(6,*)'the filename corresponding to var ',trim(vgsinamein),' is not found, stop ' call stop2(333) diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 55d9298614..8ccef7c38e 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -151,7 +151,9 @@ module gsimod readin_beta,use_localization_grid,use_gfs_ens,q_hyb_ens,i_en_perts_io, & l_ens_in_diff_time,ensemble_path,ens_fast_read,sst_staticB,limqens, & ntotensgrp,nsclgrp,naensgrp,ngvarloc,ntlevs_ens,naensloc, & - i_ensloccov4tim,i_ensloccov4var,i_ensloccov4scl,l_timloc_opt + i_ensloccov4tim,i_ensloccov4var,i_ensloccov4scl,l_timloc_opt,& + vdl_scale,vloc_varlist,& + global_spectral_filter_sd,assign_vdl_nml,parallelization_over_ensmembers use hybrid_ensemble_parameters,only : l_both_fv3sar_gfs_ens,n_ens_gfs,n_ens_fv3sar use rapidrefresh_cldsurf_mod, only: init_rapidrefresh_cldsurf, & dfi_radar_latent_heat_time_period,metar_impact_radius,& @@ -1375,6 +1377,36 @@ module gsimod ! i_ensloccov4scl - flag of cross-scale localization ! =0: cross-scale covariance is retained ! =1: cross-scale covariance is zero +! +! global_spectral_filter_sd - if true, use spectral filter function for +! scale decomposition in the global application (Huang et al. 2021) +! assign_vdl_nml - if true, vdl_scale, and vloc_varlist will be used for +! assigning variable-dependent localization upon SDL in gsiparm.anl. +! This method described in (Wang and Wang 2022, JAMES) is +! equivalent to, but different from the method associated +! with the parameter i_ensloccov4var. +! vloc_varlist - list of control variables using the same localization length, +! effective only with assign_vdl_nml=.true. For example, +! vloc_varlist(1,:) = 'sf','vp','ps','t', +! vloc_varlist(2,:) = 'q', +! vloc_varlist(3,:) = 'qr','qs','qg','dbz','w','ql','qi', +! vloc_varlist(4,:) = 'sf','vp','ps','t','q', +! vloc_varlist(5,:) = 'qr','qs','qg','dbz','w','ql','qi', +! This example indicates that 3 variable-groups will be adopted for VDL. +! 'sf','vp','ps','t' will share the same localization length of v1L1; +! 'q' will have the localization lenth of v2L1 +! 'qr','qs','qg','dbz','w','ql','qi', use the same localization length of v3L1 +! +! For L2, a different configuration of VDL can be applied: +! ~~~~~~~~~ +! 'sf','vp','ps','t','q' will share the same localization length of v2L2; +! 'qr','qs','qg','dbz','w','ql','qi', use the same localization length of v2L2 +! vdl_scale - number of variables in each variable-group, effective only with assign_vdl_nml=.true. +! if 3 variable-groups with 2 separated scale is set, +! vdl_scale = 3, 3, 3, 2, 2 +! ^ ^ ^ ^ ^ +! s_ens_h = v1L1 v2L1 v3L1 v1L2 v2L2 +! Then localization lengths will be assigned as above. ! namelist/hybrid_ensemble/l_hyb_ens,uv_hyb_ens,q_hyb_ens,aniso_a_en,generate_ens,n_ens,l_both_fv3sar_gfs_ens,n_ens_gfs,n_ens_fv3sar,nlon_ens,nlat_ens,jcap_ens,& pseudo_hybens,merge_two_grid_ensperts,regional_ensemble_option,fv3sar_bg_opt,fv3sar_ensemble_opt,full_ensemble,pwgtflg,& @@ -1382,7 +1414,9 @@ module gsimod grid_ratio_ens, & oz_univ_static,write_ens_sprd,use_localization_grid,use_gfs_ens, & i_en_perts_io,l_ens_in_diff_time,ensemble_path,ens_fast_read,sst_staticB,limqens, & - nsclgrp,l_timloc_opt,ngvarloc,naensloc,i_ensloccov4tim,i_ensloccov4var,i_ensloccov4scl + nsclgrp,l_timloc_opt,ngvarloc,naensloc,i_ensloccov4tim,i_ensloccov4var,i_ensloccov4scl,& + vdl_scale,vloc_varlist,& + global_spectral_filter_sd,assign_vdl_nml,parallelization_over_ensmembers ! rapidrefresh_cldsurf (options for cloud analysis and surface ! enhancement for RR appilcation ): @@ -1839,7 +1873,9 @@ subroutine gsimain_initialize else naensgrp=ntotensgrp endif - if(naensloc 1)then + allocate(spc_multwgt(0:jcap_ens,nsclgrp)) + allocate(spcwgt_params(4,nsclgrp)) + spc_multwgt=1.0 + + ! The below parameters are used in Huang et al. (2021, MWR) + spcwgt_params(1,1)=4000.0_r_kind + spcwgt_params(2,1)=100000000.0_r_kind + spcwgt_params(3,1)=1.0_r_kind + spcwgt_params(4,1)=3000.0_r_kind + + if( nsclgrp >=3 )then + spcwgt_params(1,3)=0.0_r_kind + spcwgt_params(2,3)=500.0_r_kind + spcwgt_params(3,3)=1.0_r_kind + spcwgt_params(4,3)=500.0_r_kind + end if + + call init_mult_spc_wgts(jcap_ens) + + end if + return end subroutine hybens_grid_setup @@ -4148,6 +4171,8 @@ subroutine hybens_localization_setup ! 2012-10-16 wu - only call setup_ens_wgt if necessary ! 2014-05-22 wu modification to allow vertically varying localization scales in regional ! 2022-09-15 yokota - add scale/variable/time-dependent localization +! 2022-12-09 Y. Wang and X. Wang - add a variable-dependent localization option (assign_vdl_nml=.true.), +! poc: xuguang.wang@ou.edu ! ! input argument list: ! @@ -4168,9 +4193,11 @@ subroutine hybens_localization_setup use hybrid_ensemble_parameters, only: readin_beta,beta_s,beta_e,beta_s0,beta_e0,sqrt_beta_s,sqrt_beta_e use hybrid_ensemble_parameters, only: readin_localization,create_hybens_localization_parameters, & vvlocal,s_ens_h,s_ens_hv,s_ens_v,s_ens_vv - use hybrid_ensemble_parameters, only: ntotensgrp,naensgrp,naensloc,ntlevs_ens,nsclgrp - use hybrid_ensemble_parameters, only: en_perts + use hybrid_ensemble_parameters, only: ntotensgrp,naensgrp,naensloc,ntlevs_ens,nsclgrp,assign_vdl_nml + use hybrid_ensemble_parameters, only: en_perts,vdl_scale,vloc_varlist,global_spectral_filter_sd + use hybrid_ensemble_parameters, only: ngvarloc use gsi_io, only: verbose + use string_utility, only: StrLowCase implicit none @@ -4183,8 +4210,11 @@ subroutine hybens_localization_setup real(r_kind),allocatable:: s_ens_h_gu_x(:,:),s_ens_h_gu_y(:,:) logical :: l_read_success type(gsi_bundle) :: a_en(n_ens) + type(gsi_bundle) :: en_pertstmp(n_ens,ntlevs_ens) + type(gsi_bundle) :: en_pertstmp1(n_ens,ntlevs_ens) type(gsi_grid) :: grid_ens real(r_kind), pointer :: values(:) => NULL() + integer(i_kind) :: iscl, iv, smooth_scales_num character(len=*),parameter::myname_=myname//'*hybens_localization_setup' l_read_success=.false. @@ -4312,7 +4342,7 @@ subroutine hybens_localization_setup call init_sf_xy(jcap_ens) endif - if(ntotensgrp>1) then + if(ntotensgrp>1 .and. (.not. global_spectral_filter_sd)) then call gsi_bundlegetpointer(en_perts(1,1,1),cvars3d,ipc3d,istatus) if(istatus/=0) then write(6,*) myname_,': cannot find 3d pointers' @@ -4326,64 +4356,159 @@ subroutine hybens_localization_setup if(nsclgrp>1) then call gsi_gridcreate(grid_ens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig) allocate(values(grd_ens%latlon11*grd_ens%nsig*n_ens)) - do ig=1,nsclgrp-1 - ii=0 - do n=1,n_ens - a_en(n)%values => values(ii+1:ii+grd_ens%latlon11*grd_ens%nsig) - call gsi_bundleset(a_en(n),grid_ens,'Ensemble Bundle',istatus,names3d=(/'a_en'/),bundle_kind=r_kind) - if (istatus/=0) then - write(6,*) myname_,': error alloc(ensemble bundle)' - call stop2(999) - endif - ii=ii+grd_ens%latlon11*grd_ens%nsig - enddo - do m=1,ntlevs_ens + if( .not. assign_vdl_nml )then + do ig=1,nsclgrp-1 + ii=0 do n=1,n_ens - en_perts(n,ig+1,m)%valuesr4=en_perts(n,ig,m)%valuesr4 + a_en(n)%values => values(ii+1:ii+grd_ens%latlon11*grd_ens%nsig) + call gsi_bundleset(a_en(n),grid_ens,'Ensemble Bundle',istatus,names3d=(/'a_en'/),bundle_kind=r_kind) + if (istatus/=0) then + write(6,*) myname_,': error alloc(ensemble bundle)' + call stop2(999) + endif + ii=ii+grd_ens%latlon11*grd_ens%nsig enddo - do ic3=1,nc3d - ipic=ipc3d(ic3) + do m=1,ntlevs_ens do n=1,n_ens - do k=1,grd_ens%nsig - a_en(n)%r3(1)%q(:,:,k)=en_perts(n,ig,m)%r3(ipic)%qr4(:,:,k) + en_perts(n,ig+1,m)%valuesr4=en_perts(n,ig,m)%valuesr4 + enddo + do ic3=1,nc3d + ipic=ipc3d(ic3) + do n=1,n_ens + do k=1,grd_ens%nsig + a_en(n)%r3(1)%q(:,:,k)=en_perts(n,ig,m)%r3(ipic)%qr4(:,:,k) + enddo + enddo + call bkgcov_a_en_new_factorization(naensgrp+ig,a_en) + do n=1,n_ens + do k=1,grd_ens%nsig + en_perts(n,ig,m)%r3(ipic)%qr4(:,:,k)=a_en(n)%r3(1)%q(:,:,k) + enddo enddo enddo - call bkgcov_a_en_new_factorization(naensgrp+ig,a_en) - do n=1,n_ens - do k=1,grd_ens%nsig - en_perts(n,ig,m)%r3(ipic)%qr4(:,:,k)=a_en(n)%r3(1)%q(:,:,k) + do ic2=1,nc2d + ipic=ipc2d(ic2) + do n=1,n_ens + do k=1,grd_ens%nsig + a_en(n)%r3(1)%q(:,:,k)=en_perts(n,ig,m)%r2(ipic)%qr4(:,:) + enddo + enddo + call bkgcov_a_en_new_factorization(naensgrp+ig,a_en) + do n=1,n_ens + en_perts(n,ig,m)%r2(ipic)%qr4(:,:)=a_en(n)%r3(1)%q(:,:,1) enddo enddo - enddo - do ic2=1,nc2d - ipic=ipc2d(ic2) do n=1,n_ens - do k=1,grd_ens%nsig - a_en(n)%r3(1)%q(:,:,k)=en_perts(n,ig,m)%r2(ipic)%qr4(:,:) + en_perts(n,ig+1,m)%valuesr4=en_perts(n,ig+1,m)%valuesr4-en_perts(n,ig,m)%valuesr4 + enddo + enddo + do n=1,n_ens + call gsi_bundleunset(a_en(n),istatus) + enddo + enddo + else ! assign_vdl_nml + smooth_scales_num = naensloc - naensgrp + ngvarloc = 1 ! forced to 1 in this option + do n = 1, n_ens + do m = 1, ntlevs_ens + call gsi_bundlecreate(en_pertstmp(n,m),grid_ens,'ensemble2',istatus,names2d=cvars2d,names3d=cvars3d,bundle_kind=r_single) + call gsi_bundlecreate(en_pertstmp1(n,m),grid_ens,'ensemble1',istatus,names2d=cvars2d,names3d=cvars3d,bundle_kind=r_single) + end do + end do + ig = 1 + do iscl=1,smooth_scales_num + 1 + ii=0 + do n=1,n_ens + a_en(n)%values => values(ii+1:ii+grd_ens%latlon11*grd_ens%nsig) + call gsi_bundleset(a_en(n),grid_ens,'Ensemble Bundle',istatus,names3d=(/'a_en'/),bundle_kind=r_kind) + if (istatus/=0) then + write(6,*) myname_,': error alloc(ensemble bundle)' + call stop2(999) + endif + ii=ii+grd_ens%latlon11*grd_ens%nsig + enddo + + do m=1,ntlevs_ens + if( ig == 1 )then + do n=1,n_ens + en_pertstmp(n,m)%valuesr4=en_perts(n,ig,m)%valuesr4 + enddo + end if + do ic3=1,nc3d + ipic=ipc3d(ic3) + do n=1,n_ens + do k=1,grd_ens%nsig + a_en(n)%r3(1)%q(:,:,k)=en_pertstmp(n,m)%r3(ipic)%qr4(:,:,k) + enddo + enddo + if(iscl <= smooth_scales_num) call bkgcov_a_en_new_factorization(naensgrp+iscl,a_en) + do n=1,n_ens + do k=1,grd_ens%nsig + en_pertstmp1(n,m)%r3(ipic)%qr4(:,:,k)=a_en(n)%r3(1)%q(:,:,k) + if( vdl_scale(ig) == 0 )then + en_perts(n,ig,m)%r3(ipic)%qr4(:,:,k)=a_en(n)%r3(1)%q(:,:,k) + else ! VDL is activated + do iv = 1, vdl_scale(ig) + en_perts(n,ig+iv-1,m)%r3(ipic)%qr4(:,:,k)=0.0_r_single + if( any( trim(StrLowCase(cvars3d(ic3))) == vloc_varlist(ig+iv-1,:) ) ) then + en_perts(n,ig+iv-1,m)%r3(ipic)%qr4(:,:,k)=a_en(n)%r3(1)%q(:,:,k) + end if + end do + end if + enddo + enddo + enddo + do ic2=1,nc2d + ipic=ipc2d(ic2) + do n=1,n_ens + do k=1,grd_ens%nsig + a_en(n)%r3(1)%q(:,:,k)=en_pertstmp(n,m)%r2(ipic)%qr4(:,:) + enddo + enddo + if(iscl <= smooth_scales_num) call bkgcov_a_en_new_factorization(naensgrp+iscl,a_en) + do n=1,n_ens + en_pertstmp1(n,m)%r2(ipic)%qr4(:,:)=a_en(n)%r3(1)%q(:,:,1) + if( vdl_scale(ig) == 0 )then + en_perts(n,ig,m)%r2(ipic)%qr4(:,:)=a_en(n)%r3(1)%q(:,:,1) + else ! VDL is activated + do iv = 1, vdl_scale(ig) + en_perts(n,ig+iv-1,m)%r2(ipic)%qr4(:,:)=0.0_r_single + if( any( trim(StrLowCase(cvars2d(ic2))) == vloc_varlist(ig+iv-1,:) ) ) then + en_perts(n,ig+iv-1,m)%r2(ipic)%qr4(:,:)=a_en(n)%r3(1)%q(:,:,1) + end if + end do + end if enddo enddo - call bkgcov_a_en_new_factorization(naensgrp+ig,a_en) do n=1,n_ens - en_perts(n,ig,m)%r2(ipic)%qr4(:,:)=a_en(n)%r3(1)%q(:,:,1) + en_pertstmp(n,m)%valuesr4=en_pertstmp(n,m)%valuesr4-en_pertstmp1(n,m)%valuesr4 enddo enddo do n=1,n_ens - en_perts(n,ig+1,m)%valuesr4=en_perts(n,ig+1,m)%valuesr4-en_perts(n,ig,m)%valuesr4 + call gsi_bundleunset(a_en(n),istatus) enddo + if( vdl_scale(ig) == 0 )then + ig = ig + 1 + else + ig = ig + vdl_scale(ig) + end if enddo do n=1,n_ens - call gsi_bundleunset(a_en(n),istatus) - enddo - enddo - deallocate(values) - endif - do ig=nsclgrp+1,ntotensgrp - do m=1,ntlevs_ens - do n=1,n_ens - en_perts(n,ig,m)%valuesr4=en_perts(n,ig-nsclgrp,m)%valuesr4 - enddo - enddo - enddo + do m=1,ntlevs_ens + call gsi_bundledestroy(en_pertstmp(n,m),istatus) + call gsi_bundledestroy(en_pertstmp1(n,m),istatus) + end do + end do + end if + deallocate(values) + endif + do ig=nsclgrp+1,ntotensgrp + do m=1,ntlevs_ens + do n=1,n_ens + en_perts(n,ig,m)%valuesr4=en_perts(n,ig-nsclgrp,m)%valuesr4 + enddo + enddo + enddo endif !!!!!!!! setup beta_s, beta_e!!!!!!!!!!!! diff --git a/src/gsi/hybrid_ensemble_parameters.f90 b/src/gsi/hybrid_ensemble_parameters.f90 index 7b1c963764..17416f68fb 100644 --- a/src/gsi/hybrid_ensemble_parameters.f90 +++ b/src/gsi/hybrid_ensemble_parameters.f90 @@ -128,6 +128,7 @@ module hybrid_ensemble_parameters ! function of z, default = .false. ! ensemble_path: path to ensemble members; default './' ! ens_fast_read: read ensemble in parallel; default '.false.' +! parallelization_over_ensmembers: parallelly read ensemble members for FV3-LAM; default '.false' ! sst_staticB: if .true. (default) uses only static part of B error covariance for SST ! nsclgrp: number of scale-dependent localization lengths ! l_timloc_opt: if true, then turn on time-dependent localization @@ -327,10 +328,17 @@ module hybrid_ensemble_parameters public :: i_ensloccov4tim,i_ensloccov4var,i_ensloccov4scl public :: idaen3d,idaen2d public :: ens_fast_read + public :: parallelization_over_ensmembers public :: l_both_fv3sar_gfs_ens public :: sst_staticB public :: limqens + public :: spc_multwgt + public :: spcwgt_params + public :: vdl_scale,vloc_varlist + public :: global_spectral_filter_sd + public :: assign_vdl_nml + logical l_hyb_ens,uv_hyb_ens,q_hyb_ens,oz_univ_static,sst_staticB logical l_timloc_opt logical aniso_a_en @@ -348,6 +356,7 @@ module hybrid_ensemble_parameters logical vvlocal logical l_ens_in_diff_time logical ens_fast_read + logical parallelization_over_ensmembers logical l_both_fv3sar_gfs_ens integer(i_kind) i_en_perts_io integer(i_kind) n_ens,nlon_ens,nlat_ens,jcap_ens,jcap_ens_test @@ -389,6 +398,13 @@ module hybrid_ensemble_parameters integer(i_kind) :: i_ensloccov4scl=0 integer(i_kind),allocatable,dimension(:) :: idaen3d,idaen2d + real(r_kind),allocatable,dimension(:,:) :: spc_multwgt + real(r_kind),allocatable,dimension(:,:) :: spcwgt_params + character(len=3) vloc_varlist(max_naensloc,max_nvars) + integer(i_kind) vdl_scale(max_naensloc) + logical :: global_spectral_filter_sd + logical :: assign_vdl_nml + ! following is for storage of ensemble perturbations: ! def en_perts - array of ensemble perturbations @@ -476,10 +492,15 @@ subroutine init_hybrid_ensemble_parameters i_en_perts_io=0 ! default for en_pert IO. 0 is no IO ensemble_path = './' ! default for path to ensemble members ens_fast_read=.false. + parallelization_over_ensmembers=.false. limqens=1.0_r_single ! default for limiting ensemble RH (+/-) l_both_fv3sar_gfs_ens=.false. n_ens_gfs=0 n_ens_fv3sar=0 + vdl_scale = 0 + vloc_varlist = 'aaa' + global_spectral_filter_sd=.false. + assign_vdl_nml=.false. end subroutine init_hybrid_ensemble_parameters diff --git a/src/gsi/obsmod.F90 b/src/gsi/obsmod.F90 index 1fb03d0940..3dd936d94e 100644 --- a/src/gsi/obsmod.F90 +++ b/src/gsi/obsmod.F90 @@ -747,7 +747,7 @@ subroutine init_obsmod_dflts if_vterminal=.false. l2rwthin =.false. if_vrobs_raw=.false. - if_model_dbz=.true. + if_model_dbz=.false. inflate_obserr=.false. whichradar="KKKK" diff --git a/src/gsi/read_dbz_nc.f90 b/src/gsi/read_dbz_nc.f90 index 89eebde8b6..ee1d3cb2e4 100644 --- a/src/gsi/read_dbz_nc.f90 +++ b/src/gsi/read_dbz_nc.f90 @@ -71,7 +71,7 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no one_tenth,r1000,r60,r60inv,r100,r400,grav_equator, & eccentricity,somigliana,grav_ratio,grav,semi_major_axis,flattening use gridmod, only: tll2xy,nsig,nlat,nlon - use obsmod, only: iadate,doradaroneob, & + use obsmod, only: iadate,doradaroneob,oneoblat,oneoblon,oneobheight,oneobradid, & mintiltdbz,maxtiltdbz,minobrangedbz,maxobrangedbz,& static_gsi_nopcp_dbz,rmesh_dbz,zmesh_dbz use hybrid_ensemble_parameters,only : l_hyb_ens @@ -380,6 +380,12 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no thislon = lon(i,j) thislat = lat(i,j) + + if(doradaroneob) then + thislat=oneoblat + thislon=oneoblon + thishgt=oneobheight + endif !-Check format of longitude and correct if necessary diff --git a/src/gsi/setupdbz.f90 b/src/gsi/setupdbz.f90 index 9bbf5ed34b..1e158de9ea 100644 --- a/src/gsi/setupdbz.f90 +++ b/src/gsi/setupdbz.f90 @@ -1451,31 +1451,37 @@ subroutine setupdbz(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,radardbz_d ! Write information to diagnostic file if(radardbz_diagsave .and. ii>0 )then - write(string,600) jiter -600 format('radardbz_',i2.2) - diag_file=trim(dirname) // trim(string) - if(init_pass) then - open(newunit=lu_diag,file=trim(diag_file),form='unformatted',status='unknown',position='rewind') + if( .not. l_use_dbz_directDA )then + write(7)'dbz',nchar,nreal,ii,mype,ioff0 + write(7)cdiagbuf(1:ii),rdiagbuf(:,1:ii) + deallocate(cdiagbuf,rdiagbuf) else - inquire(file=trim(diag_file),exist=diagexist) - if (diagexist) then - open(lu_diag,file=trim(diag_file),form='unformatted',status='old',position='append') + write(string,600) jiter +600 format('radardbz_',i2.2) + diag_file=trim(dirname) // trim(string) + if(init_pass) then + open(newunit=lu_diag,file=trim(diag_file),form='unformatted',status='unknown',position='rewind') else - open(lu_diag,file=trim(diag_file),form='unformatted',status='unknown',position='rewind') + inquire(file=trim(diag_file),exist=diagexist) + if (diagexist) then + open(lu_diag,file=trim(diag_file),form='unformatted',status='old',position='append') + else + open(lu_diag,file=trim(diag_file),form='unformatted',status='unknown',position='rewind') + endif + endif + if(init_pass .and. mype == 0) then + if ( .not. l_use_dbz_directDA ) then ! EnKF uses these diagnostics and EnKF uses single OBS file for now. + write(lu_diag) ianldate ! So do not write analysis date for binary in case of using direct reflectivity DA. + end if + write(6,*)'SETUPDBZ: write time record to file ',& + trim(diag_file), ' ',ianldate endif - endif - if(init_pass .and. mype == 0) then - if ( .not. l_use_dbz_directDA ) then ! EnKF uses these diagnostics and EnKF uses single OBS file for now. - write(lu_diag) ianldate ! So do not write analysis date for binary in case of using direct reflectivity DA. - end if - write(6,*)'SETUPDBZ: write time record to file ',& - trim(diag_file), ' ',ianldate - endif - write(lu_diag)'dbz',nchar,nreal,ii,mype,ioff0 - write(lu_diag)cdiagbuf(1:ii),rdiagbuf(:,1:ii) - deallocate(cdiagbuf,rdiagbuf) - close(lu_diag) + write(lu_diag)'dbz',nchar,nreal,ii,mype,ioff0 + write(lu_diag)cdiagbuf(1:ii),rdiagbuf(:,1:ii) + deallocate(cdiagbuf,rdiagbuf) + close(lu_diag) + end if end if write(6,*)'mype, irefsmlobs,irejrefsmlobs are ',mype,' ',irefsmlobs, ' ',irejrefsmlobs ! close(52) !simulated obs diff --git a/src/gsi/update_guess.f90 b/src/gsi/update_guess.f90 index a90e7a19d6..e5a0f64245 100644 --- a/src/gsi/update_guess.f90 +++ b/src/gsi/update_guess.f90 @@ -287,7 +287,11 @@ subroutine update_guess(sval,sbias) ! since we don't know which comes first in met-guess, we ! must postpone updating tv after all other met-guess fields endif - icloud=getindex(cloud,guess(ic)) + if( allocated(cloud) )then + icloud=getindex(cloud,guess(ic)) + else + icloud=-999 + end if if ( .not. l_use_dbz_directDA ) then ! original code if(icloud>0) then ptr3dges = max(ptr3dges+ptr3dinc,zero)