Add options for Mixed Ensemble and Scale/Variable-Dependent Localizat…

…ion (NOAA-GSL#451)

* Add options for Scale/Variable/Time-Dependent Localization

* Add options for Variable-Dependent Localization of OU/MAP

* rename 'all' to 'conv_dbz' and separate configuration files to set_rrfs_config_SDL_VDL_MixEn.sh

* remove change for preventing from crashing restart files for now

---------

Co-authored-by: Sho Yokota <syokota@Orion-login-1.HPC.MsState.Edu>

fix/gsi/anavinfo.rrfs_conv_dbz

@@ -0,0 +1,81 @@
+met_guess::
+!var     level    crtm_use    desc              orig_name
+  ps        1      -1         surface_pressure     ps
+  z         1      -1         geopotential_height  phis
+  u        65       2         zonal_wind           u
+  v        65       2         meridional_wind      v
+  w        65       2         vertical velocity    w
+  tv       65       2         virtual_temperature  tv
+  q        65       2         specific_humidity    sphu
+  oz       65       2         ozone                ozone
+  delp     65      -1         fv3 del pressure     delp
+  t2m       1      -1         2-m_T                t2m
+  q2m       1      -1         2-m_Q                q2m
+  ql       65      10         cloud_condensate     QCLOUD
+  qr       65      10         rain                 QRAIN
+  qs       65      10         snow                 QSNOW
+  qi       65      10         ice                  QICE
+  qg       65      10         graupel              QGRAUP
+  dbz      65      10         reflectivity         REFL_10CM
+::
+
+state_derivatives::
+!var  level  src
+ ps   1      met_guess
+ u    65     met_guess
+ v    65     met_guess
+ tv   65     met_guess
+ q    65     met_guess
+ oz   65     met_guess
+ cw   65     met_guess
+ prse 66     met_guess
+::
+
+state_tendencies::
+!var  levels  source
+ u    65      met_guess
+ v    65      met_guess
+ tv   65      met_guess
+ q    65      met_guess
+ oz   65      met_guess
+ prse 66      met_guess
+::
+
+state_vector::  
+!var     level  itracer source     funcof
+ u        65      0     met_guess    u
+ v        65      0     met_guess    v
+ w        65      0     met_guess    w
+ tv       65      0     met_guess    tv
+ tsen     65      0     met_guess    tv,q
+ q        65      1     met_guess    q
+ oz       65      1     met_guess    oz
+ prse     66      0     met_guess    prse
+ ps        1      0     met_guess    prse
+ sst       1      0     met_guess    sst
+ qr       65      1     met_guess    qr
+ ql       65      1     met_guess    qc
+ qs       65      1     met_guess    qs
+ qi       65      1     met_guess    qi
+ qg       65      1     met_guess    qg
+ dbz      65      1     met_guess    dbz
+::
+
+control_vector::
+!var     level  itracer as/tsfc_sdv  an_amp0   source  funcof
+ sf       65      0       0.90        -1.0     state    u,v
+ vp       65      0       0.90        -1.0     state    u,v
+ ps        1      0       1.00        -1.0     state    prse
+ t        65      0       1.40        -1.0     state    tv
+ q        65      1       0.80        -1.0     state    q
+ sst       1      0       1.00        -1.0     state    sst
+ stl       1      0       1.00        -1.0     motley   sst
+ sti       1      0       1.00        -1.0     motley   sst
+ w        65     11       1.00        -1.0     state    w
+ qr       65     11       1.00        -1.0     state    qr
+ qs       65     11       1.00        -1.0     state    qs
+ qi       65     11       1.00        -1.0     state    qi
+ qg       65     11       1.00        -1.0     state    qg
+ ql       65     11       1.00        -1.0     state    qc
+ dbz      65     11       1.00        -1.0     state    dbz
+::

fix/gsi/gsiparm.anl.sh

@@ -149,6 +149,7 @@ OBS_INPUT::
    q_hyb_ens=${q_hyb_ens},
    aniso_a_en=.false.,generate_ens=.false.,
    n_ens=${nummem},
+   l_both_fv3sar_gfs_ens=${l_both_fv3sar_gfs_ens},n_ens_gfs=${nummem_gfs},n_ens_fv3sar=${nummem_fv3sar},
    beta_s0=${beta1_inv},s_ens_h=${ens_h},s_ens_v=${ens_v},
    regional_ensemble_option=${regional_ensemble_option},
    pseudo_hybens = .false.,
@@ -159,7 +160,18 @@ OBS_INPUT::
    jcap_ens=574,
    fv3sar_bg_opt=${fv3lam_bg_type},
    readin_localization=${readin_localization},
-   ens_fast_read=${ens_fast_read},
+   parallelization_over_ensmembers=${ens_fast_read},
+   nsclgrp=${nsclgrp},l_timloc_opt=.false.,ngvarloc=${ngvarloc},naensloc=${naensloc},
+   i_ensloccov4tim=${i_ensloccov4tim},i_ensloccov4var=${i_ensloccov4var},i_ensloccov4scl=${i_ensloccov4scl},
+   global_spectral_filter_sd=.false.,assign_vdl_nml=${assign_vdl_nml},vdl_scale=${vdl_scale},
+   vloc_varlist(1,1)='sf ',vloc_varlist(2,1)='w  ',vloc_varlist(3,1)='sf ',vloc_varlist(4,1)='w  ',
+   vloc_varlist(1,2)='vp ',vloc_varlist(2,2)='qr ',vloc_varlist(3,2)='vp ',vloc_varlist(4,2)='qr ',
+   vloc_varlist(1,3)='ps ',vloc_varlist(2,3)='qs ',vloc_varlist(3,3)='ps ',vloc_varlist(4,3)='qs ',
+   vloc_varlist(1,4)='t  ',vloc_varlist(2,4)='qi ',vloc_varlist(3,4)='t  ',vloc_varlist(4,4)='qi ',
+   vloc_varlist(1,5)='q  ',vloc_varlist(2,5)='qg ',vloc_varlist(3,5)='q  ',vloc_varlist(4,5)='qg ',
+   vloc_varlist(1,6)='sst',vloc_varlist(2,6)='ql ',vloc_varlist(3,6)='sst',vloc_varlist(4,6)='ql ',
+   vloc_varlist(1,7)='stl',vloc_varlist(2,7)='dbz',vloc_varlist(3,7)='stl',vloc_varlist(4,7)='dbz',
+   vloc_varlist(1,8)='sti',vloc_varlist(2,8)='aaa',vloc_varlist(3,8)='sti',vloc_varlist(4,8)='aaa',
  /
  &RAPIDREFRESH_CLDSURF
    dfi_radar_latent_heat_time_period=20.0,

scripts/exregional_nonvarcldanl.sh

@@ -269,7 +269,7 @@ fi
 if [ ${DO_ENKF_RADAR_REF} == "TRUE" ]; then
   l_qnr_from_qr=".true."
 fi
-if [ -r "${cycle_dir}/anal_radardbz_gsi${cycle_tag}/stdout" ]; then
+if [ -r "${cycle_dir}/anal_radardbz_gsi${cycle_tag}/gsi_complete_radar.txt" ] || [ -r "${cycle_dir}/anal_conv_dbz_gsi${cycle_tag}/gsi_complete_radar.txt" ]; then
   l_precip_clear_only=".true."
   l_qnr_from_qr=".true."
 fi

scripts/exregional_run_analysis.sh

@@ -188,6 +188,10 @@ else
   regional_ensemble_option=1
 fi
 
+if  [ ${ob_type} != "conv" ] || [ ${BKTYPE} -eq 1 ]; then #not using GDAS
+  l_both_fv3sar_gfs_ens=.false.
+fi
+
 #---------------------------------------------------------------------
 #
 # decide regional_ensemble_option: global ensemble (1) or FV3LAM ensemble (5)
@@ -246,10 +250,11 @@ if  [[ ${regional_ensemble_option:-1} -eq 5 ]]; then
   if [[ $ifound -ne ${NUM_ENS_MEMBERS} ]] || [[ ${BKTYPE} -eq 1 ]]; then
     print_info_msg "Not enough FV3_LAM ensembles, will fall to GDAS"
     regional_ensemble_option=1
+    l_both_fv3sar_gfs_ens=.false.
   fi
 fi
 #
-if  [[ ${regional_ensemble_option:-1} -eq 1 ]]; then #using GDAS
+if  [[ ${regional_ensemble_option:-1} -eq 1 || ${l_both_fv3sar_gfs_ens} = ".true." ]]; then #using GDAS
   #-----------------------------------------------------------------------
   # Make a list of the latest GFS EnKF ensemble
   #-----------------------------------------------------------------------
@@ -349,23 +354,32 @@ niter2=50
 lread_obs_save=.false.
 lread_obs_skip=.false.
 if_model_dbz=.false.
+nummem_gfs=0
+nummem_fv3sar=0
+anav_type=${ob_type}
 i_use_2mQ4B=2
 i_use_2mT4B=1
 
 # Determine if hybrid option is available
 memname='atmf009'
 
 if [ ${regional_ensemble_option:-1} -eq 5 ]  && [ ${BKTYPE} != 1  ]; then 
-  nummem=$NUM_ENS_MEMBERS
+  if [ ${l_both_fv3sar_gfs_ens} = ".true." ]; then
+    nummem_gfs=$(more filelist03 | wc -l)
+    nummem_gfs=$((nummem_gfs - 3 ))
+  fi
+  nummem_fv3sar=$NUM_ENS_MEMBERS
+  nummem=`expr ${nummem_gfs} + ${nummem_fv3sar}`
   print_info_msg "$VERBOSE" "Do hybrid with FV3LAM ensemble"
   ifhyb=.true.
   print_info_msg "$VERBOSE" " Cycle ${YYYYMMDDHH}: GSI hybrid uses FV3LAM ensemble with n_ens=${nummem}" 
   echo " ${YYYYMMDDHH}(${cycle_type}): GSI hybrid uses FV3LAM ensemble with n_ens=${nummem}" >> ${EXPTDIR}/log.cycles
   grid_ratio_ens="1"
   ens_fast_read=.true.
 else    
-  nummem=$(more filelist03 | wc -l)
-  nummem=$((nummem - 3 ))
+  nummem_gfs=$(more filelist03 | wc -l)
+  nummem_gfs=$((nummem_gfs - 3 ))
+  nummem=${nummem_gfs}
   if [[ ${nummem} -ge ${HYBENSMEM_NMIN} ]]; then
     print_info_msg "$VERBOSE" "Do hybrid with ${memname}"
     ifhyb=.true.
@@ -376,6 +390,9 @@ else
     print_info_msg "$VERBOSE" " Hybrid needs at least ${HYBENSMEM_NMIN} ${memname} ensembles, only ${nummem} available"
     echo " ${YYYYMMDDHH}(${cycle_type}): GSI dose pure 3DVAR" >> ${EXPTDIR}/log.cycles
   fi
+  if [[ ${anav_type} == "conv_dbz" ]]; then
+    anav_type="conv"
+  fi
 fi
 
 #
@@ -473,7 +490,7 @@ else
   esac
 fi
 
-if [[ ${gsi_type} == "OBSERVER" || ${ob_type} == "conv" ]]; then
+if [[ ${gsi_type} == "OBSERVER" || ${anav_type} == "conv" || ${anav_type} == "conv_dbz" ]]; then
 
   obs_files_source[0]=${obspath_tmp}/${obsfileprefix}.t${HH}${SUBH}z.prepbufr.tm00
   obs_files_target[0]=prepbufr
@@ -486,15 +503,19 @@ if [[ ${gsi_type} == "OBSERVER" || ${ob_type} == "conv" ]]; then
   obs_files_source[${obs_number}]=${obspath_tmp}/${obsfileprefix}.t${HH}${SUBH}z.nexrad.tm00.bufr_d
   obs_files_target[${obs_number}]=l2rwbufr
 
-  if [ ${DO_ENKF_RADAR_REF} == "TRUE" ]; then
+  if [[ ${DO_ENKF_RADAR_REF} == "TRUE" || ${anav_type} == "conv_dbz" ]]; then
     obs_number=${#obs_files_source[@]}
-    obs_files_source[${obs_number}]=${cycle_dir}/process_radarref/00/Gridded_ref.nc
+    if [ ${cycle_type} == "spinup" ]; then
+      obs_files_source[${obs_number}]=${cycle_dir}/process_radarref_spinup/00/Gridded_ref.nc
+    else
+      obs_files_source[${obs_number}]=${cycle_dir}/process_radarref/00/Gridded_ref.nc
+    fi
     obs_files_target[${obs_number}]=dbzobs.nc
   fi
 
 else
 
-  if [ ${ob_type} == "radardbz" ]; then
+  if [ ${anav_type} == "radardbz" ]; then
 
     if [ ${cycle_type} == "spinup" ]; then
       obs_files_source[0]=${cycle_dir}/process_radarref_spinup/00/Gridded_ref.nc
@@ -506,7 +527,7 @@ else
   fi
 
 
-  if [ ${ob_type} == "AERO" ]; then
+  if [ ${anav_type} == "AERO" ]; then
     obs_files_source[0]=${OBSPATH_PM}/${YYYYMMDD}/pm25.airnow.${YYYYMMDD}${HH}.bufr
     obs_files_target[0]=pm25bufr
   fi
@@ -624,27 +645,42 @@ if [ ${DO_ENKF_RADAR_REF} == "TRUE" ]; then
   beta1_inv=0.0
   if_model_dbz=.true.
 fi
-if [[ ${gsi_type} == "ANALYSIS" && ${ob_type} == "radardbz" ]]; then
+if [[ ${gsi_type} == "ANALYSIS" && ${anav_type} == "radardbz" ]]; then
   ANAVINFO=${FIX_GSI}/${ENKF_ANAVINFO_DBZ_FN}
   miter=1
   niter1=100
   niter2=0
   bkgerr_vs=0.1
   bkgerr_hzscl="0.4,0.5,0.6"
   beta1_inv=0.0
-  ens_h=4.10790
-  ens_v=-0.30125
   readin_localization=.false.
+  ens_h=${ens_h_radardbz}
+  ens_v=${ens_v_radardbz}
+  nsclgrp=1
+  ngvarloc=1
+  i_ensloccov4tim=0
+  i_ensloccov4var=0
+  i_ensloccov4scl=0
   q_hyb_ens=.true.
   if_model_dbz=.true.
 fi
+if [[ ${gsi_type} == "ANALYSIS" && ${anav_type} == "conv_dbz" ]]; then
+  ANAVINFO=${FIX_GSI}/${ANAVINFO_CONV_DBZ_FN}
+  beta1_inv=0.0
+  if_model_dbz=.true.
+fi
+naensloc=`expr ${nsclgrp} \* ${ngvarloc} + ${nsclgrp} - 1`
+if [ ${assign_vdl_nml} = ".true." ]; then
+  nsclgrp=`expr ${nsclgrp} \* ${ngvarloc}`
+  ngvarloc=1
+fi
 CONVINFO=${FIX_GSI}/${CONVINFO_FN}
 HYBENSINFO=${FIX_GSI}/${HYBENSINFO_FN}
 OBERROR=${FIX_GSI}/${OBERROR_FN}
 BERROR=${FIX_GSI}/${BERROR_FN}
 
 
-if [[ ${gsi_type} == "ANALYSIS" && ${ob_type} == "AERO" ]]; then
+if [[ ${gsi_type} == "ANALYSIS" && ${anav_type} == "AERO" ]]; then
   if [ ${BKTYPE} -eq 1 ]; then
     echo "cold start, skip GSI SD DA"
     exit 0
@@ -833,7 +869,7 @@ fi
 # skip radar reflectivity analysis if no RRFSE ensemble
 #-----------------------------------------------------------------------
 
-if [[ ${gsi_type} == "ANALYSIS" && ${ob_type} == "radardbz" ]]; then
+if [[ ${gsi_type} == "ANALYSIS" && ${anav_type} == "radardbz" ]]; then
   if  [[ ${regional_ensemble_option:-1} -eq 1 ]]; then
      echo "No RRFSE ensemble available, cannot do radar reflectivity analysis"
      exit 0
@@ -879,7 +915,7 @@ EOF
 #
 gsi_exec="${EXECDIR}/gsi.x"
 
-if [[ ${gsi_type} == "ANALYSIS" && ${ob_type} == "AERO" ]]; then
+if [[ ${gsi_type} == "ANALYSIS" && ${anav_type} == "AERO" ]]; then
   gsi_exec="${EXECDIR}/gsi.x.sd"
 fi
 
@@ -916,7 +952,7 @@ fi
 $APRUN ./gsi.x < gsiparm.anl > stdout 2>&1 || print_err_msg_exit "\
 Call to executable to run GSI returned with nonzero exit code."
 
-if [ ${ob_type} == "radardbz" ]; then
+if [ ${anav_type} == "radardbz" ]; then
   cat fort.238 > $comout/rrfs_a.t${HH}z.fits3.tm00
 else
   mv fort.207 fit_rad1
@@ -935,7 +971,7 @@ else
 
 fi
 
-cp stdout $comout/stdout_${ob_type}
+cp stdout $comout/stdout_${anav_type}
 #
 #-----------------------------------------------------------------------
 #
@@ -945,6 +981,9 @@ cp stdout $comout/stdout_${ob_type}
 #-----------------------------------------------------------------------
 #
 touch gsi_complete.txt
+if [[ ${anav_type} == "radardbz" || ${anav_type} == "conv_dbz" ]]; then
+  touch gsi_complete_radar.txt # for nonvarcldanl
+fi
 #
 #-----------------------------------------------------------------------
 #

scripts/exregional_run_enkf.sh

@@ -238,14 +238,12 @@ CONVINFO=${FIX_GSI}/convinfo.rrfs
 
 if [ ${ob_type} == "conv" ]; then
   ANAVINFO=${FIX_GSI}/${ENKF_ANAVINFO_FN}
-  CORRLENGTH=400
-  LNSIGCUTOFF=0.5
   found_ob_type=1
 fi
 if [ ${ob_type} == "radardbz" ]; then
   ANAVINFO=${FIX_GSI}/${ENKF_ANAVINFO_DBZ_FN}
-  CORRLENGTH=18
-  LNSIGCUTOFF=0.5
+  CORRLENGTH=${CORRLENGTH_radardbz}
+  LNSIGCUTOFF=${LNSIGCUTOFF_radardbz}
   found_ob_type=1
 fi
 if [ ${found_ob_type} == 0 ]; then

ush/config.sh.RRFS_CONUS_13km

@@ -8,6 +8,7 @@ EXPT_SUBDIR="rrfs_conus_13km"
 
 PREDEF_GRID_NAME=RRFS_CONUS_13km
 . set_rrfs_config_general.sh
+. set_rrfs_config_SDL_VDL_MixEn.sh
 
 DO_DACYCLE="TRUE"
 DO_SURFACE_CYCLE="TRUE"

ush/config.sh.RRFS_CONUS_3km

@@ -8,6 +8,7 @@ EXPT_SUBDIR="RRFS_CONUS_3km"
 
 PREDEF_GRID_NAME=RRFS_CONUS_3km
 . set_rrfs_config_general.sh
+. set_rrfs_config_SDL_VDL_MixEn.sh
 
 DO_DACYCLE="TRUE"
 DO_SURFACE_CYCLE="TRUE"

ush/config.sh.RRFS_CONUS_3km.2022HWT

@@ -7,6 +7,7 @@ EXPT_SUBDIR="RRFS_CONUS_3km"
 
 PREDEF_GRID_NAME=RRFS_CONUS_3km
 . set_rrfs_config_general.sh
+. set_rrfs_config_SDL_VDL_MixEn.sh
 
 DO_DACYCLE="TRUE"
 DO_SURFACE_CYCLE="TRUE"

ush/config.sh.RRFS_NA_13km

@@ -6,6 +6,7 @@ EXPT_SUBDIR="RRFS_NA_13km_dev1"
 
 PREDEF_GRID_NAME=GSD_RAP13km
 . set_rrfs_config_general.sh
+. set_rrfs_config_SDL_VDL_MixEn.sh
 
 ADDNL_OUTPUT_GRIDS=( "130" "242" )
 

ush/config.sh.RRFS_NA_3km

@@ -7,6 +7,7 @@ EXPT_SUBDIR="RRFS_NA_3km"
 
 PREDEF_GRID_NAME="RRFS_NA_3km"
 . set_rrfs_config_general.sh
+. set_rrfs_config_SDL_VDL_MixEn.sh
 
 DO_DACYCLE="TRUE"
 DO_SURFACE_CYCLE="TRUE"

ush/config.sh.rrfs_a_c13

@@ -9,6 +9,7 @@ EXPT_SUBDIR="rrfs_conus_13km"
 PREDEF_GRID_NAME=RRFS_CONUS_13km
 
 . set_rrfs_config_general.sh
+. set_rrfs_config_SDL_VDL_MixEn.sh
 
 DO_DACYCLE="TRUE"
 DO_SURFACE_CYCLE="TRUE"

ush/config.sh.rrfs_a_c3

@@ -8,6 +8,7 @@ EXPT_SUBDIR="rrfs_conus"
 PREDEF_GRID_NAME=RRFS_CONUS_3km
 
 . set_rrfs_config_general.sh
+. set_rrfs_config_SDL_VDL_MixEn.sh
 
 #DO_ENSEMBLE="TRUE"
 #DO_ENSFCST="TRUE"