diff --git a/jobs/JLANDDA_PLOT_STATS b/jobs/JLANDDA_PLOT_STATS
index f0a83a74..73fb19b6 100755
--- a/jobs/JLANDDA_PLOT_STATS
+++ b/jobs/JLANDDA_PLOT_STATS
@@ -80,6 +80,8 @@ mkdir -p ${DATA_SHARE}
 export DATA_HOFX="${DATA_HOFX:-${DATAROOT}/DATA_SHARE/hofx}"
 mkdir -p ${DATA_HOFX}
 
+# Set other dates
+export NTIME=$($NDATE 24 $PDY$cyc)
 #
 #-----------------------------------------------------------------------
 #
diff --git a/parm/land_analysis_hera.yaml b/parm/land_analysis_hera.yaml
index c58f83f9..dcebb66c 100644
--- a/parm/land_analysis_hera.yaml
+++ b/parm/land_analysis_hera.yaml
@@ -187,33 +187,6 @@ workflow:
         taskdep:
           attrs:
             task: analysis
-    task_plot_stats:
-      attrs:
-        cycledefs: cycled
-        maxtries: 2
-      envars:
-        MACHINE: "&MACHINE;"
-        SCHED: "&SCHED;"
-        ACCOUNT: "&ACCOUNT;"
-        model_ver: "&model_ver;"
-        RUN: "&RUN;"
-        HOMElandda: "&HOMElandda;"
-        COMROOT: "&COMROOT;"
-        DATAROOT: "&DATAROOT;"
-        KEEPDATA: "&KEEPDATA;"
-        PDY: "&PDY;"
-        cyc: "&cyc;"
-      account: "&ACCOUNT;"
-      command: '&HOMElandda;/parm/task_load_modules_run_jjob.sh "plot_stats" "&HOMElandda;" "&MACHINE;"'
-      jobname: plot_stats
-      cores: 1
-      walltime: 00:10:00
-      queue: batch
-      join: "&LOGDIR;/plot_stats&LOGFN_SUFFIX;"
-      dependency:
-        taskdep:
-          attrs:
-            task: analysis
     task_forecast:
       attrs:
         cycledefs: cycled
@@ -252,3 +225,31 @@ workflow:
         taskdep:
           attrs:
             task: post_anal
+    task_plot_stats:
+      attrs:
+        cycledefs: cycled
+        maxtries: 2
+      envars:
+        MACHINE: "&MACHINE;"
+        SCHED: "&SCHED;"
+        ACCOUNT: "&ACCOUNT;"
+        model_ver: "&model_ver;"
+        RUN: "&RUN;"
+        HOMElandda: "&HOMElandda;"
+        COMROOT: "&COMROOT;"
+        DATAROOT: "&DATAROOT;"
+        KEEPDATA: "&KEEPDATA;"
+        PDY: "&PDY;"
+        cyc: "&cyc;"
+      account: "&ACCOUNT;"
+      command: '&HOMElandda;/parm/task_load_modules_run_jjob.sh "plot_stats" "&HOMElandda;" "&MACHINE;"'
+      jobname: plot_stats
+      cores: 1
+      walltime: 00:10:00
+      queue: batch
+      join: "&LOGDIR;/plot_stats&LOGFN_SUFFIX;"
+      dependency:
+        taskdep:
+          attrs:
+            task: forecast
+
diff --git a/parm/land_analysis_hercules.yaml b/parm/land_analysis_hercules.yaml
index e4a47e2a..e08902d9 100644
--- a/parm/land_analysis_hercules.yaml
+++ b/parm/land_analysis_hercules.yaml
@@ -9,7 +9,7 @@ workflow:
         group: cycled
       spec: 200001030000 200001040000 24:00:00
   entities:
-    MACHINE: "orion"
+    MACHINE: "hercules"
     SCHED: "slurm"
     ACCOUNT: "epic"
     EXP_BASEDIR: "/work2/noaa/epic/{USER}/landda_test"
