diff --git a/jobs/JGLOBAL_EXTRACTVARS b/jobs/JGLOBAL_EXTRACTVARS
new file mode 100755
index 0000000000..3478ca3976
--- /dev/null
+++ b/jobs/JGLOBAL_EXTRACTVARS
@@ -0,0 +1,47 @@
+#! /usr/bin/env bash
+
+source "${HOMEgfs}/ush/preamble.sh"
+source "${HOMEgfs}/ush/jjob_header.sh" -e "extractvars" -c "base extractvars"
+
+# Set COM Paths
+for grid in '0p25' '0p50' '1p00'; do
+  prod_dir="COMIN_ATMOS_GRIB_${grid}"
+  GRID=${grid} YMD=${PDY} HH=${cyc} declare_from_tmpl -rx "${prod_dir}:COM_ATMOS_GRIB_GRID_TMPL"
+  if [[ ! -d "${!prod_dir}" ]]; then mkdir -p "${!prod_dir}"; fi
+done
+
+YMD="${PDY}" HH="${cyc}" declare_from_tmpl -rx \
+    "COMIN_OCEAN_HISTORY:COM_OCEAN_HISTORY_TMPL" \
+    "COMIN_OCEAN_GRIB:COM_OCEAN_GRIB_TMPL" \
+    "COMIN_OCEAN_NETCDF:COM_OCEAN_NETCDF_TMPL" \
+    "COMIN_ICE_HISTORY:COM_ICE_HISTORY_TMPL" \
+    "COMIN_ICE_GRIB:COM_ICE_GRIB_TMPL" \
+    "COMIN_ICE_NETCDF:COM_ICE_NETCDF_TMPL" \
+    "COMIN_WAVE_GRID:COM_WAVE_GRID_TMPL"
+
+if [[ "${DO_ATM}" == "YES" ]]; then
+  if [[ ! -d "${ARC_RFCST_PROD_ATMOS_F2D}" ]]; then mkdir -p "${ARC_RFCST_PROD_ATMOS_F2D}"; fi
+  if [[ ! -d "${ARC_RFCST_PROD_ATMOS_F3D}" ]]; then mkdir -p "${ARC_RFCST_PROD_ATMOS_F3D}"; fi
+fi
+if [[ "${DO_OCN}" == "YES" ]]; then
+  if [[ ! -d "${ARC_RFCST_PROD_OCN}" ]]; then mkdir -p "${ARC_RFCST_PROD_OCN}"; fi
+fi
+if [[ "${DO_ICE}" == "YES" ]]; then
+  if [[ ! -d "${ARC_RFCST_PROD_ICE}" ]]; then mkdir -p "${ARC_RFCST_PROD_ICE}"; fi
+fi
+if [[ "${DO_WAVE}" == "YES" ]]; then
+  if [[ ! -d "${ARC_RFCST_PROD_WAV}" ]]; then mkdir -p "${ARC_RFCST_PROD_WAV}"; fi
+fi
+
+# Execute the Script
+"${SCRgfs}/exglobal_extractvars.sh"
+status=$?
+(( status != 0 )) && exit "${status}"
+
+##########################################
+# Remove the Temporary working directory
+##########################################
+cd "${DATAROOT}" || true
+[[ "${KEEPDATA}" = "NO" ]] && rm -rf "${DATA}"
+
+exit 0
diff --git a/jobs/rocoto/extractvars.sh b/jobs/rocoto/extractvars.sh
new file mode 100755
index 0000000000..a872431358
--- /dev/null
+++ b/jobs/rocoto/extractvars.sh
@@ -0,0 +1,23 @@
+#! /usr/bin/env bash
+
+source "${HOMEgfs}/ush/preamble.sh"
+
+###############################################################
+echo
+echo "=============== START TO SOURCE FV3GFS WORKFLOW MODULES ==============="
+. "${HOMEgfs}/ush/load_fv3gfs_modules.sh"
+status=$?
+[[ "${status}" -ne 0 ]] && exit "${status}"
+
+export job="extractvars"
+export jobid="${job}.$$"
+
+###############################################################
+echo
+echo "=============== START TO RUN EXTRACTVARS ==============="
+# Execute the JJOB
+"${HOMEgfs}/jobs/JGLOBAL_EXTRACTVARS"
+status=$?
+[[ "${status}" -ne 0 ]] && exit "${status}"
+
+exit 0
diff --git a/parm/config/gefs/config.base b/parm/config/gefs/config.base
index 9808b96579..1f5c9228e2 100644
--- a/parm/config/gefs/config.base
+++ b/parm/config/gefs/config.base
@@ -137,6 +137,7 @@ export DO_WAVE="NO"
 export DO_OCN="NO"
 export DO_ICE="NO"
 export DO_AERO="NO"
