Add ability to use GEFS replay ICs

parm/config/gefs/config.base

@@ -262,7 +262,9 @@ export MEMDIR="mem${ENSMEM}"
 
 # initialize ocean ensemble members with perturbations
 #   if true, only occurs for members greater than zero
+export USE_REPLAY_ICS=@USE_REPLAY_ICS@
 export USE_OCN_PERTURB_FILES=@USE_OCN_PERTURB_FILES@
+export USE_ATM_PERTURB_FILES=@USE_ATM_PERTURB_FILES@
 
 export DOIAU="NO"  # While we are not doing IAU, we may want to warm start w/ IAU in the future
 # Check if cycle is cold starting

parm/config/gefs/config.fcst

@@ -256,6 +256,9 @@ else
   export io_layout="1,1"
 fi
 
+if [[ ${USE_REPLAY_ICS} == "true" ]]; then
+    export reforecast="YES"
+fi
 # Remember config.efcs will over-ride these values for ensemble forecasts
 # if these variables are re-defined there.
 # Otherwise, the ensemble forecast will inherit from config.fcst

parm/config/gefs/config.stage_ic

@@ -13,12 +13,14 @@ case "${CASE}" in
     export CPL_ICEIC=""
     export CPL_OCNIC=""
     export CPL_WAVIC=""
+    export CPL_MEDIC=""
     ;;
   "C48")
     export CPL_ATMIC="gefs_test"
     export CPL_ICEIC="gefs_test"
     export CPL_OCNIC="gefs_test"
     export CPL_WAVIC="gefs_test"
+    export CPL_MEDIC="gefs_test"
     ;;
   *)
     echo "FATAL ERROR Unrecognized resolution: ${CASE}"

parm/config/gefs/config.wave

@@ -81,7 +81,11 @@ export WAVEWND_FID=
 # The start time reflects the number of hindcast hours prior to the cycle initial time
 export FHMAX_WAV=${FHMAX_GFS}
 export WAVHINDH=0
-export FHMIN_WAV=0
+if [[ "${USE_REPLAY_ICS:-"false"}" == "true" ]]; then
+    export FHMIN_WAV=3
+else
+    export FHMIN_WAV=0
+fi
 export FHOUT_WAV=3
 export FHMAX_HF_WAV=120
 export FHOUT_HF_WAV=1

parm/config/gefs/yaml/defaults.yaml

@@ -6,5 +6,7 @@ base:
   DO_MERGENSST: "NO"
   KEEPDATA: "NO" 
   FHMAX_GFS: 120
+  USE_REPLAY_ICS: "false"
+  USE_ATM_PERTURB_FILES: "false"
   USE_OCN_PERTURB_FILES: "false"
 

scripts/exgfs_wave_post_pnt.sh

@@ -247,8 +247,9 @@ source "${USHgfs}/preamble.sh"
         -e "s/FORMAT/F/g" \
                                ww3_outp_spec.inp.tmpl > ww3_outp.inp
 
-    ${NLN} mod_def.$waveuoutpGRD mod_def.ww3
-    HMS="${cyc}0000"
+    ln -s mod_def.$waveuoutpGRD mod_def.ww3
+    HH=$(date --utc -d "${PDY:0:8} ${cyc} + ${FHMIN_WAV} hours" +%H)
+    HMS="${HH}0000"
     if [[ -f "${COM_WAVE_HISTORY}/${WAV_MOD_TAG}.out_pnt.${waveuoutpGRD}.${PDY}.${HMS}" ]]; then
       ${NLN} "${COM_WAVE_HISTORY}/${WAV_MOD_TAG}.out_pnt.${waveuoutpGRD}.${PDY}.${HMS}" \
         "./out_pnt.${waveuoutpGRD}"

scripts/exglobal_stage_ic.sh

@@ -33,14 +33,25 @@ for MEMDIR in "${MEMDIR_ARRAY[@]}"; do
     # Stage the FV3 restarts to ROTDIR (warm start)
     RUN=${rCDUMP} YMD=${gPDY} HH=${gcyc} declare_from_tmpl COM_ATMOS_RESTART_PREV:COM_ATMOS_RESTART_TMPL
     [[ ! -d "${COM_ATMOS_RESTART_PREV}" ]] && mkdir -p "${COM_ATMOS_RESTART_PREV}"