@@ -187,33 +187,6 @@ workflow:
         taskdep:
           attrs:
             task: analysis
-    task_plot_stats:
-      attrs:
-        cycledefs: cycled
-        maxtries: 2
-      envars:
-        MACHINE: "&MACHINE;"
-        SCHED: "&SCHED;"
-        ACCOUNT: "&ACCOUNT;"
-        model_ver: "&model_ver;"
-        RUN: "&RUN;"
-        HOMElandda: "&HOMElandda;"
-        COMROOT: "&COMROOT;"
-        DATAROOT: "&DATAROOT;"
-        KEEPDATA: "&KEEPDATA;"
-        PDY: "&PDY;"
-        cyc: "&cyc;"
-      account: "&ACCOUNT;"
-      command: '&HOMElandda;/parm/task_load_modules_run_jjob.sh "plot_stats" "&HOMElandda;" "&MACHINE;"'
-      jobname: plot_stats
-      cores: 1
-      walltime: 00:10:00
-      queue: batch
-      join: "&LOGDIR;/plot_stats&LOGFN_SUFFIX;"
-      dependency:
-        taskdep:
-          attrs:
-            task: analysis
     task_forecast:
       attrs:
         cycledefs: cycled
@@ -252,3 +225,31 @@ workflow:
         taskdep:
           attrs:
             task: post_anal
+    task_plot_stats:
+      attrs:
+        cycledefs: cycled
+        maxtries: 2
+      envars:
+        MACHINE: "&MACHINE;"
+        SCHED: "&SCHED;"
+        ACCOUNT: "&ACCOUNT;"
+        model_ver: "&model_ver;"
+        RUN: "&RUN;"
+        HOMElandda: "&HOMElandda;"
+        COMROOT: "&COMROOT;"
+        DATAROOT: "&DATAROOT;"
+        KEEPDATA: "&KEEPDATA;"
+        PDY: "&PDY;"
+        cyc: "&cyc;"
+      account: "&ACCOUNT;"
+      command: '&HOMElandda;/parm/task_load_modules_run_jjob.sh "plot_stats" "&HOMElandda;" "&MACHINE;"'
+      jobname: plot_stats
+      cores: 1
+      walltime: 00:10:00
+      queue: batch
+      join: "&LOGDIR;/plot_stats&LOGFN_SUFFIX;"
+      dependency:
+        taskdep:
+          attrs:
+            task: forecast
+
diff --git a/parm/land_analysis_orion.yaml b/parm/land_analysis_orion.yaml
index 575f9c3e..444ced28 100644
--- a/parm/land_analysis_orion.yaml
+++ b/parm/land_analysis_orion.yaml
@@ -187,33 +187,6 @@ workflow:
         taskdep:
           attrs:
             task: analysis
-    task_plot_stats:
-      attrs:
-        cycledefs: cycled
-        maxtries: 2
-      envars:
-        MACHINE: "&MACHINE;"
-        SCHED: "&SCHED;"
-        ACCOUNT: "&ACCOUNT;"
-        model_ver: "&model_ver;"
-        RUN: "&RUN;"
-        HOMElandda: "&HOMElandda;"
-        COMROOT: "&COMROOT;"
-        DATAROOT: "&DATAROOT;"
-        KEEPDATA: "&KEEPDATA;"
-        PDY: "&PDY;"
-        cyc: "&cyc;"
-      account: "&ACCOUNT;"
-      command: '&HOMElandda;/parm/task_load_modules_run_jjob.sh "plot_stats" "&HOMElandda;" "&MACHINE;"'
-      jobname: plot_stats
-      cores: 1
-      walltime: 00:10:00
-      queue: batch
-      join: "&LOGDIR;/plot_stats&LOGFN_SUFFIX;"
-      dependency:
-        taskdep:
-          attrs:
-            task: analysis
     task_forecast:
       attrs:
         cycledefs: cycled
@@ -252,3 +225,31 @@ workflow:
         taskdep:
           attrs:
             task: post_anal