+export DO_EXTRACTVARS="@DO_EXTRACTVARS@" # Option to process and extract a subset of products to save on disk
 export AERO_FCST_CDUMP="" # When to run aerosol forecast: gdas, gfs, or both
 export AERO_ANL_CDUMP="" # When to run aerosol analysis: gdas, gfs, or both
 export WAVE_CDUMP="" # When to include wave suite: gdas, gfs, or both
diff --git a/parm/config/gefs/config.extractvars b/parm/config/gefs/config.extractvars
new file mode 100644
index 0000000000..706fe18450
--- /dev/null
+++ b/parm/config/gefs/config.extractvars
@@ -0,0 +1,41 @@
+#! /usr/bin/env bash                                                                                                                                                                                                                     
+
+########## config.extractvars ##########
+# Extractvars specific
+
+echo "BEGIN: config.extractvars"
+
+. "${EXPDIR}/config.resources" extractvars
+
+export COMPRSCMD=${COMPRSCMD:-bzip2}
+
+export compress_ocn=0 #1: Compress extracted ocean product, 0: Do not compress extracted ocean product
+export compress_ice=0 #1: Compress extracted ice product, 0: Do not compress extracted ice product
+
+export ocnres="5p00" # Resolution of ocean products
+export iceres="5p00" # Resolution of ice products
+export wavres="5p00" # Resolution of wave products
+
+export depthvar_name="z_l" # Name of depth variable in NetCDF ocean products
+export zmin="0." # Minimum depth to extract from NetCDF ocean products
+export zmax="300." # Maximum depth to extract from NetCDF ocean products
+
+export FHOUT_WAV_EXTRACT=6 # Frequency of wave output to be saved on disk
+
+#Paramater Tables used
+export varlist_2d="${PARMgfs}/product/gefs_shortparmlist_2d.parm" # Parameter table for surface variables
+export varlist_3d="${PARMgfs}/product/gefs_shortparmlist_3d_h.parm" # Parameter table for upper air instantaneous variables
+export varlist_3d_d="${PARMgfs}/product/gefs_shortparmlist_3d_d.parm" # Parameter table for upper air daily-averaged variables
+export varlist_wav="${PARMgfs}/product/gefs_wav_shortparmlist.parm" # Parameter table for wave variables
+export varlist_ocn_netcdf="${PARMgfs}/product/gefs_ocn_shortparmlist.parm" # Parameter table for ocean netcdf variables
+export varlist_ice_netcdf="${PARMgfs}/product/gefs_ice_shortparmlist.parm" # Parameter table for ice netcdf variables
+
+#Directory to save extracted variables
+export ARC_RFCST_PROD="${ARCDIR}/rfcst/${PDY:0:4}/${PDY:0:6}/${PDY:0:8}/mem${ENSMEM}"
+export ARC_RFCST_PROD_ATMOS_F2D="${ARC_RFCST_PROD}/atmos/f2d"
+export ARC_RFCST_PROD_ATMOS_F3D="${ARC_RFCST_PROD}/atmos/f3d"
+export ARC_RFCST_PROD_OCN="${ARC_RFCST_PROD}/ocn"
+export ARC_RFCST_PROD_ICE="${ARC_RFCST_PROD}/ice"
+export ARC_RFCST_PROD_WAV="${ARC_RFCST_PROD}/wav"
+
+echo "END: config.extractvars"
diff --git a/parm/config/gefs/config.resources b/parm/config/gefs/config.resources
index 7c3d77de1d..3e4f05b4c1 100644
--- a/parm/config/gefs/config.resources
+++ b/parm/config/gefs/config.resources
@@ -293,6 +293,18 @@ case ${step} in
     export NTASKS=${npe_wavepostpnt}
     ;;
 
+  "extractvars")
+    export wtime_extractvars="00:30:00"
+    export npe_extractvars=1
+    export nth_extractvars=1
+    export npe_node_extractvars="${npe_extractvars}"
+    export wtime_extractvars_gfs="${wtime_extractvars}"
+    export npe_extractvars_gfs="${npe_extractvars}"
+    export nth_extractvars_gfs="${nth_extractvars}"
+    export npe_node_extractvars_gfs="${npe_node_extractvars}"
+    export is_exclusive=False
+    ;;
+
   *)
     echo "FATAL ERROR: Invalid job ${step} passed to ${BASH_SOURCE[0]}"
     exit 1