-    for ftype in coupler.res fv_core.res.nc; do
-      src="${BASE_CPLIC}/${CPL_ATMIC:-}/${PDY}${cyc}/${MEMDIR}/atmos/${PDY}.${cyc}0000.${ftype}"
-      tgt="${COM_ATMOS_RESTART_PREV}/${PDY}.${cyc}0000.${ftype}"
-      ${NCP} "${src}" "${tgt}"
-      rc=$?
-      ((rc != 0)) && error_message "${src}" "${tgt}" "${rc}"
-      err=$((err + rc))
-    done
+
+    # fv_core.res.nc file
+    src="${BASE_CPLIC}/${CPL_ATMIC:-}/${PDY}${cyc}/${MEMDIR}/atmos/${PDY}.${cyc}0000.fv_core.res.nc"
+    tgt="${COM_ATMOS_RESTART_PREV}/${PDY}.${cyc}0000.fv_core.res.nc"
+    ${NCP} "${src}" "${tgt}"
+    rc=$?
+    ((rc != 0)) && error_message "${src}" "${tgt}" "${rc}"
+    err=$((err + rc))
+
+    # coupler.res file
+    if [[ "${USE_REPLAY_ICS:-"false"}" == "false" ]]; then
+        src="${BASE_CPLIC}/${CPL_ATMIC:-}/${PDY}${cyc}/${MEMDIR}/atmos/${PDY}.${cyc}0000.coupler.res"
+        tgt="${COM_ATMOS_RESTART_PREV}/${PDY}.${cyc}0000.coupler.res"
+        ${NCP} "${src}" "${tgt}"
+        rc=$?
+        ((rc != 0)) && error_message "${src}" "${tgt}" "${rc}"
+        err=$((err + rc))
+    fi
+
     for ftype in ca_data fv_core.res fv_srf_wnd.res fv_tracer.res phy_data sfc_data; do
       for ((tt = 1; tt <= 6; tt++)); do
         src="${BASE_CPLIC}/${CPL_ATMIC:-}/${PDY}${cyc}/${MEMDIR}/atmos/${PDY}.${cyc}0000.${ftype}.tile${tt}.nc"
@@ -72,6 +83,17 @@ for MEMDIR in "${MEMDIR_ARRAY[@]}"; do
       done
     done
   fi
+
+  # Atmosphere Perturbation Files (usually used with replay ICS)
+  # Extra zero on MEMDIR ensure we have a number even if the string is empty
+  if (( 0${MEMDIR:3} > 0 )) && [[ "${USE_ATM_PERTURB_FILES:-false}" == "true" ]]; then
+      src="${BASE_CPLIC}/${CPL_OCNIC:-}/${PDY}${cyc}/${MEMDIR}/atmos/${PDY}.${cyc}0000.fv3_perturbation.nc"
+      tgt="${COM_ATMOS_RESTART_PREV}/${PDY}.${cyc}0000.fv3_perturbation.nc"
+      ${NCP} "${src}" "${tgt}"
+      rc=${?}
+      ((rc != 0)) && error_message "${src}" "${tgt}" "${rc}"
+      err=$((err + rc))
+  fi
 
   # Stage ocean initial conditions to ROTDIR (warm start)
   if [[ "${DO_OCN:-}" = "YES" ]]; then
@@ -107,8 +129,8 @@ for MEMDIR in "${MEMDIR_ARRAY[@]}"; do
     # Ocean Perturbation Files
     # Extra zero on MEMDIR ensure we have a number even if the string is empty
     if (( 0${MEMDIR:3} > 0 )) && [[ "${USE_OCN_PERTURB_FILES:-false}" == "true" ]]; then