+    task_plot_stats:
+      attrs:
+        cycledefs: cycled
+        maxtries: 2
+      envars:
+        MACHINE: "&MACHINE;"
+        SCHED: "&SCHED;"
+        ACCOUNT: "&ACCOUNT;"
+        model_ver: "&model_ver;"
+        RUN: "&RUN;"
+        HOMElandda: "&HOMElandda;"
+        COMROOT: "&COMROOT;"
+        DATAROOT: "&DATAROOT;"
+        KEEPDATA: "&KEEPDATA;"
+        PDY: "&PDY;"
+        cyc: "&cyc;"
+      account: "&ACCOUNT;"
+      command: '&HOMElandda;/parm/task_load_modules_run_jjob.sh "plot_stats" "&HOMElandda;" "&MACHINE;"'
+      jobname: plot_stats
+      cores: 1
+      walltime: 00:10:00
+      queue: batch
+      join: "&LOGDIR;/plot_stats&LOGFN_SUFFIX;"
+      dependency:
+        taskdep:
+          attrs:
+            task: forecast
+
diff --git a/parm/templates/template.plot_restart.yaml b/parm/templates/template.plot_restart.yaml
new file mode 100644
index 00000000..d1def18e
--- /dev/null
+++ b/parm/templates/template.plot_restart.yaml
@@ -0,0 +1,8 @@
+path_data: 'PATH_DATA'
+work_dir: 'WORK_DIR'
+fn_data_base: 'XXFN_DATA_BASE'
+fn_data_ext: 'XXFN_DATA_EXT'
+soil_lvl_number: 'XXSOIL_LEVEL_NUMBER'
+out_title_base: 'XXOUT_TITLE_BASE'
+out_fn_base: 'XXOUT_FN_BASE'
+machine: 'XXMACHINE'
diff --git a/scripts/exlandda_plot_stats.sh b/scripts/exlandda_plot_stats.sh
index 208f32bb..a8203fbc 100755
--- a/scripts/exlandda_plot_stats.sh
+++ b/scripts/exlandda_plot_stats.sh
@@ -2,11 +2,16 @@
 
 set -xue
 
+YYYY=${PDY:0:4}
+MM=${PDY:4:2}
+DD=${PDY:6:2}
+HH=${cyc}
+
 ############################################################
-# Scatter Plot
+# Stats Plot
 ############################################################
 
-cp ${PARMlandda}/templates/template.plot_hofx.yaml plot_hofx.yaml
+cp -p ${PARMlandda}/templates/template.plot_hofx.yaml plot_hofx.yaml
 
 # Path to the directory containing the input file
 INPUTFP="${DATA_HOFX}"
@@ -36,9 +41,48 @@ sed -i -e "s/XXMACHINE/${MACHINE}/g" plot_hofx.yaml
 
 ${USHlandda}/hofx_analysis_stats.py
 if [[ $? != 0 ]]; then
-  echo "Scatter/Histogram plots failed"
+  echo "FATAL ERROR: Scatter/Histogram plots failed"
   exit 33
 fi
 
 # Copy result files to COMOUT
 cp -p ${PREOUTFN}* ${COMOUTplot}