diff --git a/parm/config/gefs/yaml/defaults.yaml b/parm/config/gefs/yaml/defaults.yaml
index f6797758da..d2b486e7ca 100644
--- a/parm/config/gefs/yaml/defaults.yaml
+++ b/parm/config/gefs/yaml/defaults.yaml
@@ -7,7 +7,8 @@ base:
   DO_BUFRSND: "NO"
   DO_GEMPAK: "NO"
   DO_AWIPS: "NO"
-  KEEPDATA: "NO" 
+  KEEPDATA: "NO"
+  DO_EXTRACTVARS: "NO"
   FHMAX_GFS: 120
   FHMAX_HF_GFS: 0
   REPLAY_ICS: "NO"
diff --git a/parm/product/gefs_ice_shortparmlist.parm b/parm/product/gefs_ice_shortparmlist.parm
new file mode 100644
index 0000000000..07db948fe3
--- /dev/null
+++ b/parm/product/gefs_ice_shortparmlist.parm
@@ -0,0 +1,10 @@
+aice_h
+hi_h
+Tsfc_h
+uvel_h
+vvel_h
+hs_h
+albsni_h
+melts_h
+meltb_h
+frzmlt_h
diff --git a/parm/product/gefs_ocn_shortparmlist.parm b/parm/product/gefs_ocn_shortparmlist.parm
new file mode 100644
index 0000000000..6673ddb16e
--- /dev/null
+++ b/parm/product/gefs_ocn_shortparmlist.parm
@@ -0,0 +1,9 @@
+temp
+SST
+SSH
+SSS
+MLD_003
+taux
+tauy
+SSU
+SSV
diff --git a/parm/product/gefs_shortparmlist_2d.parm b/parm/product/gefs_shortparmlist_2d.parm
new file mode 100644
index 0000000000..bc13101926
--- /dev/null
+++ b/parm/product/gefs_shortparmlist_2d.parm
@@ -0,0 +1,38 @@
+:PRES:surface:
+:WEASD:surface:
+:TMP:2 m above ground:
+:TMP:surface:
+:RH:2 m above ground:
+:TMAX:2 m above ground:
+:TMIN:2 m above ground:
+:UGRD:10 m above ground:
+:VGRD:10 m above ground:
+:APCP:surface:
+:CSNOW:surface:
+:CICEP:surface:
+:CFRZR:surface:
+:CRAIN:surface:
+:PWAT:entire atmosphere (considered as a single layer):
+:TCDC:entire atmosphere (considered as a single layer):
+:DSWRF:surface:
+:DLWRF:surface:
+:ULWRF:top of atmosphere:
+:HLCY:3000-0 m above ground:
+:CAPE:180-0 mb above ground:
+:CIN:180-0 mb above ground:
+:PRMSL:mean sea level:
+:USWRF:surface:
+:ULWRF:surface:
+:TSOIL:0-0.1 m below ground:
+:TSOIL:0.1-0.4 m below ground:
+:SOILW:0-0.1 m below ground:
+:SOILW:0.1-0.4 m below ground:
+:SOILW:0.4-1 m below ground:
+:SOILW:1-2 m below ground:
+:PEVPR:surface:
+:LHTFL:surface:
+:SHTFL:surface:
+:WATR:surface:
+:TSNOWP:surface:
+:FDNSSTMP:surface:
+:HGT:highest tropospheric freezing level:
diff --git a/parm/product/gefs_shortparmlist_3d_d.parm b/parm/product/gefs_shortparmlist_3d_d.parm
new file mode 100644
index 0000000000..37a2678826
--- /dev/null
+++ b/parm/product/gefs_shortparmlist_3d_d.parm
@@ -0,0 +1,34 @@
+:UGRD:1 mb:
+:UGRD:2 mb:
+:UGRD:3 mb:
+:UGRD:5 mb:
+:UGRD:7 mb:
+:UGRD:10 mb:
+:UGRD:20 mb:
+:UGRD:30 mb:
+:UGRD:50 mb:
+:UGRD:70 mb:
+:VGRD:1 mb:
+:VGRD:2 mb:
+:VGRD:3 mb:
+:VGRD:5 mb:
+:VGRD:7 mb:
+:VGRD:10 mb:
+:VGRD:20 mb:
+:VGRD:30 mb:
+:VGRD:50 mb:
+:VGRD:70 mb:
+:TMP:1 mb:
+:TMP:2 mb:
+:TMP:3 mb:
+:TMP:5 mb:
+:TMP:7 mb:
+:TMP:10 mb:
+:TMP:20 mb:
+:TMP:30 mb:
+:TMP:50 mb:
+:TMP:70 mb:
+:HGT:10 mb:
+:HGT:50 mb:
+:O3MR:10 mb:
+:O3MR:50 mb:
diff --git a/parm/product/gefs_shortparmlist_3d_h.parm b/parm/product/gefs_shortparmlist_3d_h.parm
new file mode 100644
index 0000000000..d7241f633c
--- /dev/null
+++ b/parm/product/gefs_shortparmlist_3d_h.parm
@@ -0,0 +1,45 @@
+:HGT:100 mb:
+:TMP:100 mb:
+:UGRD:100 mb:
+:VGRD:100 mb:
+:O3MR:100 mb:
+:HGT:200 mb:
+:TMP:200 mb:
+:RH:200 mb:
+:UGRD:200 mb:
+:VGRD:200 mb:
+:HGT:250 mb:
+:TMP:250 mb:
+:RH:250 mb:
+:UGRD:250 mb:
+:VGRD:250 mb:
+:HGT:500 mb:
+:TMP:500 mb:
+:RH:500 mb:
+:UGRD:500 mb:
+:VGRD:500 mb:
+:HGT:700 mb:
+:TMP:700 mb:
+:RH:700 mb:
+:UGRD:700 mb:
+:VGRD:700 mb:
+:HGT:850 mb:
+:TMP:850 mb:
+:RH:850 mb:
+:VVEL:850 mb:
+:UGRD:850 mb:
+:VGRD:850 mb:
+:HGT:925 mb:
+:TMP:925 mb:
+:RH:925 mb:
+:UGRD:925 mb:
+:VGRD:925 mb:
+:TMP:1000 mb:
+:RH:1000 mb:
+:UGRD:1000 mb:
+:VGRD:1000 mb:
+:HGT:1000 mb:
+:TMP:0.995 sigma level:
+:RH:0.995 sigma level:
+:UGRD:0.995 sigma level:
+:VGRD:0.995 sigma level:
diff --git a/parm/product/gefs_wav_shortparmlist.parm b/parm/product/gefs_wav_shortparmlist.parm
new file mode 100644
index 0000000000..a45e023c85
--- /dev/null
+++ b/parm/product/gefs_wav_shortparmlist.parm
@@ -0,0 +1,3 @@
+:UGRD:surface:
+:VGRD:surface:
+:HTSGW:surface:
diff --git a/scripts/exglobal_extractvars.sh b/scripts/exglobal_extractvars.sh
new file mode 100755
index 0000000000..a124667679
--- /dev/null
+++ b/scripts/exglobal_extractvars.sh
@@ -0,0 +1,53 @@
+#! /usr/bin/env bash                                                                                                                                                                    
+
+################################################################################
+## UNIX Script Documentation Block
+## Script name:         exglobal_extractvars.sh
+## Script description:  Extracts variables from atmosphere, ocean, ice and wave
+##                      products and saves these variables in arcdir
+#######################
+# Main body starts here
+#######################
+
+source "${USHgfs}/preamble.sh"
+source "${USHgfs}/extractvars_tools.sh"
+
+# Scripts used
+EXTRCTVARA="${USHgfs}/atmos_extractvars.sh"
+EXTRCTVARO="${USHgfs}/ocnice_extractvars.sh"
+EXTRCTVARW="${USHgfs}/wave_extractvars.sh"
+
+# Set FHMAX_HF_GFS equal to FHMAX_GFS if FHMAX_HF_GFS is greater than FHMAX_GFS
+if (( FHMAX_GFS < FHMAX_HF_GFS )); then
+  export FHMAX_HF_GFS=${FHMAX_GFS}
+fi
+
+# Set FHOUT_WAV_EXTRACT equal to FHOUT_WAV if FHOUT_WAV is not a factor of FHOUT_WAV_EXTRACT
+if (( FHOUT_WAV_EXTRACT % FHOUT_WAV != 0 )); then
+  FHOUT_WAV_EXTRACT=${FHOUT_WAV}
+fi
+
+# Extract variables for atmosphere
+if [[ "${DO_ATM}" == "YES" ]]; then
+  ${EXTRCTVARA} "${DATA}/atmos"
+fi
+
+# Extract variables for ocean
+if [[ "${DO_OCN}" == "YES" ]]; then
+  export component_name="ocn"
+  ${EXTRCTVARO} "${DATA}/ocn" "${varlist_ocn_netcdf}" "${ocnres}" "${compress_ocn}" "${FHOUT_OCN_GFS}" "${ARC_RFCST_PROD_OCN}"
+fi
+
+# Extract variables for ice
+if [[ "${DO_ICE}" == "YES" ]]; then
+  export component_name="ice"
+  ${EXTRCTVARO} "${DATA}/ice" "${varlist_ice_netcdf}" "${iceres}" "${compress_ice}" "${FHOUT_ICE_GFS}" "${ARC_RFCST_PROD_ICE}"
+fi
+
+# Extract variables for wave
+if [[ "${DO_WAVE}" == "YES" ]]; then
+  export component_name="wav"
+  ${EXTRCTVARW} "${DATA}/wav"
+fi
+
+exit 0
diff --git a/ush/atmos_extractvars.sh b/ush/atmos_extractvars.sh
new file mode 100755
index 0000000000..70e86b2f4e
--- /dev/null
+++ b/ush/atmos_extractvars.sh
@@ -0,0 +1,98 @@
+#! /usr/bin/env bash
+
+################################################################################
+## UNIX Script Documentation Block
+## Script name:         atmos_extractvars.sh
+## Script description:  Extracts and calculates 24-hr averages of variables
+##                      from atmosphere products and saves these variables in arcdir
+#######################
+# Main body starts here
+#######################
+
+source "${USHgfs}/preamble.sh"
+
+fcnt=1 # 1 is 1st quarter, 2 is 2nd quarter and 3 is 3rd quarter of the day
+dcnt=1 # lead day
+subdata=${1}
+
+[[ -d "${subdata}" ]] || mkdir -p "${subdata}"
+
+for outtype in "f2d" "f3d"; do
+
+  if [[ "${outtype}" == "f2d" ]]; then
+    varlist=${varlist_2d}
+    ARC_RFCST_PROD_ATMOS="${ARC_RFCST_PROD_ATMOS_F2D}"
+  elif [[ "${outtype}" == "f3d" ]]; then
+    varlist=${varlist_3d}
+    varlist_d=${varlist_3d_d}
+    ARC_RFCST_PROD_ATMOS="${ARC_RFCST_PROD_ATMOS_F3D}"
+  fi
+
+  outdirpre="${subdata}/${outtype}"
+  [[ -d "${outdirpre}" ]] || mkdir -p "${outdirpre}"
+
+  nh=${FHMIN}
+  while (( nh <= FHMAX_GFS )); do
+    fnh=$(printf "%3.3d" "${nh}")
+
+    if [[ "${outtype}" == "f2d" ]]; then
+      if (( nh < FHMAX_HF_GFS )); then
+        outres="0p25"
+      else
+        outres="0p50"
+      fi
+    elif [[ "${outtype}" == "f3d" ]]; then
+      outres="1p00"
+    fi
+
+    if (( nh <= FHMAX_HF_GFS )); then
+      outfreq=${FHOUT_HF_GFS}
+    else
+      outfreq=${FHOUT_GFS}
+    fi                                      
+
+    com_var="COMIN_ATMOS_GRIB_${outres}"
+    infile1="${!com_var}/${RUN}.t${cyc}z.pgrb2.${outres}.f${fnh}"
+    infile2="${!com_var}/${RUN}.t${cyc}z.pgrb2b.${outres}.f${fnh}"
+    outfile="${outdirpre}/${RUN}.t${cyc}z.pgrb2.${outres}.f${fnh}"
+    rm -f "${outfile}" #remove outfile if it already exists before extraction
+
+    for infile in "${infile1}" "${infile2}"; do
+      if [[ -f "${infile}" ]]; then # check if input file exists before extraction
+        # shellcheck disable=SC2312
+        ${WGRIB2} "${infile}" | grep -F -f "${varlist}" | ${WGRIB2} -i "${infile}" -append -grib "${outfile}"
+      else
+        echo "WARNING: ${infile} does not exist."
+      fi
+    done
+
+    check_atmos "${infile1}" "${infile2}" "${varlist}" "${fnh}"
+    copy_to_comout "${outfile}" "${ARC_RFCST_PROD_ATMOS}"
+
+    # Compute daily average for a subset of variables
+    if (( nh % 6 == 0 )) && (( nh != 0 )) && [[ "${outtype}" == "f3d" ]]; then
+      outfile=${subdata}/vartmp_raw_vari_ldy${dcnt}.grib2
+      for infile in "${infile1}" "${infile2}"; do
+        if [[ -f "${infile}" ]]; then # check if input file exists before extraction
+          # shellcheck disable=SC2312
+          ${WGRIB2} "${infile}" | grep -F -f "${varlist_d}" | ${WGRIB2} -i "${infile}" -append -grib "${outfile}"
+        else
+          echo "WARNING: ${infile} does not exist."
+        fi
+      done
+      if [[ ${fcnt} -eq 4 ]]; then
+        daily_avg_atmos "${outfile}" "${dcnt}" "${outres}"
+        copy_to_comout "${davg_file}" "${ARC_RFCST_PROD_ATMOS}"
+        fcnt=1
+        dcnt=$(( dcnt + 1 ))
+      else
+        fcnt=$(( fcnt + 1 ))
+      fi # If at final lead hour of a given day
+    fi # if lead hour is divisible by 6 and outtype is f3d
+
+    nh=$(( nh + outfreq ))
+  done # nh
+
+done # f2d,f3d
+
+exit 0
diff --git a/ush/extractvars_tools.sh b/ush/extractvars_tools.sh
new file mode 100644
index 0000000000..daf61a3d2e
--- /dev/null
+++ b/ush/extractvars_tools.sh
@@ -0,0 +1,60 @@
+#! /usr/bin/env bash
+
+check_atmos() {
+  # Function to check if there are any missing parm variables in any of the input product grib2 files
+  # A warning will be displayed if there is a parm variable that cannot be found in any of the given input product grib2 files
+  infile1p=$1
+  infile2p=$2
+  varlistl=$3
+  fnhl=$4
+  requestedvar_in_allgrb2file="${subdata}/parmvarsingribfil.txt"
+  rm -f "${requestedvar_in_allgrb2file}"
+  touch "${requestedvar_in_allgrb2file}"
+  for infilep in "${infile1p}" "${infile2p}"; do
+    # It is permitted for an empty string to return if no parmlist vars are in infilep, therefore do not return exit 1 error
+    # shellcheck disable=SC2312
+    ${WGRIB2} "${infilep}" | grep -F -f "${varlist}" >> "${requestedvar_in_allgrb2file}" || true
+  done
+  mapfile -t requestedvar_in_allgrb2file_arr < "${requestedvar_in_allgrb2file}"
+  while read -r vari; do
+    if [[ ! ${requestedvar_in_allgrb2file_arr[*]} =~ ${vari} ]] ;then
+      echo "WARNING: PARM VARIABLE (${vari}) is not available in pgrb and pgrb2b for f${fnhl}."
+    fi
+  done <"${varlistl}"
+}
+
+daily_avg_atmos() {
+  # Function to calculate the 24-hr average of a grib2 file with atmospheric fields
+  # The input grib2 file must contain all the time records to be averaged (e.g. 6hr, 12hr, 18hr and 24hr record in one grib2 file)
+  outfile_p=$1
+  dcnt_p=$2
+  outres_p=$3
+  fnd=$(printf "%2.2d" "${dcnt_p}")
+  davg_file=${outdirpre}/${RUN}.t${cyc}z.pgrb2.${outres_p}.24hr_avg.ldy${fnd}
+  vcnt=1 #count variables in varlist_d
+  while read -r vari; do
+    davgtmp=${subdata}/atmos_tmp.ldy${fnd}.${vcnt}
+    # shellcheck disable=SC2312
+    ${WGRIB2} "${outfile_p}" | grep "${vari}" | ${WGRIB2} -i "${outfile_p}" -fcst_ave 6hr "${davgtmp}"
+    # shellcheck disable=SC2312
+    ${WGRIB2} "${davgtmp}" | ${WGRIB2} -i "${davgtmp}" -append -grib "${davg_file}"
+    rm -f "${davgtmp}"
+    vcnt=$(( vcnt + 1 ))
+  done <"${varlist_d}" # variable
+}
+
+copy_to_comout() {
+  # Function to copy the output file with the extracted product variables to a user-defined destination directory
+  rundir_outfile=$1 # output data file generated in RUNDIR
+  comout_dir=$2 # destination directory to which to copy the data file
+  if [[ -f "${rundir_outfile}" ]]; then
+    cpfs "${rundir_outfile}" "${comout_dir}"
+  else
+    echo "FATAL ERROR: Output file (${rundir_outfile}) does not exist."
+    export err=1; err_chk
+  fi
+}
+
+declare -xf check_atmos
+declare -xf daily_avg_atmos
+declare -xf copy_to_comout
diff --git a/ush/ocnice_extractvars.sh b/ush/ocnice_extractvars.sh
new file mode 100755
index 0000000000..f0660bb6ec
--- /dev/null
+++ b/ush/ocnice_extractvars.sh
@@ -0,0 +1,66 @@
+#! /usr/bin/env bash
+
+################################################################################
+## UNIX Script Documentation Block
+## Script name:         ocnice_extractvars.sh
+## Script description:  Extracts and optionally compresses variables
+##                      from ocean and ice products
+##                      and saves these variables in arcdir
+#######################
+# Main body starts here
+#######################
+
+source "${USHgfs}/preamble.sh"
+
+subdata=${1}
+varlist=${2}
+datares=${3}
+datacompress=${4}
+fhout_ocnice=${5}
+comout_rfcst_prod_ocnice=${6}
+
+[[ -d "${subdata}" ]] || mkdir -p "${subdata}"
+
+for (( nh = FHMIN_GFS; nh <= FHMAX_GFS; nh = nh + fhout_ocnice )); do
+  fnh=$(printf "%3.3d" "${nh}")
+
+  if [[ ${component_name} == "ocn" ]]; then
+    infile=${COMIN_OCEAN_NETCDF}/${RUN}.ocean.t${cyc}z.${datares}.f${fnh}.nc
+    # For ocean products, add an argument to extract a subset of levels
+    otherargs=(-d "${depthvar_name},""${zmin},""${zmax}")
+  elif [[ ${component_name} == "ice" ]]; then
+    infile=${COMIN_ICE_NETCDF}/${RUN}.ice.t${cyc}z.${datares}.f${fnh}.nc
+    otherargs=()
+  fi
+  outfile=${subdata}/${RUN}.${component_name}.t${cyc}z.${datares}.f${fnh}.nc
+
+  if [[ -f "${infile}" ]]; then #check if input file exists before extraction
+    varsrequested=$(paste -s "${varlist}")
+    varsinfile=$(cdo -showname "${infile}")
+    varsavailable=""
+    for i in ${varsrequested}; do
+      # Check if variable from parm file is available in netcdf file. If variable is not in netcdf file, do not try to extract that variable.
+      if [[ ${varsinfile} == *"${i}"* ]]; then
+        varsavailable+="${i},"
+      else
+        echo "WARNING: ${i} is not available in ${infile}."
+      fi
+    done
+    if [[ -z "${varsavailable}" ]]; then
+      echo "WARNING: No variables from parm file ${varlist} are available in netcdf file ${infile}."
+    else
+      ocnice_vars=${varsavailable::-1}
+      ncks -v "${ocnice_vars}" "${otherargs[@]}" "${infile}" "${outfile}"
+    fi
+    if [[ ${datacompress} -eq 1 ]]; then
+      ${COMPRSCMD} "${outfile}"
+      copy_to_comout "${outfile}.bz2" "${comout_rfcst_prod_ocnice}"
+    else
+      copy_to_comout "${outfile}" "${comout_rfcst_prod_ocnice}"
+    fi 
+  else
+    echo "WARNING: ${infile} does not exist."
+  fi
+done # nh
+
+exit 0                                                                                                                                                                                        
diff --git a/ush/wave_extractvars.sh b/ush/wave_extractvars.sh
new file mode 100755
index 0000000000..32ee44986b
--- /dev/null
+++ b/ush/wave_extractvars.sh
@@ -0,0 +1,34 @@
+#! /usr/bin/env bash                                                                                                                                                                          
+
+################################################################################
+## UNIX Script Documentation Block
+## Script name:         wave_extractvars.sh
+## Script description:  Extracts variables from wave products
+##                      and saves these variables in arcdir
+#######################
+# Main body starts here
+#######################
+
+source "${USHgfs}/preamble.sh"
+
+subdata=${1}
+
+[[ -d "${subdata}" ]] || mkdir -p "${subdata}"
+
+for (( nh = FHOUT_WAV_EXTRACT; nh <= FHMAX_WAV; nh = nh + FHOUT_WAV_EXTRACT )); do
+  fnh=$(printf "%3.3d" "${nh}")
+
+  infile=${COMIN_WAVE_GRID}/${RUN}wave.t${cyc}z.global.${wavres}.f${fnh}.grib2
+  outfile=${subdata}/${RUN}wave.t${cyc}z.global.${wavres}.f${fnh}.grib2
+  rm -f "${outfile}" # Remove outfile if it already exists before extraction
+
+  if [[ -f "${infile}" ]]; then # Check if input file exists before extraction
+    # shellcheck disable=SC2312 
+    ${WGRIB2} "${infile}" | grep -F -f "${varlist_wav}" | ${WGRIB2} -i "${infile}" -append -grib "${outfile}"
+  else
+    echo "WARNING: ${infile} does not exist."
+  fi 
+  copy_to_comout "${outfile}" "${ARC_RFCST_PROD_WAV}"
+done # nh
+
+exit 0                                                                                                                                                                                        
diff --git a/workflow/applications/applications.py b/workflow/applications/applications.py
index 64440830c1..95d7cd9bb3 100644
--- a/workflow/applications/applications.py
+++ b/workflow/applications/applications.py
@@ -65,6 +65,7 @@ def __init__(self, conf: Configuration) -> None:
         self.do_upp = not _base.get('WRITE_DOPOST', True)
         self.do_goes = _base.get('DO_GOES', False)
         self.do_mos = _base.get('DO_MOS', False)
+        self.do_extractvars = _base.get('DO_EXTRACTVARS', False)
 
         self.do_hpssarch = _base.get('HPSSARCH', False)
 
diff --git a/workflow/applications/gefs.py b/workflow/applications/gefs.py
index c165f9d1ca..bdff2186d0 100644
--- a/workflow/applications/gefs.py
+++ b/workflow/applications/gefs.py
@@ -30,6 +30,9 @@ def _get_app_configs(self):
         if self.do_aero:
             configs += ['prep_emissions']
 
+        if self.do_extractvars:
+            configs += ['extractvars']
+
         return configs
 
     @staticmethod
@@ -73,4 +76,7 @@ def get_task_names(self):
                 tasks += ['wavepostbndpnt', 'wavepostbndpntbll']
             tasks += ['wavepostpnt']
 
+        if self.do_extractvars:
+            tasks += ['extractvars']
+
         return {f"{self._base['CDUMP']}": tasks}
diff --git a/workflow/rocoto/gefs_tasks.py b/workflow/rocoto/gefs_tasks.py
index cfd8fa7093..f6d1ad6cdc 100644
--- a/workflow/rocoto/gefs_tasks.py
+++ b/workflow/rocoto/gefs_tasks.py
@@ -483,3 +483,48 @@ def wavepostpnt(self):
         task = rocoto.create_task(member_metatask_dict)
 
         return task