-        src="${BASE_CPLIC}/${CPL_OCNIC:-}/${PDY}${cyc}/${MEMDIR}/ocean/${PDY}.${cyc}0000.mom6_increment.nc"
-        tgt="${COM_OCEAN_RESTART_PREV}/${PDY}.${cyc}0000.mom6_increment.nc"
+        src="${BASE_CPLIC}/${CPL_OCNIC:-}/${PDY}${cyc}/${MEMDIR}/ocean/${PDY}.${cyc}0000.mom6_perturbation.nc"
+        tgt="${COM_OCEAN_RESTART_PREV}/${PDY}.${cyc}0000.mom6_perturbation.nc"
         ${NCP} "${src}" "${tgt}"
         rc=${?}
         ((rc != 0)) && error_message "${src}" "${tgt}" "${rc}"
@@ -149,7 +171,7 @@ for MEMDIR in "${MEMDIR_ARRAY[@]}"; do
 
   # Stage the WW3 initial conditions to ROTDIR (warm start; TODO: these should be placed in $RUN.$gPDY/$gcyc)
   if [[ "${DO_WAVE:-}" = "YES" ]]; then
-    YMD=${PDY} HH=${cyc} declare_from_tmpl COM_WAVE_RESTART
+    YMD=${gPDY} HH=${gcyc} declare_from_tmpl COM_WAVE_RESTART
     [[ ! -d "${COM_WAVE_RESTART}" ]] && mkdir -p "${COM_WAVE_RESTART}"
     for grdID in ${waveGRD}; do # TODO: check if this is a bash array; if so adjust
       src="${BASE_CPLIC}/${CPL_WAVIC:-}/${PDY}${cyc}/${MEMDIR}/wave/${PDY}.${cyc}0000.restart.${grdID}"

ush/forecast_det.sh

