diff --git a/parm/soca/berror/saber_blocks.yaml b/parm/soca/berror/saber_blocks.yaml
index 9487ddae3..b8965cf05 100644
--- a/parm/soca/berror/saber_blocks.yaml
+++ b/parm/soca/berror/saber_blocks.yaml
@@ -60,7 +60,7 @@ components:
         ice_filename: 'ice.%mem%.nc'
         state variables: [tocn, socn, ssh, uocn, vocn, cicen, hicen, hsnon]
       pattern: '%mem%'
-      nmembers: ${CLIM_ENS_SIZE}
+      nmembers: ${ENS_SIZE}
     localization:
       localization method: SABER
       saber central block:
diff --git a/parm/soca/berror/soca_ensb.yaml b/parm/soca/berror/soca_ensb.yaml
index 447d5b5fd..42e054531 100644
--- a/parm/soca/berror/soca_ensb.yaml
+++ b/parm/soca/berror/soca_ensb.yaml
@@ -17,7 +17,7 @@ vertical geometry:
   read_from_file: 3
 
 soca increments:
-  number of increments: ${CLIM_ENS_SIZE}
+  number of increments: ${ENS_SIZE}
   pattern: '%mem%'
   template:
     date: '{{ATM_WINDOW_BEGIN}}'
diff --git a/scripts/exgdas_global_marine_analysis_bmat_vrfy.sh b/scripts/exgdas_global_marine_analysis_bmat_vrfy.sh
deleted file mode 100755
index 4903d7d40..000000000
--- a/scripts/exgdas_global_marine_analysis_bmat_vrfy.sh
+++ /dev/null
@@ -1,80 +0,0 @@
-#!/bin/bash
-################################################################################
-####  UNIX Script Documentation Block
-#                      .                                             .
-# Script name:         exgdas_global_marine_analysis_bmat_vrfy.sh
-# Script description:  Diagnose the JEDI/SOCA B-matrix
-#
-# Author: Guillaume Vernieres        Org: NCEP/EMC     Date: 2023-01-12
-#
-# Abstract: This script runs the JEDI/SOCA dirac application with the B-matrix
-#           used in the variational step
-# $Id$
-#
-# Attributes:
-#   Language: POSIX shell
-#   Machine: Orion
-#
-################################################################################
-
-#  Set environment.
-export VERBOSE=${VERBOSE:-"YES"}
-if [ $VERBOSE = "YES" ]; then
-   echo $(date) EXECUTING $0 $* >&2
-   set -x
-fi
-
-#  Directories
-pwd=$(pwd)
-
-#  Utilities
-export NLN=${NLN:-"/bin/ln -sf"}
-export DIRAC="${HOMEgfs}/sorc/gdas.cd/ush/ufsda/dirac_yaml.py"
-
-################################################################################
-# Compute the impulse response of the B-matrix from the variational application
-# field = 1: tocn
-# field = 2: socn
-# field = 3: ssh
-# field = 4: cicen
-# field = 5: hicen
-
-arr=("tocn" "socn" "ssh" "cicen" "hicen")
-level=1
-ndiracs=100
-
-for i in "${!arr[@]}"
-do
-    var=${arr[i]}
-    ifield=$((i+1))
-    # generate dirac yamls
-cat > dirac_output.yaml << EOL
-datadir: ./Data
-exp: dirac_${var}_${level}
-type: an
-EOL
-
-    ${DIRAC} --varyaml 'var.yaml' \
-             --fields 'soca_gridspec.nc' \
-             --dim1 'xaxis_1' \
-             --dim2 'yaxis_1' \
-             --diracyaml 'dirac.yaml' \
-             --ndiracs ${ndiracs} \
-             --level ${level} \
-             --fieldindex ${ifield} \
-             --diracoutput dirac_output.yaml
-    export err=$?
-
-    # run the dirac application
-    $APRUN_OCNANAL $JEDI_BIN/soca_dirac.x dirac.yaml
-    export err=$?
-    if [ $err -gt 0  ]; then
-        exit $err
-    fi
-done
-
-################################################################################
-
-exit ${err}
-
-################################################################################
diff --git a/scripts/exgdas_global_marine_analysis_prep.py b/scripts/exgdas_global_marine_analysis_prep.py
index 4e0635be2..aaeadc612 100755
--- a/scripts/exgdas_global_marine_analysis_prep.py
+++ b/scripts/exgdas_global_marine_analysis_prep.py
@@ -89,6 +89,10 @@ def cice_hist2fms(input_filename, output_filename):
     # open the CICE history file
     ds = xr.open_dataset(input_filename)
 
+    if 'aicen' in ds.variables and 'hicen' in ds.variables and 'hsnon' in ds.variables:
+        logging.info(f"*** Already reformatted, skipping.")
+        return
+
     # rename the dimensions to xaxis_1 and yaxis_1
     ds = ds.rename({'ni': 'xaxis_1', 'nj': 'yaxis_1'})
 
@@ -224,6 +228,11 @@ def find_clim_ens(input_date):
 comin_obs = os.getenv('COMIN_OBS')  # R2D2 DB for now
 anl_dir = os.getenv('DATA')
 staticsoca_dir = os.getenv('SOCA_INPUT_FIX_DIR')
+if os.getenv('DOHYBVAR') == "YES":
+    dohybvar = True
+    nmem_ens = int(os.getenv('NMEM_ENS'))
+else:
+    dohybvar = False
 
 # create analysis directories
 diags = os.path.join(anl_dir, 'diags')            # output dir for soca DA obs space
@@ -242,6 +251,7 @@ def find_clim_ens(input_date):
 fcst_begin = datetime.strptime(os.getenv('PDY')+os.getenv('cyc'), '%Y%m%d%H')
 RUN = os.getenv('RUN')
 cyc = os.getenv('cyc')
+gcyc = os.getenv('gcyc')
 PDY = os.getenv('PDY')
 
 ################################################################################
@@ -310,21 +320,43 @@ def find_clim_ens(input_date):
 FileHandler({'copy': bkgerr_list}).sync()
 
 ################################################################################
-# stage static ensemble
-
-logging.info(f"---------------- Stage climatological ensemble")
-clim_ens_member_list = []
-clim_ens_dir = find_clim_ens(pytz.utc.localize(window_begin, is_dst=None))
-clim_ens_size = len(glob.glob(os.path.abspath(os.path.join(clim_ens_dir, 'ocn.*.nc'))))
-os.environ['CLIM_ENS_SIZE'] = str(clim_ens_size)
-
-for domain in ['ocn', 'ice']:
-    for mem in range(1, clim_ens_size+1):
-        fname = domain+"."+str(mem)+".nc"
-        fname_in = os.path.abspath(os.path.join(clim_ens_dir, fname))
-        fname_out = os.path.abspath(os.path.join(static_ens, fname))
-        clim_ens_member_list.append([fname_in, fname_out])
-FileHandler({'copy': clim_ens_member_list}).sync()
+# stage ensemble members
+if dohybvar:
+    logging.info("---------------- Stage ensemble members")
+    nmem_ens = int(os.getenv('NMEM_ENS'))
+    longname = {'ocn': 'ocean', 'ice': 'ice'}
+    ens_member_list = []
+    for mem in range(1, nmem_ens+1):
+        for domain in ['ocn', 'ice']:
+            # TODO(Guillaume): make use and define ensemble COM in the j-job
+            ensdir = os.path.join(os.getenv('COM_OCEAN_HISTORY_PREV'), '..', '..', '..', '..', '..',
+                                  f'enkf{RUN}.{PDY}', f'{gcyc}', f'mem{str(mem).zfill(3)}',
+                                  'model_data', longname[domain], 'history')
+            ensdir = os.path.normpath(ensdir)
+            f009 = f'enkfgdas.t{gcyc}z.{domain}f009.nc'
+
+            fname_in = os.path.abspath(os.path.join(ensdir, f009))
+            fname_out = os.path.abspath(os.path.join(static_ens, domain+"."+str(mem)+".nc"))
+            ens_member_list.append([fname_in, fname_out])
+    FileHandler({'copy': ens_member_list}).sync()
+
+    # reformat the cice history output
+    for mem in range(1, nmem_ens+1):
+        cice_fname = os.path.abspath(os.path.join(static_ens, "ice."+str(mem)+".nc"))
+        cice_hist2fms(cice_fname, cice_fname)
+else:
+    logging.info("---------------- Stage offline ensemble members")
+    ens_member_list = []
+    clim_ens_dir = find_clim_ens(pytz.utc.localize(window_begin, is_dst=None))
+    nmem_ens = len(glob.glob(os.path.abspath(os.path.join(clim_ens_dir, 'ocn.*.nc'))))
+    for domain in ['ocn', 'ice']:
+        for mem in range(1, nmem_ens+1):
+            fname = domain+"."+str(mem)+".nc"
+            fname_in = os.path.abspath(os.path.join(clim_ens_dir, fname))
+            fname_out = os.path.abspath(os.path.join(static_ens, fname))
+            ens_member_list.append([fname_in, fname_out])
+    FileHandler({'copy': ens_member_list}).sync()
+os.environ['ENS_SIZE'] = str(nmem_ens)
 
 ################################################################################
 # prepare JEDI yamls