+
+
+############################################################
+# Restart Plot
+############################################################
+
+nYYYY=${NTIME:0:4}
+nMM=${NTIME:4:2}
+nDD=${NTIME:6:2}
+nHH=${NTIME:8:2}
+
+cp -p ${PARMlandda}/templates/template.plot_restart.yaml plot_restart.yaml
+
+PATH_DATA="${COMIN}"
+WORK_DIR="${DATA}"
+FN_DATA_BASE="ufs_land_restart.${nYYYY}-${nMM}-${nDD}_${nHH}-00-00.tile"
+FN_DATA_EXT=".nc"
+SOIL_LEVEL_NUMBER="1"
+OUT_TITLE_BASE="Land-DA::restart::${nYYYY}-${nMM}-${nDD}_${nHH}::"
+OUT_FN_BASE="landda_out_restart_${nYYYY}-${nMM}-${nDD}_${nHH}_"
+
+sed -i "s|PATH_DATA|${PATH_DATA}|g" plot_restart.yaml
+sed -i "s|WORK_DIR|${WORK_DIR}|g" plot_restart.yaml
+sed -i -e "s/XXFN_DATA_BASE/${FN_DATA_BASE}/g" plot_restart.yaml
+sed -i -e "s/XXFN_DATA_EXT/${FN_DATA_EXT}/g" plot_restart.yaml
+sed -i -e "s/XXSOIL_LEVEL_NUMBER/${SOIL_LEVEL_NUMBER}/g" plot_restart.yaml
+sed -i -e "s/XXOUT_TITLE_BASE/${OUT_TITLE_BASE}/g" plot_restart.yaml
+sed -i -e "s/XXOUT_FN_BASE/${OUT_FN_BASE}/g" plot_restart.yaml
+sed -i -e "s/XXMACHINE/${MACHINE}/g" plot_restart.yaml
+
+${USHlandda}/plot_forecast_restart.py
+if [[ $? != 0 ]]; then
+  echo "FATAL ERROR: Forecast restart plots failed"
+  exit 44
+fi
+
+# Copy result files to COMOUT
+cp -p ${OUT_FN_BASE}* ${COMOUTplot}
+
diff --git a/ush/plot_forecast_restart.py b/ush/plot_forecast_restart.py
new file mode 100755
index 00000000..9f0e91b9
--- /dev/null
+++ b/ush/plot_forecast_restart.py
@@ -0,0 +1,240 @@
+#!/usr/bin/env python3
+
+###################################################################### CHJ #####
+## Name		: plot_forecast_restart.py
+## Usage	: Plot restart output file of land-DA workflow
+## Input files  : ufs_land_restart.tile#.nc
+## NOAA/EPIC
+## History ===============================
+## V000: 2024/09/26: Chan-Hoo Jeon : Preliminary version
+## V001: 2024/10/05: Chan-Hoo Jeon : Add to land-DA workflow
+###################################################################### CHJ #####
+
+import os, sys
+import yaml
+import numpy as np
+import netCDF4 as nc
+import cartopy
+import cartopy.crs as ccrs
+import cartopy.feature as cfeature
+import matplotlib.pyplot as plt
+import matplotlib.colors as colors
+import matplotlib.ticker
+import matplotlib as mpl
+from matplotlib.colors import ListedColormap
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+
+# Main part (will be called at the end) ============================= CHJ =====
+def main():
+# =================================================================== CHJ =====
+
+    global num_tiles
+
+    yaml_file="plot_restart.yaml"
+    with open(yaml_file, 'r') as f:
+        yaml_data=yaml.load(f, Loader=yaml.FullLoader)
+    f.close()
+    print("YAML_DATA:",yaml_data)
+
+    path_data=yaml_data['path_data']
+    work_dir=yaml_data['work_dir']
+    fn_data_base=yaml_data['fn_data_base']
+    fn_data_ext=yaml_data['fn_data_ext']
+    soil_lvl_num=yaml_data['soil_lvl_number']
+    out_title_base=yaml_data['out_title_base']
+    out_fn_base=yaml_data['out_fn_base']
+    machine=yaml_data['machine']
+
+    if yaml_data['machine']=='hera':
+        cartopy.config['data_dir']='/scratch2/NAGAPE/epic/UFS_Land-DA_Dev/inputs/NaturalEarth'
+    elif yaml_data['machine']=='orion' or yaml_data['machine']=='hercules':
+        cartopy.config['data_dir']='/work/noaa/epic/UFS_Land-DA_Dev/inputs/NaturalEarth'
+
+    #var_list=["snwdph","smc"]
+    var_list=["snwdph"]
+    # Number of tiles
+    num_tiles=6
+
+    # get lon, lat
+    get_geo(path_data,fn_data_base,fn_data_ext)
+    # plot restart file
+    for var_nm in var_list:
+        plot_data(path_data,fn_data_base,fn_data_ext,var_nm,soil_lvl_num,out_title_base,out_fn_base,work_dir)
+       
+
+# geo lon/lat from orography ======================================== CHJ =====
+def get_geo(path_data,fn_data_base,fn_data_ext):
+# =================================================================== CHJ =====
+
+    global glon,glat
+
+    print(' ===== geo data files ====================================')
+
+    # open the data file
+    for it in range(num_tiles):
+        itp=it+1
+        fn_data=fn_data_base+str(itp)+fn_data_ext
+        fp_data=os.path.join(path_data,fn_data)
+        try: data_raw=nc.Dataset(fp_data)
+        except: raise Exception('Could NOT find the file',fp_data)
+        if itp == 1:
+            print(data_raw)
+        # Extract geo data
+        glon_data=np.ma.masked_invalid(data_raw.variables['grid_xt'])
+        print('Dimension of glon(grid_xt)=',glon_data.shape)
+        print('Tile',itp,',max:',np.max(glon_data))
+        print('Tile',itp,',min:',np.min(glon_data))
+
+        glat_data=np.ma.masked_invalid(data_raw.variables['grid_yt'])
+        print('Dimension of glat(grid_yt)=',glat_data.shape)
+        print('Tile',itp,',max:',np.max(glat_data))
+        print('Tile',itp,',min:',np.min(glat_data))
+
+        if itp == 1:
+            ny,nx=glon_data.shape
+            glon=np.zeros((num_tiles,ny,nx))
+            glat=np.zeros((num_tiles,ny,nx))
+
+        data_raw.close()        
+        glon[it,:,:]=glon_data[:,:]
+        glat[it,:,:]=glat_data[:,:]
+
+    print('Dimension of glon=',glon.shape)
+    print('Dimension of glon=',glat.shape)
+
+
+# Get sfc_data from files and plot ================================== CHJ =====
+def plot_data(path_data,fn_data_base,fn_data_ext,var_nm,soil_lvl_num,out_title_base,out_fn_base,work_dir):
+# =================================================================== CHJ =====
+
+    # center of map
+    c_lon=-77.0369
+
+    print(' ===== data file: '+var_nm+' ========================')
+    # open the data file
+    for it in range(num_tiles):
+        itp=it+1
+        fn_data=fn_data_base+str(itp)+fn_data_ext
+        fp_data=os.path.join(path_data,fn_data)
+        try: data_raw=nc.Dataset(fp_data)
+        except: raise Exception('Could NOT find the file',fp_data)
+        # Extract valid variable
+        var_data=np.ma.masked_invalid(data_raw.variables[var_nm])
+        if var_nm == 'stc' or var_nm == 'smc' or var_nm == 'slc':
+            print('Dimension of original data=',var_data.shape)
+            var_data_2d=var_data[:,soil_lvl_num-1,:,:]
+        else:
+            var_data_2d=var_data                
+ 
+        print('Dimension of data=',var_data_2d.shape)
+        print('Tile',itp,',max:',np.max(var_data_2d))
+        print('Tile',itp,',min:',np.min(var_data_2d))
+
+        if itp == 1:
+            plt_var=var_data_2d
+        else:
+            plt_var=np.ma.concatenate((plt_var,var_data_2d),axis=0)
+        data_raw.close()
+
+#    plt_var=np.vstack(var_data_all)
+    print('Dimension of data set=',plt_var.shape)
+
+    cs_max=np.max(plt_var)
+    cs_min=np.min(plt_var)
+    print('cs_max=',cs_max)
+    print('cs_min=',cs_min)
+
+    cs_cmap='gist_ncar_r'
+    cbar_extend='neither'
+
+    for it in range(num_tiles):
+        itp=it+1
+        glon_tile=np.squeeze(glon[it,:,:])
+        if itp == 1:
+            glon_tile=(glon_tile+180)%360-180
+        glat_tile=np.squeeze(glat[it,:,:])
+        var_tile=np.squeeze(plt_var[it,:,:])
+        c_glon=np.round(np.mean(glon_tile),decimals=2)
+        c_glat=np.round(np.mean(glat_tile),decimals=2)
+        print("c_glon, c_glat for tile",str(it+1),"=",c_glon,c_glat)
+
+        out_title=out_title_base+var_nm+'::Tile'+str(itp)
+        out_fn=out_fn_base+var_nm+'_tile'+str(itp)
+
+        fig,ax=plt.subplots(1,1,subplot_kw=dict(projection=ccrs.Orthographic(c_glon,c_glat)))
+        ax.set_title(out_title, fontsize=6)
+        # Call background plot
+        back_plot(ax)
+
+        cs=ax.pcolormesh(glon_tile,glat_tile,var_tile,cmap=cs_cmap,
+            rasterized=True,vmin=cs_min,vmax=cs_max,transform=ccrs.PlateCarree())
+        divider=make_axes_locatable(ax)
+        ax_cb=divider.new_horizontal(size="3%",pad=0.1,axes_class=plt.Axes)
+        fig.add_axes(ax_cb)
+        cbar=plt.colorbar(cs,cax=ax_cb,extend='neither')
+        cbar.ax.tick_params(labelsize=6)
+        cbar.set_label(var_nm,fontsize=6)
+
+#    for it in range(num_tiles):
+#        cs=ax.pcolormesh(glon[it,:,:],glat[it,:,:],plt_var[it,:,:],cmap=cs_cmap,rasterized=True,
+#            vmin=cs_min,vmax=cs_max,transform=ccrs.PlateCarree())
+#    divider=make_axes_locatable(ax)
+#    ax_cb=divider.new_horizontal(size="3%",pad=0.1,axes_class=plt.Axes)
+#    fig.add_axes(ax_cb)
+#    cbar=plt.colorbar(cs,cax=ax_cb,extend=cbar_extend)
+#    cbar.ax.tick_params(labelsize=6)
+#    cbar.set_label(var_nm,fontsize=6)
+
+        # Output figure
+        ndpi=300
+        out_file(work_dir,out_fn,ndpi)
+
+
+# Background plot ==================================================== CHJ =====
+def back_plot(ax):
+# ==================================================================== CHJ =====
+
+    # Resolution of background natural earth data ('50m' or '110m')
+    back_res='50m'
+
+    fline_wd=0.5  # line width
+    falpha=0.7 # transparency
+
+    # natural_earth
+    land=cfeature.NaturalEarthFeature('physical','land',back_res,
+                      edgecolor='face',facecolor=cfeature.COLORS['land'],
+                      alpha=falpha)
+    lakes=cfeature.NaturalEarthFeature('physical','lakes',back_res,
+                      edgecolor='blue',facecolor='none',
+                      linewidth=fline_wd,alpha=falpha)
+    coastline=cfeature.NaturalEarthFeature('physical','coastline',
+                      back_res,edgecolor='black',facecolor='none',
+                      linewidth=fline_wd,alpha=falpha)
+    states=cfeature.NaturalEarthFeature('cultural','admin_1_states_provinces',
+                      back_res,edgecolor='green',facecolor='none',
+                      linewidth=fline_wd,linestyle=':',alpha=falpha)
+    borders=cfeature.NaturalEarthFeature('cultural','admin_0_countries',
+                      back_res,edgecolor='red',facecolor='none',
+                      linewidth=fline_wd,alpha=falpha)
+
+#    ax.add_feature(land)
+#    ax.add_feature(lakes)
+#    ax.add_feature(states)
+#    ax.add_feature(borders)
+    ax.add_feature(coastline)
+
+
+# Output file ======================================================= CHJ =====
+def out_file(work_dir,out_file,ndpi):
+# =================================================================== CHJ =====
+    # Output figure
+    fp_out=os.path.join(work_dir,out_file)
+    plt.savefig(fp_out+'.png',dpi=ndpi,bbox_inches='tight')
+    plt.close('all')
+
+
+# Main call ========================================================= CHJ =====
+if __name__=='__main__':
+    main()
+