@@ -6,7 +6,10 @@ UFS_det(){
   echo "SUB ${FUNCNAME[0]}: Run type determination for UFS"
 
   # Determine if the current cycle is a warm start (based on the availability of restarts)
-  if [[ "${DOIAU:-}" == "YES" ]]; then
+  USE_REPLAY_ICS=${USE_REPLAY_ICS:-'false'}
+  if [[ "${USE_REPLAY_ICS}" == "true" ]]; then
+      warm_start=".true."
+  elif [[ "${DOIAU:-}" == "YES" ]]; then
     if [[ -f "${COM_ATMOS_RESTART_PREV}/${current_cycle_begin:0:8}.${current_cycle_begin:8:2}0000.coupler.res" ]]; then
       warm_start=".true."
     fi

ush/forecast_postdet.sh

@@ -33,7 +33,10 @@ FV3_postdet() {
 
     # Create an array of FV3 restart files
     local fv3_restart_files tile_files fv3_restart_file restart_file
-    fv3_restart_files=(coupler.res fv_core.res.nc)
+    fv3_restart_files=(fv_core.res.nc)
+    if [[ "${USE_REPLAY_ICS}" == "false" ]]; then
+        fv_restart_files+=(coupler.res)
+    fi
     tile_files=(fv_core.res fv_srf_wnd.res fv_tracer.res phy_data sfc_data ca_data)
     local nn tt
     for (( nn = 1; nn <= ntiles; nn++ )); do
@@ -97,16 +100,28 @@ FV3_postdet() {
         local model_start_time="${current_cycle}"
         local model_current_time="${current_cycle}"
       fi
+      if [[ "${USE_REPLAY_ICS}" != "true" ]]; then
       rm -f "${DATA}/INPUT/coupler.res"
       cat >> "${DATA}/INPUT/coupler.res" << EOF
       3        (Calendar: no_calendar=0, thirty_day_months=1, julian=2, gregorian=3, noleap=4)
       ${model_start_time:0:4}  ${model_start_time:4:2}  ${model_start_time:6:2}  ${model_start_time:8:2}  0  0        Model start time: year, month, day, hour, minute, second
       ${model_current_time:0:4}  ${model_current_time:4:2}  ${model_current_time:6:2}  ${model_current_time:8:2}  0  0        Current model time: year, month, day, hour, minute, second
 EOF
+      fi
 
       # Create a array of increment files
       local inc_files inc_file iaufhrs iaufhr
-      if [[ "${DOIAU}" == "YES" ]]; then
+      if [[ "${USE_REPLAY_ICS}" == "true" ]]; then
+        if (( MEMBER > 0 )) && [[ "${USE_ATM_PERTURB_FILES:-false}" == "true" ]]; then
+            increment_file="${COM_ATMOS_RESTART_PREV}/${current_cycle:0:8}.${current_cycle:8:2}0000.fv3_perturbation.nc"
+            ${NCP} "${increment_file}" "${DATA}/INPUT/fv3_perturbation.nc" \
+            || ( echo "FATAL ERROR: Unable to copy FV3 perturbation file, ABORT!"; exit 1 )
+            read_increment=".true."
+            res_latlon_dynamics="fv3_perturbation.nc"
+        fi
+      elif [[ "${RUN}" == "gefs" ]]; then
+        inc_files=()
+      elif [[ "${DOIAU}" == "YES" ]]; then
         # create an array of inc_files for each IAU hour
         IFS=',' read -ra iaufhrs <<< "${IAUFHRS}"
         inc_files=()
@@ -407,7 +422,7 @@ MOM6_postdet() {
     # TODO if [[ $RUN} == "gefs" ]] block maybe be needed
     #     to ensure it does not interfere with the GFS when ensemble is updated in the GFS
     if (( MEMBER > 0 )) && [[ "${ODA_INCUPD:-False}" == "True" ]]; then
-      ${NCP} "${COM_OCEAN_RESTART_PREV}/${restart_date:0:8}.${restart_date:0:8}0000.mom6_increment.nc" "${DATA}/INPUT/mom6_increment.nc" \
+      ${NCP} "${COM_OCEAN_RESTART_PREV}/${restart_date:0:8}.${restart_date:8:2}0000.mom6_perturbation.nc" "${DATA}/INPUT/mom6_increment.nc" \
       || ( echo "FATAL ERROR: Unable to copy ensemble MOM6 increment, ABORT!"; exit 1 )
     fi
   fi  # if [[ "${RERUN}" == "NO" ]]; then
@@ -490,6 +505,7 @@ MOM6_out() {
   esac
 
   # Copy MOM6 restarts at the end of the forecast segment to COM for RUN=gfs|gefs
+  if [[ "${USE_REPLAY_ICS}" != "true" ]]; then
   local restart_file
   if [[ "${RUN}" == "gfs" || "${RUN}" == "gefs" ]]; then
     echo "Copying MOM6 restarts for 'RUN=${RUN}' at ${forecast_end_cycle}"
@@ -499,6 +515,7 @@ MOM6_out() {
              "${COM_OCEAN_RESTART}/${restart_file}"
     done
   fi
+  fi
 
   # Copy restarts at the beginning/middle of the next assimilation cycle to COM for RUN=gdas|enkfgdas|enkfgfs
   if [[ "${RUN}" =~ "gdas" || "${RUN}" == "enkfgfs" ]]; then
@@ -590,6 +607,7 @@ CICE_out() {
   ${NCP} "${DATA}/ice_in" "${COM_CONF}/ufs.ice_in"
 
   # Copy CICE restarts at the end of the forecast segment to COM for RUN=gfs|gefs
+  if [[ "${USE_REPLAY_ICS}" != "true" ]]; then
   local seconds source_file target_file
   if [[ "${RUN}" == "gfs" || "${RUN}" == "gefs" ]]; then
     echo "Copying CICE restarts for 'RUN=${RUN}' at ${forecast_end_cycle}"
@@ -599,6 +617,7 @@ CICE_out() {
     ${NCP} "${DATArestart}/CICE_RESTART/${source_file}" \
            "${COM_ICE_RESTART}/${target_file}"
   fi
+  fi
 
   # Copy restarts at the beginning/middle of the next assimilation cycle to COM for RUN=gdas|enkfgdas|enkfgfs
   if [[ "${RUN}" =~ "gdas" || "${RUN}" == "enkfgfs" ]]; then
@@ -729,6 +748,7 @@ CMEPS_out() {
   echo "SUB ${FUNCNAME[0]}: Copying output data for CMEPS mediator"
 
   # Copy mediator restarts at the end of the forecast segment to COM for RUN=gfs|gefs
+  if [[ "${USE_REPLAY_ICS}" != "true" ]]; then
   echo "Copying mediator restarts for 'RUN=${RUN}' at ${forecast_end_cycle}"
   local seconds source_file target_file
   seconds=$(to_seconds "${forecast_end_cycle:8:2}"0000)
@@ -740,6 +760,7 @@ CMEPS_out() {
   else
     echo "Mediator restart '${DATArestart}/CMEPS_RESTART/${source_file}' not found."
   fi
+  fi
 
   # Copy restarts at the beginning/middle of the next assimilation cycle to COM for RUN=gdas|enkfgdas|enkfgfs
   if [[ "${RUN}" =~ "gdas" || "${RUN}" == "enkfgfs" ]]; then

ush/parsing_model_configure_FV3.sh

@@ -22,7 +22,11 @@ local SMONTH=${model_start_date:4:2}
 local SDAY=${model_start_date:6:2}
 local SHOUR=${model_start_date:8:2}
 # FHMAX
-local FHROT=${IAU_FHROT:-0}
+if [[ "${USE_REPLAY_ICS:-"false"}" == "true" ]]; then
+    local FHROT=3
+else
+    local FHROT=${IAU_FHROT:-0}
+fi
 local DT_ATMOS=${DELTIM}
 local RESTART_INTERVAL="${FV3_RESTART_FH[*]}"
 # QUILTING

ush/parsing_namelists_WW3.sh

@@ -18,7 +18,7 @@ WW3_namelists(){
   fi
   # Set time stamps for model start and output
   # For special case when IAU is on but this is an initial half cycle 
-  if [ $IAU_OFFSET = 0 ]; then
+  if [ ${IAU_OFFSET:-0} = 0 ]; then
     ymdh_beg=$YMDH
   else
     ymdh_beg=$($NDATE -$WAVHINDH $YMDH)

workflow/rocoto/gefs_tasks.py

@@ -9,21 +9,31 @@ def __init__(self, app_config: AppConfig, cdump: str) -> None:
         super().__init__(app_config, cdump)
 
     def stage_ic(self):
-
-        cpl_ic = self._configs['stage_ic']
-
-        deps = []
-
+        cpl_ic = self._configs['stage_ic'] 
+        deps = [] 
         # Atm ICs
         if self.app_config.do_atm:
             prefix = f"{cpl_ic['BASE_CPLIC']}/{cpl_ic['CPL_ATMIC']}/@Y@m@d@H/mem000/atmos"
-            for file in ['gfs_ctrl.nc'] + \
-                        [f'{datatype}_data.tile{tile}.nc'
-                         for datatype in ['gfs', 'sfc']
-                         for tile in range(1, self.n_tiles + 1)]:
-                data = f"{prefix}/{file}"
+            if self._base['EXP_WARM_START']:
+                for file in ['fv_core.res.nc'] + \
+                            [f'{datatype}.tile{tile}.nc'
+                             for datatype in ['ca_data', 'fv_core.res', 'fv_srf_wnd.res', 'fv_tracer.res', 'phy_data', 'sfc_data']
+                             for tile in range(1, self.n_tiles + 1)]:
+                    data = f"{prefix}/@Y@m@d.@H0000.{file}"
+                    dep_dict = {'type': 'data', 'data': data}
+                    deps.append(rocoto.add_dependency(dep_dict))
+                prefix = f"{cpl_ic['BASE_CPLIC']}/{cpl_ic['CPL_ATMIC']}/@Y@m@d@H/mem000/med"
+                data = f"{prefix}/@Y@m@d.@H0000.ufs.cpld.cpl.r.nc"
                 dep_dict = {'type': 'data', 'data': data}
                 deps.append(rocoto.add_dependency(dep_dict))
+            else:
+                for file in ['gfs_ctrl.nc'] + \
+                            [f'{datatype}_data.tile{tile}.nc'
+                             for datatype in ['gfs', 'sfc']
+                             for tile in range(1, self.n_tiles + 1)]:
+                    data = f"{prefix}/{file}"
+                    dep_dict = {'type': 'data', 'data': data}
+                    deps.append(rocoto.add_dependency(dep_dict))
 
         # Ocean ICs
         if self.app_config.do_ocean: