From 5f5b15284be47c790bb6936cff077b29a3c08516 Mon Sep 17 00:00:00 2001
From: "Chan-Hoo.Jeon-NOAA" <60152248+chan-hoo@users.noreply.github.com>
Date: Tue, 22 Oct 2024 13:59:10 -0400
Subject: [PATCH] Add time-history plots for stats (#151)

* add python for plotting time-history

* update scripts

* fix wrong parms

* fix incorrect var names

* fix missing env

* remove legend from plot

* add nobs plot

* change we2e_vav to we2e_test

* remove templates for plot and set yaml with cat eof
---
 modulefiles/tasks/hera/task.plot_stats.lua |   1 +
 parm/land_analysis_hera.yaml               |  11 +-
 parm/land_analysis_hercules.yaml           |  11 +-
 parm/land_analysis_orion.yaml              |  11 +-
 parm/templates/template.plot_hofx.yaml     |   8 -
 parm/templates/template.plot_restart.yaml  |   8 -
 scripts/exlandda_analysis.sh               |   4 +-
 scripts/exlandda_forecast.sh               |   4 +-
 scripts/exlandda_plot_stats.sh             |  90 +++++--
 scripts/exlandda_post_anal.sh              |   4 +-
 ush/plot_analysis_timehistory.py           | 294 +++++++++++++++++++++
 versions/run.ver_hera                      |   1 +
 versions/run.ver_hercules                  |   1 +
 versions/run.ver_orion                     |   1 +
 14 files changed, 388 insertions(+), 61 deletions(-)
 delete mode 100644 parm/templates/template.plot_hofx.yaml
 delete mode 100644 parm/templates/template.plot_restart.yaml
 create mode 100755 ush/plot_analysis_timehistory.py

diff --git a/modulefiles/tasks/hera/task.plot_stats.lua b/modulefiles/tasks/hera/task.plot_stats.lua
index 1b1dd4e0..612429da 100644
--- a/modulefiles/tasks/hera/task.plot_stats.lua
+++ b/modulefiles/tasks/hera/task.plot_stats.lua
@@ -9,6 +9,7 @@ load(pathJoin("py-cartopy", py_cartopy_ver))
 load(pathJoin("py-matplotlib", py_matplotlib_ver))
 load(pathJoin("py-netcdf4", py_netcdf4_ver))
 load(pathJoin("py-numpy", py_numpy_ver))
+load(pathJoin("py-pandas", py_pandas_ver))
 load(pathJoin("py-pyyaml", py_pyyaml_ver))
 load(pathJoin("py-scipy", py_scipy_ver))
 load(pathJoin("py-xarray", py_xarray_ver))
diff --git a/parm/land_analysis_hera.yaml b/parm/land_analysis_hera.yaml
index f9de5635..2abffbb2 100644
--- a/parm/land_analysis_hera.yaml
+++ b/parm/land_analysis_hera.yaml
@@ -34,7 +34,7 @@ workflow:
     OBS_TYPES: "GHCN"
     DAtype: "letkfoi_snow"
     TSTUB: "oro_C96.mx100"
-    WE2E_VAV: "YES"
+    WE2E_TEST: "NO"
     WE2E_ATOL: "1e-7"
     WE2E_LOG_FN: "we2e.log"
     NET: "landda"
@@ -137,7 +137,7 @@ workflow:
         ACCOUNT: "&ACCOUNT;"
         RES: "&RES;"
         TSTUB: "&TSTUB;"
-        WE2E_VAV: "&WE2E_VAV;"
+        WE2E_TEST: "&WE2E_TEST;"
         WE2E_ATOL: "&WE2E_ATOL;"
         WE2E_LOG_FN: "&WE2E_LOG_FN;"
         LOGDIR: "&LOGDIR;"
@@ -172,7 +172,7 @@ workflow:
         ACCOUNT: "&ACCOUNT;"
         RES: "&RES;"
         TSTUB: "&TSTUB;"
-        WE2E_VAV: "&WE2E_VAV;"
+        WE2E_TEST: "&WE2E_TEST;"
         WE2E_ATOL: "&WE2E_ATOL;"
         WE2E_LOG_FN: "&WE2E_LOG_FN;"
         LOGDIR: "&LOGDIR;"
@@ -208,7 +208,7 @@ workflow:
         ATMOS_FORC: "&ATMOS_FORC;"
         RES: "&RES;"
         WARMSTART_DIR: "&WARMSTART_DIR;"
-        WE2E_VAV: "&WE2E_VAV;"
+        WE2E_TEST: "&WE2E_TEST;"
         WE2E_ATOL: "&WE2E_ATOL;"
         WE2E_LOG_FN: "&WE2E_LOG_FN;"
         LOGDIR: "&LOGDIR;"
@@ -259,6 +259,9 @@ workflow:
         KEEPDATA: "&KEEPDATA;"
         PDY: "&PDY;"
         cyc: "&cyc;"
+        LOGDIR: "&LOGDIR;"
+        NPROCS_ANALYSIS: "&NPROCS_ANALYSIS;"
+        NPROCS_FORECAST: "&NPROCS_FORECAST;"
       account: "&ACCOUNT;"
       command: '&HOMElandda;/parm/task_load_modules_run_jjob.sh "plot_stats" "&HOMElandda;" "&MACHINE;"'
       jobname: plot_stats
diff --git a/parm/land_analysis_hercules.yaml b/parm/land_analysis_hercules.yaml
index a39e22ad..27a729cc 100644
--- a/parm/land_analysis_hercules.yaml
+++ b/parm/land_analysis_hercules.yaml
@@ -34,7 +34,7 @@ workflow:
     OBS_TYPES: "GHCN"
     DAtype: "letkfoi_snow"
     TSTUB: "oro_C96.mx100"
-    WE2E_VAV: "YES"
+    WE2E_TEST: "NO"
     WE2E_ATOL: "1e-7"
     WE2E_LOG_FN: "we2e.log"
     NET: "landda"
@@ -137,7 +137,7 @@ workflow:
         ACCOUNT: "&ACCOUNT;"
         RES: "&RES;"
         TSTUB: "&TSTUB;"
-        WE2E_VAV: "&WE2E_VAV;"
+        WE2E_TEST: "&WE2E_TEST;"
         WE2E_ATOL: "&WE2E_ATOL;"
         WE2E_LOG_FN: "&WE2E_LOG_FN;"
         LOGDIR: "&LOGDIR;"
@@ -172,7 +172,7 @@ workflow:
         ACCOUNT: "&ACCOUNT;"
         RES: "&RES;"
         TSTUB: "&TSTUB;"
-        WE2E_VAV: "&WE2E_VAV;"
+        WE2E_TEST: "&WE2E_TEST;"
         WE2E_ATOL: "&WE2E_ATOL;"
         WE2E_LOG_FN: "&WE2E_LOG_FN;"
         LOGDIR: "&LOGDIR;"
@@ -208,7 +208,7 @@ workflow:
         ATMOS_FORC: "&ATMOS_FORC;"
         RES: "&RES;"
         WARMSTART_DIR: "&WARMSTART_DIR;"
-        WE2E_VAV: "&WE2E_VAV;"
+        WE2E_TEST: "&WE2E_TEST;"
         WE2E_ATOL: "&WE2E_ATOL;"
         WE2E_LOG_FN: "&WE2E_LOG_FN;"
         LOGDIR: "&LOGDIR;"
@@ -259,6 +259,9 @@ workflow:
         KEEPDATA: "&KEEPDATA;"
         PDY: "&PDY;"
         cyc: "&cyc;"
+        LOGDIR: "&LOGDIR;"
+        NPROCS_ANALYSIS: "&NPROCS_ANALYSIS;"
+        NPROCS_FORECAST: "&NPROCS_FORECAST;"
       account: "&ACCOUNT;"
       command: '&HOMElandda;/parm/task_load_modules_run_jjob.sh "plot_stats" "&HOMElandda;" "&MACHINE;"'
       jobname: plot_stats
diff --git a/parm/land_analysis_orion.yaml b/parm/land_analysis_orion.yaml
index 1b8b1f61..05d60a1f 100644
--- a/parm/land_analysis_orion.yaml
+++ b/parm/land_analysis_orion.yaml
@@ -34,7 +34,7 @@ workflow:
     OBS_TYPES: "GHCN"
     DAtype: "letkfoi_snow"
     TSTUB: "oro_C96.mx100"
-    WE2E_VAV: "YES"
+    WE2E_TEST: "NO"
     WE2E_ATOL: "1e-7"
     WE2E_LOG_FN: "we2e.log"
     NET: "landda"
@@ -137,7 +137,7 @@ workflow:
         ACCOUNT: "&ACCOUNT;"
         RES: "&RES;"
         TSTUB: "&TSTUB;"
-        WE2E_VAV: "&WE2E_VAV;"
+        WE2E_TEST: "&WE2E_TEST;"
         WE2E_ATOL: "&WE2E_ATOL;"
         WE2E_LOG_FN: "&WE2E_LOG_FN;"
         LOGDIR: "&LOGDIR;"
@@ -172,7 +172,7 @@ workflow:
         ACCOUNT: "&ACCOUNT;"
         RES: "&RES;"
         TSTUB: "&TSTUB;"
-        WE2E_VAV: "&WE2E_VAV;"
+        WE2E_TEST: "&WE2E_TEST;"
         WE2E_ATOL: "&WE2E_ATOL;"
         WE2E_LOG_FN: "&WE2E_LOG_FN;"
         LOGDIR: "&LOGDIR;"
@@ -208,7 +208,7 @@ workflow:
         ATMOS_FORC: "&ATMOS_FORC;"
         RES: "&RES;"
         WARMSTART_DIR: "&WARMSTART_DIR;"
-        WE2E_VAV: "&WE2E_VAV;"
+        WE2E_TEST: "&WE2E_TEST;"
         WE2E_ATOL: "&WE2E_ATOL;"
         WE2E_LOG_FN: "&WE2E_LOG_FN;"
         LOGDIR: "&LOGDIR;"
@@ -259,6 +259,9 @@ workflow:
         KEEPDATA: "&KEEPDATA;"
         PDY: "&PDY;"
         cyc: "&cyc;"
+        LOGDIR: "&LOGDIR;"
+        NPROCS_ANALYSIS: "&NPROCS_ANALYSIS;"
+        NPROCS_FORECAST: "&NPROCS_FORECAST;"
       account: "&ACCOUNT;"
       command: '&HOMElandda;/parm/task_load_modules_run_jjob.sh "plot_stats" "&HOMElandda;" "&MACHINE;"'
       jobname: plot_stats
diff --git a/parm/templates/template.plot_hofx.yaml b/parm/templates/template.plot_hofx.yaml
deleted file mode 100644
index fbf9b0ee..00000000
--- a/parm/templates/template.plot_hofx.yaml
+++ /dev/null
@@ -1,8 +0,0 @@
-hofx_files: 'INPUTFP'
-field_var: 'XXFIELDVAR'
-field_range: [XXFRLOW,XXFRHIGH]
-nbins: XXNBINS
-plottype: 'XXPLOTTYPE'
-title_fig: 'XXFIGTITLE'
-output_prefix: 'XXPREOUTFN'
-machine: 'XXMACHINE'
diff --git a/parm/templates/template.plot_restart.yaml b/parm/templates/template.plot_restart.yaml
deleted file mode 100644
index d1def18e..00000000
--- a/parm/templates/template.plot_restart.yaml
+++ /dev/null
@@ -1,8 +0,0 @@
-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_analysis.sh b/scripts/exlandda_analysis.sh
index 6f76cb98..6967ee7f 100755
--- a/scripts/exlandda_analysis.sh
+++ b/scripts/exlandda_analysis.sh
@@ -251,8 +251,8 @@ if [[ -d output/DA/hofx ]]; then
   ln -nsf ${COMOUThofx}/* ${DATA_HOFX}
 fi
 
-# WE2E V&V
-if [[ "${WE2E_VAV}" == "YES" ]]; then
+# WE2E test
+if [[ "${WE2E_TEST}" == "YES" ]]; then
   path_fbase="${FIXlandda}/test_base/we2e_com/${RUN}.${PDY}"
   fn_sfc="${FILEDATE}.sfc_data.tile"
   fn_inc="${FILEDATE}.xainc.sfc_data.tile"
diff --git a/scripts/exlandda_forecast.sh b/scripts/exlandda_forecast.sh
index a660d935..78bf4f98 100755
--- a/scripts/exlandda_forecast.sh
+++ b/scripts/exlandda_forecast.sh
@@ -155,8 +155,8 @@ do
   ln -nsf ${COMOUT}/ufs_land_restart.${nYYYY}-${nMM}-${nDD}_${nHH}-00-00.tile${itile}.nc ${DATA_RESTART}
 done
 
-# WE2E V&V
-if [[ "${WE2E_VAV}" == "YES" ]]; then
+# WE2E test
+if [[ "${WE2E_TEST}" == "YES" ]]; then
   path_fbase="${FIXlandda}/test_base/we2e_com/${RUN}.${PDY}"
   fn_res="ufs_land_restart.${nYYYY}-${nMM}-${nDD}_${nHH}-00-00.tile"
   we2e_log_fp="${LOGDIR}/${WE2E_LOG_FN}"
diff --git a/scripts/exlandda_plot_stats.sh b/scripts/exlandda_plot_stats.sh
index a8203fbc..8a026f27 100755
--- a/scripts/exlandda_plot_stats.sh
+++ b/scripts/exlandda_plot_stats.sh
@@ -7,12 +7,15 @@ MM=${PDY:4:2}
 DD=${PDY:6:2}
 HH=${cyc}
 
+nYYYY=${NTIME:0:4}
+nMM=${NTIME:4:2}
+nDD=${NTIME:6:2}
+nHH=${NTIME:8:2}
+
 ############################################################
 # Stats Plot
 ############################################################
 
-cp -p ${PARMlandda}/templates/template.plot_hofx.yaml plot_hofx.yaml
-
 # Path to the directory containing the input file
 INPUTFP="${DATA_HOFX}"
 # Field variable
@@ -29,15 +32,16 @@ FIGTITLE="GHCN Snow Depth (mm)::Obs-Ana::${PDY}"
 # Prefix of output file name
 PREOUTFN="hofx_oma_${PDY}"
 
-sed -i "s|INPUTFP|${INPUTFP}|g" plot_hofx.yaml
-sed -i -e "s/XXFIELDVAR/${FIELDVAR}/g" plot_hofx.yaml
-sed -i -e "s/XXFRLOW/${FRLOW}/g" plot_hofx.yaml
-sed -i -e "s/XXFRHIGH/${FRHIGH}/g" plot_hofx.yaml
-sed -i -e "s/XXNBINS/${NBINS}/g" plot_hofx.yaml
-sed -i -e "s/XXPLOTTYPE/${PLOTTYPE}/g" plot_hofx.yaml
-sed -i -e "s/XXFIGTITLE/${FIGTITLE}/g" plot_hofx.yaml
-sed -i -e "s/XXPREOUTFN/${PREOUTFN}/g" plot_hofx.yaml
-sed -i -e "s/XXMACHINE/${MACHINE}/g" plot_hofx.yaml
+cat > plot_hofx.yaml <<EOF
+hofx_files: '${INPUTFP}'
+field_var: '${FIELDVAR}'
+field_range: [${FRLOW},${FRHIGH}]
+nbins: ${NBINS}
+plottype: '${PLOTTYPE}'
+title_fig: '${FIGTITLE}'
+output_prefix: '${PREOUTFN}'
+machine: '${MACHINE}'
+EOF
 
 ${USHlandda}/hofx_analysis_stats.py
 if [[ $? != 0 ]]; then
@@ -50,32 +54,64 @@ cp -p ${PREOUTFN}* ${COMOUTplot}
 
 
 ############################################################
-# Restart Plot
+# Time-history Plot
 ############################################################
 
-nYYYY=${NTIME:0:4}
-nMM=${NTIME:4:2}
-nDD=${NTIME:6:2}
-nHH=${NTIME:8:2}
+FN_DATA_ANAL_PREFIX="analysis_"
+FN_DATA_ANAL_SUFFIX=".log"
+FN_DATA_FCST_PREFIX="forecast_"
+FN_DATA_FCST_SUFFIX=".log"
+OUT_TITLE_ANAL_BASE="Land-DA::Analysis::QC SnowDepthGHCN::"
+OUT_FN_ANAL_BASE="landda_timehistory_"
+OUT_TITLE_TIME="Land-DA::Wall-clock time"
+OUT_FN_TIME="landda_timehistory_wtime"
+
+cat > plot_timehistory.yaml <<EOF
+path_data: '${LOGDIR}'
+work_dir: '${DATA}'
+fn_data_anal_prefix: '${FN_DATA_ANAL_PREFIX}'
+fn_data_anal_suffix: '${FN_DATA_ANAL_SUFFIX}'
+fn_data_fcst_prefix: '${FN_DATA_FCST_PREFIX}'
+fn_data_fcst_suffix: '${FN_DATA_FCST_SUFFIX}'
+nprocs_anal: '${NPROCS_ANALYSIS}'
+nprocs_fcst: '${NPROCS_FORECAST}'
+out_title_anal_base: '${OUT_TITLE_ANAL_BASE}'
+out_fn_anal_base: '${OUT_FN_ANAL_BASE}'
+out_title_time: '${OUT_TITLE_TIME}'
+out_fn_time: '${OUT_FN_TIME}'
+EOF
+
+${USHlandda}/plot_analysis_timehistory.py
+if [[ $? != 0 ]]; then
+  echo "FATAL ERROR: Time-history plots failed"
+  exit 44
+fi
+
+# Copy result files to COMOUT
+cp -p ${OUT_FN_ANAL_BASE}* ${COMOUTplot}
+cp -p ${OUT_FN_TIME}* ${COMOUTplot}
+
 
-cp -p ${PARMlandda}/templates/template.plot_restart.yaml plot_restart.yaml
+############################################################
+# Restart Plot
+############################################################
 
-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
+cat > plot_restart.yaml <<EOF
+path_data: '${COMIN}'
+work_dir: '${DATA}'
+fn_data_base: '${FN_DATA_BASE}'
+fn_data_ext: '${FN_DATA_EXT}'
+soil_lvl_number: '${SOIL_LEVEL_NUMBER}'
+out_title_base: '${OUT_TITLE_BASE}'
+out_fn_base: '${OUT_FN_BASE}'
+machine: '${MACHINE}'
+EOF
 
 ${USHlandda}/plot_forecast_restart.py
 if [[ $? != 0 ]]; then
diff --git a/scripts/exlandda_post_anal.sh b/scripts/exlandda_post_anal.sh
index 7cd62552..0b38e1bd 100755
--- a/scripts/exlandda_post_anal.sh
+++ b/scripts/exlandda_post_anal.sh
@@ -76,8 +76,8 @@ do
   cp -p ${DATA}/ufs_land_restart.${YYYY}-${MM}-${DD}_${HH}-00-00.tile${itile}.nc ${COMOUT}/ufs_land_restart.anal.${YYYY}-${MM}-${DD}_${HH}-00-00.tile${itile}.nc
 done
 
-# WE2E V&V
-if [[ "${WE2E_VAV}" == "YES" ]]; then
+# WE2E test
+if [[ "${WE2E_TEST}" == "YES" ]]; then
   path_fbase="${FIXlandda}/test_base/we2e_com/${RUN}.${PDY}"
   fn_res="ufs_land_restart.anal.${YYYY}-${MM}-${DD}_${HH}-00-00.tile"
   we2e_log_fp="${LOGDIR}/${WE2E_LOG_FN}"
diff --git a/ush/plot_analysis_timehistory.py b/ush/plot_analysis_timehistory.py
new file mode 100755
index 00000000..daaa4fda
--- /dev/null
+++ b/ush/plot_analysis_timehistory.py
@@ -0,0 +1,294 @@
+#!/usr/bin/env python3
+
+###################################################################### CHJ #####
+## Name		: plot_analysis_timehistory.py
+## Usage	: Plot timehistory of analysis output
+## Input files  : analysis_YYYYMMDDHH.log
+## NOAA/EPIC
+## History ===============================
+## V000: 2024/10/14: Chan-Hoo Jeon : Preliminary version
+## V001: 2024/10/15: Chan-Hoo Jeon : Add wall-clock time plot
+###################################################################### CHJ #####
+
+import os, sys
+import pathlib
+import yaml
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import matplotlib.colors as colors
+import matplotlib.ticker
+import matplotlib as mpl
+
+
+# Main part (will be called at the end) ============================= CHJ =====
+def main():
+# =================================================================== CHJ =====
+
+    yaml_file="plot_timehistory.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_anal_prefix = yaml_data['fn_data_anal_prefix']
+    fn_data_anal_suffix = yaml_data['fn_data_anal_suffix']
+    fn_data_fcst_prefix = yaml_data['fn_data_fcst_prefix']
+    fn_data_fcst_suffix = yaml_data['fn_data_fcst_suffix']
+    nprocs_anal = yaml_data['nprocs_anal']
+    nprocs_fcst = yaml_data['nprocs_fcst']
+    out_title_anal_base = yaml_data['out_title_anal_base']
+    out_fn_anal_base = yaml_data['out_fn_anal_base']
+    out_title_time = yaml_data['out_title_time']
+    out_fn_time = yaml_data['out_fn_time']
+
+    var_list = ["totalSnowDepth"]
+    nprocs_anal = int(nprocs_anal)
+    nprocs_fcst = int(nprocs_fcst)
+
+    # plot time-history
+    for var_nm in var_list:
+        var_dict_anal = get_data_analysis(path_data,fn_data_anal_prefix,fn_data_anal_suffix,nprocs_anal,var_nm)
+        var_dict_fcst = get_data_forecast(path_data,fn_data_fcst_prefix,fn_data_fcst_suffix,nprocs_fcst)
+        plot_data(var_dict_anal,var_dict_fcst,out_title_anal_base,out_fn_anal_base,out_title_time,out_fn_time,work_dir,var_nm)
+
+
+# Get data from files =============================================== CHJ =====
+def get_data_analysis(path_data,fn_data_anal_prefix,fn_data_anal_suffix,nprocs_anal,var_nm):
+# =================================================================== CHJ =====
+
+    print(' ===== var name: '+var_nm+' ========================')
+    # Find files with the sampe prefix
+    fp_data_anal_prefix = os.path.join(path_data,fn_data_anal_prefix)
+    files = []
+    for entry in os.scandir(path_data):
+        if entry.is_file() and entry.name.startswith(fn_data_anal_prefix):
+            files.append(entry.path)
+
+    files.sort()
+    #print("Files=",files)
+
+    nobs_qc_prefix = "QC SnowDepthGHCN totalSnowDepth"
+    wtime_oops_prefix = "OOPS_STATS util::Timers::Total"
+
+    file_date = []
+    min_val_final = []
+    max_val_final = []
+    rms_val_final = []
+    nobs_qc_final = []
+    nobs_in_final = []
+    wtime_oops = []
+    for file_fp in files:
+        file_date_raw = file_fp.removeprefix(fp_data_anal_prefix)
+        file_date_raw = file_date_raw.removesuffix(fn_data_anal_suffix)
+        file_date_tmp = f"{file_date_raw[0:4]}-{file_date_raw[4:6]}-{file_date_raw[6:8]}-{file_date_raw[8:10]}"
+        file_date.append(file_date_tmp)
+        print("File date=",file_date_tmp)
+
+        min_val_file = []
+        max_val_file = []
+        rms_val_file = []
+        nobs_qc_file = []
+        nobs_in_file = []
+        with open(file_fp, 'r') as file:
+            for line in file:
+                if line.startswith(var_nm):
+                    line_data_raw = line
+                    #print("Line data=",line_data_raw)
+                    line_split = line.split('| ')[1].split(' ')
+                    #print("Line split=",line_split)
+                    min_var = line_split[0].split(':')[0]
+                    min_val = line_split[0].split(':')[1]
+                    min_val = float(min_val)
+                    min_val_file.append(min_val)
+                    #print(min_var,"=",min_val,type(min_val))
+                    max_var = line_split[1].split(':')[0]
+                    max_val = line_split[1].split(':')[1]
+                    max_val = float(max_val)
+                    max_val_file.append(max_val)
+                    #print(max_var,"=",max_val,type(max_val))
+                    rms_var = line_split[2].split(':')[0]
+                    rms_val = line_split[2].split(':')[1]
+                    rms_val = float(rms_val)
+                    rms_val_file.append(rms_val)
+                    #print(rms_var,"=",rms_val,type(rms_val))
+
+                if line.startswith(nobs_qc_prefix):
+                    line_data_raw = line
+                    line_split = line.split(': ')[1].split(' ')
+                    #print("QC split=",line_split)
+                    if len(line_split) == 6 and line_split[4] != "of":
+                        nobs_qc_val = int(line_split[0])
+                        nobs_qc_file.append(nobs_qc_val)
+                        nobs_in_val = int(line_split[4])
+                        nobs_in_file.append(nobs_in_val)
+                        #print("NOBS ini=",nobs_in_file,", QC=",nobs_qc_file)
+
+                if line.startswith(wtime_oops_prefix):
+                    line_wtime_raw = line
+                    line_split = line.split(' : ')[1].split(' ')
+                    line_split = list(filter(None, line_split))
+                    if len(line_split) == 5:
+                        wtime_oops_file = float(line_split[2])
+                        #print("WTIME OOPS AVG=",wtime_oops_file)
+
+        min_val_final.append(min_val_file[-1])
+        max_val_final.append(max_val_file[-1])
+        rms_val_final.append(rms_val_file[-1])
+        nobs_qc_final.append(nobs_qc_file[-1])
+        nobs_in_final.append(nobs_in_file[-1])
+        wtime_oops.append(wtime_oops_file)
+    # ms to sec 
+    wtime_oops = [x * 0.001 for x in wtime_oops]
+    tcpu_oops = [x * nprocs_anal for x in wtime_oops]
+
+    # Create dictionary
+    var_dict_anal = {
+        "Date": file_date,
+        "Min": min_val_final,
+        "Max": max_val_final,
+        "RMS": rms_val_final,
+        "nobs_QC": nobs_qc_final,
+        "nobs_in": nobs_in_final,
+        "wtime_oops": wtime_oops,
+        "tcpu_oops": tcpu_oops
+    }
+    print("DICT=",var_dict_anal)
+
+    return var_dict_anal
+
+
+# Get data from files =============================================== CHJ =====
+def get_data_forecast(path_data,fn_data_fcst_prefix,fn_data_fcst_suffix,nprocs_fcst):
+# =================================================================== CHJ =====
+
+    # Find files with the sampe prefix
+    fp_data_fcst_prefix = os.path.join(path_data,fn_data_fcst_prefix)
+    files = []
+    for entry in os.scandir(path_data):
+        if entry.is_file() and entry.name.startswith(fn_data_fcst_prefix):
+            files.append(entry.path)
+
+    files.sort()
+    #print("Files=",files)
+
+    wtime_uwm_prefix = "The total amount of wall time"
+
+    file_date = []
+    wtime_uwm = []
+    for file_fp in files:
+        file_date_raw = file_fp.removeprefix(fp_data_fcst_prefix)
+        file_date_raw = file_date_raw.removesuffix(fn_data_fcst_suffix)
+        file_date_tmp = f"{file_date_raw[0:4]}-{file_date_raw[4:6]}-{file_date_raw[6:8]}-{file_date_raw[8:10]}"
+        file_date.append(file_date_tmp)
+        print("File date=",file_date_tmp)
+
+        with open(file_fp, 'r') as file:
+            for line in file:
+                if line.startswith(wtime_uwm_prefix):
+                    line_wtime_raw = line
+                    line_split = line.split(' = ')[1]
+                    wtime_uwm_file = float(line_split)
+                    #print("WTIME UFS Weather Model=",wtime_uwm_file)
+
+        wtime_uwm.append(wtime_uwm_file)
+
+    tcpu_uwm = [x * nprocs_fcst for x in wtime_uwm]
+
+    # Create dictionary
+    var_dict_fcst = {
+        "Date": file_date,
+        "wtime_uwm": wtime_uwm,
+        "tcpu_uwm": tcpu_uwm
+    }
+    print("DICT=",var_dict_fcst)
+
+    return var_dict_fcst
+
+
+# Plot data ========================================================= CHJ =====
+def plot_data(var_dict_anal,var_dict_fcst,out_title_anal_base,out_fn_anal_base,out_title_time,out_fn_time,work_dir,var_nm):
+# =================================================================== CHJ =====
+
+    out_title_anal = out_title_anal_base+var_nm
+    out_fn_anal = out_fn_anal_base+var_nm
+
+    dfa = pd.DataFrame(var_dict_anal)
+    dff = pd.DataFrame(var_dict_fcst)
+
+    txt_fnt=7
+    ln_wdth=0.75
+    mk_sz=3
+    
+    # PLOT 1
+    # figsize=(width,height) in inches
+    fig, axes = plt.subplots(nrows=3, ncols=1, sharex=True, figsize=(6,6))
+    fig.suptitle(out_title_anal,fontsize=txt_fnt+1,y=0.97)
+
+    axes[0].plot(dfa['Date'],dfa['Min'],'o-',color='blue',linewidth=ln_wdth,markersize=mk_sz,label='Min')
+    axes[0].plot(dfa['Date'],dfa['Max'],'s-.',color='red',mfc='none',linewidth=ln_wdth,markersize=mk_sz,label='Max')
+    axes[0].set_ylabel('Min / Max', fontsize=txt_fnt-1)
+    axes[0].tick_params(axis="y",labelsize=txt_fnt-2)
+    axes[0].legend(fontsize=txt_fnt-1, loc='center right')
+    axes[0].grid(linewidth=0.2)
+
+    axes[1].plot(dfa['Date'],dfa['RMS'],'o-',color='blue',linewidth=ln_wdth,markersize=mk_sz)
+    axes[1].set_ylabel('RMS', fontsize=txt_fnt-1)
+    axes[1].tick_params(axis="y",labelsize=txt_fnt-2)
+    axes[1].grid(linewidth=0.2)
+
+    axes[2].plot(dfa['Date'],dfa['nobs_in'],'o-',color='blue',linewidth=ln_wdth,markersize=mk_sz,label='N_obs:raw')
+    axes[2].plot(dfa['Date'],dfa['nobs_QC'],'s-.',color='red',mfc='none',linewidth=ln_wdth,markersize=mk_sz,label='N_obs:QC')
+    axes[2].set_xlabel('Date', fontsize=txt_fnt-1)
+    axes[2].set_ylabel('Number of observations', fontsize=txt_fnt-1)
+    axes[2].tick_params(axis="x",labelsize=txt_fnt-2)
+    axes[2].tick_params(axis="y",labelsize=txt_fnt-2)
+    axes[2].legend(fontsize=txt_fnt-1, loc='center right')
+    axes[2].grid(linewidth=0.2)
+
+    plt.xticks(rotation=30, ha='right')
+    plt.tight_layout()
+    # Output figure
+    ndpi = 300
+    out_file(work_dir,out_fn_anal,ndpi)
+
+    # PLOT 2
+    # figsize=(width,height) in inches
+    fig, axes = plt.subplots(nrows=2, ncols=1, sharex=True, figsize=(6,4))
+    fig.suptitle(out_title_time,fontsize=txt_fnt+1,y=0.95)
+
+    axes[0].plot(dfa['Date'],dfa['wtime_oops'],'o-',color='blue',linewidth=ln_wdth,markersize=mk_sz,label='Wall-clock')
+    axes[0].set_ylabel('Wall-clock time: OOPS (s)', fontsize=txt_fnt-1)
+    axes[0].tick_params(axis="y",labelsize=txt_fnt-2)
+    #axes[0].legend(fontsize=txt_fnt-1)
+    axes[0].grid(linewidth=0.2)
+
+    axes[1].plot(dff['Date'],dff['wtime_uwm'],'o-',color='blue',linewidth=ln_wdth,markersize=mk_sz,label='Wall-clock')
+    axes[1].set_xlabel('Date', fontsize=txt_fnt-1)
+    axes[1].set_ylabel('Wall-clock time: ufs_model (s)', fontsize=txt_fnt-1)
+    axes[1].tick_params(axis="x",labelsize=txt_fnt-2)
+    axes[1].tick_params(axis="y",labelsize=txt_fnt-2)
+    axes[1].grid(linewidth=0.2)
+
+    plt.xticks(rotation=30, ha='right')
+    plt.tight_layout()
+    # Output figure
+    ndpi = 300
+    out_file(work_dir,out_fn_time,ndpi)
+
+
+# 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()
+
diff --git a/versions/run.ver_hera b/versions/run.ver_hera
index 5d095668..7e5ae889 100644
--- a/versions/run.ver_hera
+++ b/versions/run.ver_hera
@@ -26,6 +26,7 @@ export py_cartopy_ver="0.21.1"
 export py_matplotlib_ver="3.7.3"
 export py_netcdf4_ver="1.5.8"
 export py_numpy_ver="1.22.3"
+export py_pandas_ver="1.5.3"
 export py_pyyaml_ver="6.0"
 export py_scipy_ver="1.11.3"
 export py_xarray_ver="2023.7.0"
diff --git a/versions/run.ver_hercules b/versions/run.ver_hercules
index 771fc8f2..72af1154 100644
--- a/versions/run.ver_hercules
+++ b/versions/run.ver_hercules
@@ -26,6 +26,7 @@ export py_cartopy_ver="0.21.1"
 export py_matplotlib_ver="3.7.3"
 export py_netcdf4_ver="1.5.8"
 export py_numpy_ver="1.22.3"
+export py_pandas_ver="1.5.3"
 export py_pyyaml_ver="6.0"
 export py_scipy_ver="1.11.3"
 export py_xarray_ver="2023.7.0"
diff --git a/versions/run.ver_orion b/versions/run.ver_orion
index a1b97600..a0ef6e99 100644
--- a/versions/run.ver_orion
+++ b/versions/run.ver_orion
@@ -26,6 +26,7 @@ export py_cartopy_ver="0.21.1"
 export py_matplotlib_ver="3.7.3"
 export py_netcdf4_ver="1.5.8"
 export py_numpy_ver="1.22.3"
+export py_pandas_ver="1.5.3"
 export py_pyyaml_ver="6.0"
 export py_scipy_ver="1.11.3"
 export py_xarray_ver="2023.7.0"