+
+    def extractvars(self):
+        deps = []
+        if self.app_config.do_wave:
+            dep_dict = {'type': 'task', 'name': 'wave_post_grid_mem#member#'}
+            deps.append(rocoto.add_dependency(dep_dict))
+        if self.app_config.do_ocean:
+            dep_dict = {'type': 'metatask', 'name': 'ocean_prod_#member#'}
+            deps.append(rocoto.add_dependency(dep_dict))
+        if self.app_config.do_ice:
+            dep_dict = {'type': 'metatask', 'name': 'ice_prod_#member#'}
+            deps.append(rocoto.add_dependency(dep_dict))
+        if self.app_config.do_atm:
+            dep_dict = {'type': 'metatask', 'name': 'atmos_prod_#member#'}
+            deps.append(rocoto.add_dependency(dep_dict))
+        dependencies = rocoto.create_dependency(dep_condition='and', dep=deps)
+        extractvars_envars = self.envars.copy()
+        extractvars_dict = {'ENSMEM': '#member#',
+                            'MEMDIR': 'mem#member#',
+                            }
+        for key, value in extractvars_dict.items():
+            extractvars_envars.append(rocoto.create_envar(name=key, value=str(value)))
+
+        resources = self.get_resource('extractvars')
+        task_name = f'extractvars_mem#member#'
+        task_dict = {'task_name': task_name,
+                     'resources': resources,
+                     'dependency': dependencies,
+                     'envars': extractvars_envars,
+                     'cycledef': 'gefs',
+                     'command': f'{self.HOMEgfs}/jobs/rocoto/extractvars.sh',
+                     'job_name': f'{self.pslot}_{task_name}_@H',
+                     'log': f'{self.rotdir}/logs/@Y@m@d@H/{task_name}.log',
+                     'maxtries': '&MAXTRIES;'
+                     }
+
+        member_var_dict = {'member': ' '.join([str(mem).zfill(3) for mem in range(0, self.nmem + 1)])}
+        member_metatask_dict = {'task_name': 'extractvars',
+                                'task_dict': task_dict,
+                                'var_dict': member_var_dict
+                                }
+
+        task = rocoto.create_task(member_metatask_dict)
+
+        return task