@@ -456,7 +488,6 @@ def find_clim_ens(input_date):
 else:
     logging.info(f"using default SABER blocks yaml")
     os.environ['SABER_BLOCKS_YAML'] = os.path.join(gdas_home, 'parm', 'soca', 'berror', 'saber_blocks.yaml')
-os.environ['CLIM_ENS_SIZE'] = str(clim_ens_size)
 
 # substitute templated variables in the var config
 logging.info(f"{config}")
diff --git a/test/soca/gw/CMakeLists.txt b/test/soca/gw/CMakeLists.txt
index 575bc7f2b..c5a3a26ca 100644
--- a/test/soca/gw/CMakeLists.txt
+++ b/test/soca/gw/CMakeLists.txt
@@ -39,7 +39,6 @@ file(MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/test/soca/gw/apps_scratch)
 # Test JGDAS_GLOBAL_OCEAN_ANALYSIS_*
 set(jjob_list "JGDAS_GLOBAL_OCEAN_ANALYSIS_PREP"
               "JGDAS_GLOBAL_OCEAN_ANALYSIS_BMAT"
-#              "JGDAS_GLOBAL_OCEAN_ANALYSIS_BMAT_VRFY"
               "JGDAS_GLOBAL_OCEAN_ANALYSIS_RUN"
               "JGDAS_GLOBAL_OCEAN_ANALYSIS_CHKPT"
               "JGDAS_GLOBAL_OCEAN_ANALYSIS_POST"
diff --git a/test/soca/gw/prep.sh b/test/soca/gw/prep.sh
index a3590752b..2aeadf224 100755
--- a/test/soca/gw/prep.sh
+++ b/test/soca/gw/prep.sh
@@ -60,3 +60,21 @@ ice_files=$(ls ${clim_ens_dir}/ice.*.nc)
 for ice_file in ${ice_files}; do
     ncrename -O -d ni,xaxis_1 -d nj,yaxis_1 -v aice_h,aicen -v hi_h,hicen -v hs_h,hsnon ${ice_file}
 done
+
+# Invent ensemble forecast for f009
+fake_ocean_members=`ls`
+fake_ice_members=`ls ${clim_ens_dir/ice.*.nc}`
+COMENS=${project_binary_dir}/test/soca/gw/COM/enkfgdas.20180415
+for mem in  {1..4}
+do
+    echo "member: $mem"
+    # ocean member
+    oceandir=${COMENS}/06/mem00${mem}/model_data/ocean/history
+    mkdir -p $oceandir
+    cp ${clim_ens_dir}/ocn.${mem}.nc $oceandir/enkfgdas.t06z.ocnf009.nc
+    echo ${clim_ens_dir}/ocn.${mem}.nc $oceandir/enkfgdas.t06z.ocnf009.nc
+    # ice member
+    icedir=${COMENS}/06/mem00${mem}/model_data/ice/history
+    mkdir -p $icedir
+    cp ${clim_ens_dir}/ice.${mem}.nc $icedir/enkfgdas.t06z.icef009.nc
+done
diff --git a/test/soca/gw/run_jjobs.yaml.test b/test/soca/gw/run_jjobs.yaml.test
index f908c841f..3270d820a 100644
--- a/test/soca/gw/run_jjobs.yaml.test
+++ b/test/soca/gw/run_jjobs.yaml.test
@@ -43,6 +43,8 @@ setup_expt config:
     DO_JEDIOCNVAR: "YES"
     DO_JEDILANDDA: "NO"
     DO_MERGENSST: "NO"
+    NMEM_ENS: "4"
+    DOHYBVAR: "YES"
   ocnanal:
     SOCA_INPUT_FIX_DIR: @HOMEgfs@/sorc/gdas.cd/build/soca_static
     CASE_ANL: C48
diff --git a/ush/soca/run_jjobs.py b/ush/soca/run_jjobs.py
index d4149633c..80d8a3c54 100755
--- a/ush/soca/run_jjobs.py
+++ b/ush/soca/run_jjobs.py
@@ -190,6 +190,14 @@ def copy_bkgs(self):
         self.f.write(f"cp {com_ice_history_src}/*icef*.nc $COM_ICE_HISTORY_PREV \n")
         self.f.write(f"cp {com_ice_restart_src}/*cice_model*.nc $COM_ICE_RESTART_PREV \n")
 
+        # copy ensemble members
+        ensbkgs = os.path.join(self.com_src, f"enkf{self.RUN}.{self.gPDY}")
+        if os.path.exists(os.path.join(ensbkgs, self.gcyc)):
+            self.f.write(f"cp -r {ensbkgs} $ROTDIR \n")
+        else:
+            print('Aborting, ensemble backgrounds not found')
+            sys.exit()
+
     def fixconfigs(self):
         """
         Replace cone of the env. var. in the configs