diff --git a/scripts/exgdas_global_marine_analysis_prep.py b/scripts/exgdas_global_marine_analysis_prep.py
index 6b216c011..f66bdb09e 100755
--- a/scripts/exgdas_global_marine_analysis_prep.py
+++ b/scripts/exgdas_global_marine_analysis_prep.py
@@ -162,9 +162,8 @@ def gen_bkg_list(bkg_path, out_path, window_begin=' ', file_type='gdas.t*.ocnf00
 ufsda.mkdir(bkg_dir)
 
 # create output directory for soca DA
-anl_out = os.path.join(comout, 'ocnanal_' + os.getenv('cyc'), 'Data')
+anl_out = os.path.join(anl_dir, 'Data')
 ufsda.mkdir(anl_out)
-ufsda.symlink(anl_out, os.path.join(anl_dir, 'Data'), remove=False)
 
 ################################################################################
 # fetch observations
diff --git a/scripts/exgdas_global_marine_analysis_run.sh b/scripts/exgdas_global_marine_analysis_run.sh
index e313de089..463e03cc4 100755
--- a/scripts/exgdas_global_marine_analysis_run.sh
+++ b/scripts/exgdas_global_marine_analysis_run.sh
@@ -84,7 +84,7 @@ export err=$?; err_chk
 # Note: ${DATA}/INPUT/MOM.res.nc points to the MOM6 history file from the start of the window
 #       and is used to get the valid vertical geometry of the increment
 socaincr=$(ls -t ${DATA}/Data/ocn.*3dvar*.incr* | head -1)
-mom6incr=${COMOUT}/inc.nc
+mom6incr=${DATA}/inc.nc
 ( socaincr2mom6 ${socaincr} ${DATA}/INPUT/MOM.res.nc ${DATA}/soca_gridspec.nc ${mom6incr} )
 
 
diff --git a/scripts/exgdas_global_marine_analysis_vrfy.py b/scripts/exgdas_global_marine_analysis_vrfy.py
new file mode 100755
index 000000000..1de770cdf
--- /dev/null
+++ b/scripts/exgdas_global_marine_analysis_vrfy.py
@@ -0,0 +1,191 @@
+#!/usr/bin/env python3
+################################################################################
+#  UNIX Script Documentation Block
+#                      .                                             .
+# Script name:         exufsda_global_marine_analysis_vrfy.py
+# Script description:  State and observation space verification for the
+#                      UFS Global Marine Analysis
+#
+# Author: Guillaume Vernieres      Org: NCEP/EMC     Date: 2023-01-23
+#
+# Abstract: This script produces figures relevant to the marine DA cycle
+#
+# $Id$
+#
+# Attributes:
+#   Language: Python3
+#
+################################################################################
+
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+import xarray as xr
+import cartopy
+import cartopy.crs as ccrs
+
+
+def plot_config(grid_file=[], data_file=[],
+                variable=[], levels=[], bounds=[], colormap=[], comout=[], lats=[]):
+    """
+    Prepares the configuration for the plotting functions below
+    """
+    config = {}
+    config['grid file'] = grid_file
+    config['fields file'] = data_file
+    config['variable'] = variable
+    config['levels'] = levels
+    config['bounds'] = bounds
+    config['colormap'] = colormap
+    config['lats'] = lats
+    config['comout'] = comout
+    config['max depth'] = 5000.0
+    return config
+
+
+def plot_horizontal_slice(config):
+    """
+    pcolormesh of a horizontal slice of an ocean field
+    """
+    level = config['levels'][0]
+    grid = xr.open_dataset(config['grid file'])
+    data = xr.open_dataset(config['fields file'])
+    if config['variable'] in ['Temp', 'Salt', 'u', 'v']:
+        slice_data = np.squeeze(data[config['variable']])[level, :, :]
+    else:
+        slice_data = np.squeeze(data[config['variable']])
+    bounds = config['bounds']
+    fig, ax = plt.subplots(figsize=(8, 5), subplot_kw={'projection': ccrs.PlateCarree()})
+    plt.pcolormesh(np.squeeze(grid.lon),
+                   np.squeeze(grid.lat),
+                   slice_data,
+                   vmin=bounds[0], vmax=bounds[1],
+                   transform=ccrs.PlateCarree(),
+                   cmap=config['colormap'])
+    plt.colorbar(label=config['variable']+' Level '+str(level), shrink=0.5, orientation='horizontal')
+    # ax.coastlines()  # TODO: make this work on hpc
+    ax.gridlines(draw_labels=True)
+    # ax.add_feature(cartopy.feature.LAND)  # TODO: make this work on hpc
+    dirname = os.path.join(config['comout'], config['variable'])
+    os.makedirs(dirname, exist_ok=True)
+    figname = os.path.join(dirname, config['variable']+'_Level_'+str(level))
+    plt.savefig(figname, bbox_inches='tight')
+
+
+def plot_zonal_slice(config):
+    """
+    pcolormesh of a zonal slice of an ocean field
+    """
+    lat = float(config['lats'][0])
+    grid = xr.open_dataset(config['grid file'])
+    data = xr.open_dataset(config['fields file'])
+    lat_index = np.argmin(np.array(np.abs(np.squeeze(grid.lat)[:, 0]-lat)))
+    slice_data = np.squeeze(np.array(data[config['variable']].sel(yaxis_1=lat_index)))
+    depth = np.squeeze(np.array(grid['h'].sel(yaxis_1=lat_index)))
+    depth[np.where(np.abs(depth) > 10000.0)] = 0.0
+    depth = np.cumsum(depth, axis=0)
+    bounds = config['bounds']
+    x = np.tile(np.squeeze(grid.lon[:, lat_index]), (np.shape(depth)[0], 1))
+    fig, ax = plt.subplots(figsize=(8, 5))
+    plt.pcolormesh(x, -depth, slice_data,
+                   vmin=bounds[0], vmax=bounds[1],
+                   cmap=config['colormap'])
+    plt.colorbar(label=config['variable']+' Lat '+str(lat), shrink=0.5, orientation='horizontal')
+    ax.set_ylim(-config['max depth'], 0)
+    dirname = os.path.join(config['comout'], config['variable'])
+    os.makedirs(dirname, exist_ok=True)
+    figname = os.path.join(dirname, config['variable'] +
+                           'zonal_lat_'+str(int(lat)) + '_' + str(int(config['max depth'])) + 'm')
+    plt.savefig(figname, bbox_inches='tight')
+
+
+comout = os.getenv('COMOUT')
+cyc = os.getenv('cyc')
+bcyc = str((int(cyc) - 3) % 24)
+
+# TODO: do not write to COM, instead dump figures in DATA and copy to COM.
+grid_file = os.path.join(comout, 'gdas.t'+bcyc+'z.ocngrid.nc')
+
+#######################################
+# INCREMENT
+#######################################
+incr_cmap = 'RdBu'
+data_file = os.path.join(comout, 'gdas.t'+cyc+'z.ocninc.nc')
+config = plot_config(grid_file=grid_file,
+                     data_file=data_file,
+                     colormap=incr_cmap,
+                     comout=os.path.join(comout, 'vrfy', 'incr'))
+
+#######################################
+# zonal slices
+
+for lat in np.arange(-60, 60, 10):
+
+    for max_depth in [700.0, 5000.0]:
+        config['lats'] = [lat]
+        config['max depth'] = max_depth
+
+        # Temperature
+        config.update({'variable': 'Temp', 'levels': [1], 'bounds': [-.5, .5]})
+        plot_zonal_slice(config)
+
+        # Salinity
+        config.update({'variable': 'Salt', 'levels': [1], 'bounds': [-.1, .1]})
+        plot_zonal_slice(config)
+
+#######################################
+# Horizontal slices
+
+# Temperature
+config.update({'variable': 'Temp', 'levels': [1], 'bounds': [-1, 1]})
+plot_horizontal_slice(config)
+
+# Salinity
+config.update({'variable': 'Salt', 'bounds': [-0.1, 0.1]})
+plot_horizontal_slice(config)
+
+# Sea surface height
+config.update({'variable': 'ave_ssh', 'bounds': [-0.1, 0.1]})
+plot_horizontal_slice(config)
+
+#######################################
+# Std Bkg. Error
+#######################################
+bmat_cmap = 'jet'
+data_file = os.path.join(comout, 'gdas.t'+cyc+'z.ocn.bkgerr_stddev.nc')
+config = plot_config(grid_file=grid_file,
+                     data_file=data_file,
+                     colormap=bmat_cmap,
+                     comout=os.path.join(comout, 'vrfy', 'bkgerr'))
+
+#######################################
+# zonal slices
+
+for lat in np.arange(-60, 60, 10):
+
+    for max_depth in [700.0, 5000.0]:
+        config['lats'] = [lat]
+        config['max depth'] = max_depth
+
+        # Temperature
+        config.update({'variable': 'Temp', 'levels': [1], 'bounds': [0, 1.5]})
+        plot_zonal_slice(config)
+
+        # Salinity
+        config.update({'variable': 'Salt', 'levels': [1], 'bounds': [0, .2]})
+        plot_zonal_slice(config)
+
+#######################################
+# Horizontal slices
+
+# Temperature
+config.update({'variable': 'Temp', 'levels': [1], 'bounds': [0, 2]})
+plot_horizontal_slice(config)
+
+# Salinity
+config.update({'variable': 'Salt', 'bounds': [0, 0.2]})
+plot_horizontal_slice(config)
+
+# Sea surface height
+config.update({'variable': 'ave_ssh', 'bounds': [0, 0.1]})
+plot_horizontal_slice(config)
diff --git a/test/soca/gw/CMakeLists.txt b/test/soca/gw/CMakeLists.txt
index 00e232e83..2c89adc22 100644
--- a/test/soca/gw/CMakeLists.txt
+++ b/test/soca/gw/CMakeLists.txt
@@ -21,7 +21,8 @@ add_test(NAME test_gdasapp_soca_run_clean
 set(jjob_list "JGDAS_GLOBAL_OCEAN_ANALYSIS_PREP"
               "JGDAS_GLOBAL_OCEAN_ANALYSIS_BMAT"
               "JGDAS_GLOBAL_OCEAN_ANALYSIS_RUN"
-              "JGDAS_GLOBAL_OCEAN_ANALYSIS_POST")
+              "JGDAS_GLOBAL_OCEAN_ANALYSIS_POST"
+              "JGDAS_GLOBAL_OCEAN_ANALYSIS_VRFY")
 
 set(setup "")
 foreach(jjob ${jjob_list})
diff --git a/test/soca/test_run.sh b/test/soca/test_run.sh
index 7c9dd6b9a..ba653433b 100755
--- a/test/soca/test_run.sh
+++ b/test/soca/test_run.sh
@@ -35,4 +35,4 @@ test_file  ${DATA}/Data/ocn.iter1.incr.2018-04-15T09:00:00Z.nc
 test_file  ${DATA}/Data/ocn.3dvarfgat_pseudo.an.2018-04-15T09:00:00Z.nc
 test_file  ${DATA}/Data/ocn.3dvarfgat_pseudo.an.2018-04-15T12:00:00Z.nc
 test_file  ${DATA}/Data/ocn.3dvarfgat_pseudo.an.2018-04-15T15:00:00Z.nc
-test_file  ${COMOUT}/inc.nc
+test_file  ${COMOUT}/ocnanal_12/inc.nc
diff --git a/ush/soca/run_jjobs.py b/ush/soca/run_jjobs.py
index dcb5db6ca..d69dfb348 100755
--- a/ush/soca/run_jjobs.py
+++ b/ush/soca/run_jjobs.py
@@ -7,6 +7,11 @@
 
 
 machines = {"container", "hera", "orion"}
+MODS = {'JGDAS_GLOBAL_OCEAN_ANALYSIS_PREP': 'GDAS',
+        'JGDAS_GLOBAL_OCEAN_ANALYSIS_BMAT': 'GDAS',
+        'JGDAS_GLOBAL_OCEAN_ANALYSIS_RUN': 'GDAS',
+        'JGDAS_GLOBAL_OCEAN_ANALYSIS_POST': 'GDAS',
+        'JGDAS_GLOBAL_OCEAN_ANALYSIS_VRFY': 'EVA'}
 
 
 class JobCard:
@@ -112,14 +117,14 @@ def close(self):
         self.f.close()
         subprocess.run(["chmod", "+x", self.name])
 
-    def modules(self):
+    def _modules(self, jjob):
         """
         Write a section that will load the machine dependent modules
         """
         if self.machine != "container":
             self.f.write("module purge \n")
             self.f.write("module use ${HOMEgfs}/sorc/gdas.cd/modulefiles \n")
-            self.f.write(f"module load GDAS/{self.machine} \n")
+            self.f.write(f"module load {MODS[jjob]}/{self.machine} \n")
 
     def copy_bkgs(self):
         """
@@ -161,6 +166,7 @@ def jjobs(self):
         """
         runjobs = "# Run jjobs\n"
         for job in self.config['jjobs']:
+            self._modules(job)  # Add module's jjob
             thejob = "${HOMEgfs}/jobs/"+job
             runjobs += f"{thejob} &>{job}.out\n"
         self.f.write(runjobs)
@@ -208,7 +214,6 @@ def main():
     run_card.fixconfigs()                # over-write some of the config variables
     run_card.header()                    # prepare a machine dependent header (SLURM or nothing)
     run_card.export_env_vars_script()
-    run_card.modules()
     run_card.copy_bkgs()
     run_card.jjobs()
     run_card.close()