diff --git a/.gitignore b/.gitignore index fe60ca11..fde5ace8 100644 --- a/.gitignore +++ b/.gitignore @@ -19,3 +19,6 @@ sorc/rtofs_hycom.fd/esmf_4_0_0rp2/obj sorc/rtofs_hycom.fd/esmf_4_0_0rp2/make.out sorc/*/*.log +# Ignore python cache +ush/__pycache__/* + diff --git a/scripts/exrtofs_glo_ncoda_inc.sh b/scripts/exrtofs_glo_ncoda_inc.sh index 70960c8c..5563900b 100755 --- a/scripts/exrtofs_glo_ncoda_inc.sh +++ b/scripts/exrtofs_glo_ncoda_inc.sh @@ -1,22 +1,5 @@ #!/bin/sh set -xa -############################################################################### -#### UNIX Script Documentation Block # -# # -# Script name: exrtofs_glo_ncoda_inc.sh # -# Script description: # -# # -# Author: Dan Iredell Org: NP23 Date: 2020-07-30 # -# # -# Abstract: Remap an archive file to an NCODA analysis, new layer depths. # -# # -# Sub-scripts called: # -# # -# Script history log: # -# 2020-07-30 Dan Iredell # -# 2023-02-08 Dmitry Dukhovskoy modified for updated ncoda_archv_lyrinc # -# # -############################################################################### export PS4='$SECONDS + ' @@ -70,7 +53,6 @@ types=salint_lyr_1o${SIZN} typeu=uucurr_lyr_1o${SIZN} typev=vvcurr_lyr_1o${SIZN} typep=lyrprs_lyr_1o${SIZN} -typec=icecov_sfc_1o${SIZN} ln -sf $COMINm1/rtofs_glo.t00z.n00.archv.a archv.${archday}.a ln -sf $COMINm1/rtofs_glo.t00z.n00.archv.b archv.${archday}.b @@ -118,25 +100,9 @@ else msg="$COMIN/ncoda/hycom_var/restart/${typep}_${dtg}_0000_analinc is missing" err_exit $msg fi -# Ice Coverage -if [ -e $COMIN/ncoda/hycom_var/restart/${typec}_${dtg}_0000_analfld ]; then - ln -sf $COMIN/ncoda/hycom_var/restart/${typec}_${dtg}_0000_analfld ./icecov_${dtg}_analfld -else - msg="$COMIN/ncoda/hycom_var/restart/${typec}_${dtg}_0000_analfld is missing" - err_exit $msg -fi -#create ssmi.r file -rm -f ssmi1.a ssmi2.a ssmi.$dtg.r -$EXECrtofs/rtofs_raw2hycom icecov_${dtg}_analfld $IDM $JDM 999.00 ssmi1.a > ssmi1.b -err=$?; export err ; err_chk -echo " error from rtofs_raw2hycom=",$err -$EXECrtofs/rtofs_hycom_expr ssmi1.a "ONE" $IDM $JDM 0.01 0 ssmi2.a > ssmi2.b -err=$?; export err ; err_chk -echo " error from rtofs_hycom_expr=",$err -$EXECrtofs/rtofs_hycom2raw8 ssmi2.a $IDM $JDM 1 1 $IDM $JDMA ssmi.$dtg.r -err=$?; export err ; err_chk -echo " error from rtofs_hycom2raw8=",$err +# Create sea ice concentration file: sic.nc +$USHrtofs/rtofs_glo_ice_update.sh ar=archv_1_inc.${archday} rm -f $ar.[a,b] diff --git a/scripts/load_py_modules.sh b/scripts/load_py_modules.sh new file mode 100755 index 00000000..18800ab0 --- /dev/null +++ b/scripts/load_py_modules.sh @@ -0,0 +1,7 @@ +#!/bin/bash + +export py_mod="ve/evs/2.0_py312" + +module use /apps/dev/modulefiles/ +module load ${py_mod} +module list diff --git a/ush/convert_bin_inc_to_nc.py b/ush/convert_bin_inc_to_nc.py new file mode 100755 index 00000000..47345ba9 --- /dev/null +++ b/ush/convert_bin_inc_to_nc.py @@ -0,0 +1,89 @@ +#!/usr/bin/env python3 + +""" +Convert analysis increments or (2d) fields from NCODA analysis +from binary to netcdf format. +""" + +import os +import sys +import argparse +import numpy as np +import xarray as xr + +from utils import bin_to_nc_2d + +# --------------------------------------------------------------------------- # +# 1. User Inputs +# --------------------------------------------------------------------------- # +get_inputs = argparse.ArgumentParser( + prog='convert_bin_inc_to_nc.py', + description='To convert RTOFS-DA output from binary to netcdf format.', + formatter_class=argparse.ArgumentDefaultsHelpFormatter +) + +get_inputs.add_argument('--topog_file', type=str, + help='Path to topography file (NetCDF)', + required=True) + +get_inputs.add_argument('--input_file', type=str, + help='Name of NCODA binary output, including path', + required=True) + +get_inputs.add_argument('--var_name', type=str, + help='2-d variable name: icecov, icethk, icetmp, or mixlyr', + default="icecov") + +get_inputs.add_argument('--output_path', type=str, + help='Path to where output file is to be saved', + required=True) + +get_inputs.add_argument('--plot', action='store_true', + help='If present, generate a diagnostic plot of the output file', + default=False) + +args = get_inputs.parse_args() + +# Assign variables +topo_file = args.topog_file +input_file = args.input_file +var_name = args.var_name +output_path = args.output_path + +# --------------------------------------------------------------------------- # +# 2. Safety Checks +# --------------------------------------------------------------------------- # + +if not os.path.isfile(topo_file): + sys.exit(f"ERROR: Topography file not found at: {topo_file}") + +if not os.path.isfile(input_file): + sys.exit(f"ERROR: NCODA binary file not found at: {input_file}") + +allowed_vars = ["icecov", "icethk", "icetmp", "mixlyr"] +if var_name not in allowed_vars: + # Exit code 2 to distinguish from path errors + print(f"ERROR: Variable '{var_name}' is not in allowed list: {allowed_vars}") + sys.exit(2) + +if not os.path.isdir(output_path): + sys.exit(f"ERROR: Output path does not exist: {output_path}") + +if not os.access(output_path, os.W_OK): + sys.exit(f"ERROR: Output path is not writable: {output_path}") + +print(f"STATUS: All inputs validated. Beginning conversion for {var_name}...") + +# --------------------------------------------------------------------------- # +# 3. Execution +# --------------------------------------------------------------------------- # +bin_to_nc_2d(input_file, var_name, output_path, topo_file) + +# --------------------------------------------------------------------------- # +# 4. Plotting (Optional) +# --------------------------------------------------------------------------- # +if args.plot: + from utils import plot_var + # Use the name generated by bin_to_nc_2d + out_nc_file = os.path.join(output_path, f"{os.path.basename(input_file)}.nc") + plot_var(out_nc_file, var_name, output_path) diff --git a/ush/rtofs_glo2d_ice.sh b/ush/rtofs_glo2d_ice.sh index c545bd45..b719e889 100755 --- a/ush/rtofs_glo2d_ice.sh +++ b/ush/rtofs_glo2d_ice.sh @@ -1,47 +1,103 @@ #!/bin/sh -# -############################################################################### -#### UNIX Script Documentation Block # -# # -# Script name: rtofs_glo2d_1hrly.sh # -# Script description: # -# # -# Authors: Bhavani Rajan & Ilya Rivin Org: NP23 Date: 2020-07-02 # -# # -# Abstract: This script creates the surface ice fields for RTOFS-Global # -# every 1 or 3 hours in netCDF format # -# # -# Sub-scripts called: # -# # -# Executables called: # -# rtofs_archv2ncdf2d # -# # -# # -# Imported variables: # -# RUN # -# modID # -# fhr # -# EXECrtofs # -# PARMglobal # -# mode # -# # -# # -# Script history log: # -# XXXX-XX-XXX Joe Dow # -# # -############################################################################### - -set -x -typeset -Z3 fhr - -echo "*** Started script $0 on hostname "$(hostname)' at time '$(date) - -export CDF021=${RUN}_${modID}_2ds_${mode}${fhr}_ice.nc - -export pgm="rtofs_archv2ncdf2d" -. prep_step -startmsg -${EXECrtofs}/rtofs_field2ncdf2d < ${PARMrtofs}/${RUN}_${modID}.${inputgrid}.field2ncdf2d.in >> $pgmout 2>>errfile -export err=$? ; err_chk - -echo "*** Finished script $0 on hostname "$(hostname)' at time '$(date) + +if [[ $# -ne 4 ]]; then + echo " " + echo "Usage: " + echo "$0" "topog_file ncoda_file variable_name output_path" + echo " " + echo "Example inputs: " + echo "topog_file: depth_GLBb0.08_09m11ob2_mom6.nc" + echo "ncoda_file: icecov_20251215_analfld" + echo "variable_name: options: icecov, icethk, icetmp, mixlyr" + echo "output_path: /lfs/h2/emc/stmp/santha.akella/test/" + echo " " + exit 1 +fi +echo " " + +set -ux + +msg="$(basename -- "$0") JOB has begun on $(hostname) at $(date)" +postmsg "$msg" + +# Allowed variable names +allowed_variables=("icecov" "icethk" "icetmp" "mixlyr") + +is_not_allowed=true +# --------------------------------------------------------------------------- # + +# Check for validity of inputs +topog_file=$1 +ncoda_file=$2 +var_name=$3 +output_path=$4 + +# Is there is a topography file in the "fix/", it is needed. +if [[ ! -f "${topog_file}" ]]; then + echo "A topography file is needed for this script to work." + echo "Check in ${FIXrtofs} for depth_*.nc" + exit 1 +fi + +# Ice coverage file +if [[ ! -f "${ncoda_file}" ]]; then + echo "A NCODA binary file is needed for this script to work." + echo "For example: icecov_yyyymmdd00_analfld" + exit 1 +fi + +# Variable name +if [[ -z "${var_name}" ]]; then + echo "Variable name (string) is empty. Valid options are: icecov, icethk, icetmp, mixlyr." + exit 1 +fi + +# Loop through the allowed variables to check for a match +for val in "${allowed_variables[@]}"; do + if [ "${var_name}" = "$val" ]; then + is_not_allowed=false + break # Exit loop early if a match is found + fi +done + +# Check the flag set within the loop +if [ "$is_not_allowed" = true ]; then + echo "Error: variable name '${var_name}' is not an allowed value." + echo "Valid options are: icecov, icethk, icetmp, mixlyr." + exit 2 +fi + +# Output path check +if [[ ! -d "${output_path}" ]]; then + echo "Error: Output path '${output_path}' does not exist or is not a directory." + exit 1 +fi + +if [[ ! -w "${output_path}" ]]; then + echo "Error: Output path '${output_path}' is not writable. Check permissions." + exit 1 +fi +# --------------------------------------------------------------------------- # + +# Load module(s) that provide python packages (such as xarray) +#source ${HOMErtofs}/scripts/load_py_modules.sh + +echo "Converting format for ${var_name}..." + +set +x # Turn off tracing to keep the log clean +$USHrtofs/convert_bin_inc_to_nc.py --topog_file "${topog_file}" \ + --input_file "${ncoda_file}" \ + --var_name "${var_name}" \ + --output_path "${output_path}" 2>&1 | tee output.log + +# Get the exit code of the Python script, NOT the tee command +err=${PIPESTATUS[0]} +export err +set -x # Turn tracing back on + +# Run error check +err_chk +echo " error from convert_bin_inc_to_nc.py = $err" + +msg="THE $(basename -- "$0") JOB HAS ENDED NORMALLY on $(hostname) at $(date)" +postmsg "$msg" diff --git a/ush/rtofs_glo_ice_update.sh b/ush/rtofs_glo_ice_update.sh new file mode 100755 index 00000000..efb92b5d --- /dev/null +++ b/ush/rtofs_glo_ice_update.sh @@ -0,0 +1,56 @@ +#!/bin/sh + +# +# Abstract: +# - Input: NCODA (data assimilation) output sea ice coverage (binary format). +# - Output: Converts to netcdf format (used by the sea ice model). +# +set -xa + +export PS4='$SECONDS + ' + +cd $DATA + +msg="$(basename -- "$0") JOB has begun on $(hostname) at $(date)" +postmsg "$msg" +# --------------------------------------------------------------------------- # + +# 1. Set up inputs for run + +dtg=${PDYm1}00 + +echo $dtg + +IDM=$(cat ${BLKDATA_FILE} | grep idm | cut -d' ' -f1 | tr -d '[:space:]') +JDM=$(cat ${BLKDATA_FILE} | grep jdm | cut -d' ' -f1 | tr -d '[:space:]') +SIZN="${IDM}x${JDM}" + +# 2. Link to ncoda var ice coverage restart file + +typec=icecov_sfc_1o${SIZN} + +if [ -e $COMIN/ncoda/hycom_var/restart/${typec}_${dtg}_0000_analfld ]; then + ln -sf $COMIN/ncoda/hycom_var/restart/${typec}_${dtg}_0000_analfld ./icecov_${dtg}_analfld +else + msg="$COMIN/ncoda/hycom_var/restart/${typec}_${dtg}_0000_analfld is missing" + err_exit $msg +fi + +# 3. Link in topo files + +ln -f -s ${FIXrtofs}/depth_GLBb0.08_09m11ob2_mom6.nc depth_GLBb0.08_09m11ob2_mom6.nc + +# 4. Create sic.nc file + +$USHrtofs/rtofs_glo2d_ice.sh depth_GLBb0.08_09m11ob2_mom6.nc icecov_${dtg}_analfld icecov ${DATA} +err=$?; export err ; err_chk +echo " error from rtofs_glo2d_ice.sh=",$err + +if [ -f "icecov_${dtg}_analfld.nc" ]; then + mv "icecov_${dtg}_analfld.nc" sic.nc +else + echo "WARNING: icecov_${dtg}_analfld.nc not found. Skipping rename." +fi + +msg="THE $(basename -- "$0") JOB HAS ENDED NORMALLY on $(hostname) at $(date)" +postmsg "$msg" diff --git a/ush/utils.py b/ush/utils.py new file mode 100644 index 00000000..1764feb5 --- /dev/null +++ b/ush/utils.py @@ -0,0 +1,155 @@ +#!/usr/bin/env python3 + +""" +Helper utilities to process NCODA output +""" + +import os +import sys +import xarray as xr +import numpy as np + +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt + +# --------------------------------------------------------------------------- # +# GLOBAL CONFIGURATION +# --------------------------------------------------------------------------- # +FILL_VALUE = -9.99e+33 +NCODA_SPVAL_THRESHOLD = -100.0 +PLOT_DPI = 60 # Low DPI for faster diagnostic previews + +def var_atts(vName): + """ + Returns (scale_factor, units, long_name, cmap) based on variable name. + """ + scale = 1.0 + units = "m" + long_name = vName + cmap = "viridis" + + if vName == "icecov": + scale = 0.01 + units = "fraction" + long_name = "Ice Concentration" + cmap = "Blues" + elif vName == "icethk": + units = "m" + long_name = "Ice Thickness" + cmap = "YlGnBu_r" + elif vName == "icetmp": + units = "degC" + long_name = "Ice Temperature" + cmap = "RdYlBu_r" + elif vName == "mixlyr": + units = "m" + long_name = "Mixed Layer Depth" + cmap = "plasma" + + return scale, units, long_name, cmap + +def read_ncoda_2d(inFile_path, ni, nj): + fName = os.path.basename(inFile_path) + parts = fName.split('_') + + if len(parts) < 2: + sys.exit(f"ERROR: Filename {fName} does not match expected format.") + + date_str = parts[1] + fDate = f"{date_str[0:4]}-{date_str[4:6]}-{date_str[6:8]}" + fTime = np.array([fDate], dtype='datetime64[D]') + + try: + with open(inFile_path, 'rb') as f: + vals = np.fromfile(f, dtype='>f4', count=ni*nj) + if vals.size != (ni * nj): + sys.exit(f"ERROR: File {fName} size mismatch. Expected {ni*nj}, got {vals.size}") + vals = vals.reshape((nj, ni)) + except Exception as e: + sys.exit(f"ERROR: Failed to read or reshape binary file: {e}") + + return fTime, vals + +def plot_var(nc_file, vName, outPath): + """ + Diagnostic plotting. Note: vName here is the output name (e.g., 'sic'). + """ + if not os.path.exists(nc_file): + return + + ds = xr.open_dataset(nc_file) + plot_data = ds[vName].squeeze() + plot_data = plot_data.where(plot_data != FILL_VALUE) + + # Map back to original name if needed for cmap lookup + lookup_name = "icecov" if vName == "sic" else vName + _, _, _, v_cmap = var_atts(lookup_name) + + plt.figure(figsize=(10, 6)) + plot_data.plot(cmap=v_cmap, robust=True) + #plt.title(f"Diagnostic Plot: {vName}\n{os.path.basename(nc_file)}") + + plot_name = nc_file.replace('.nc', '.png') + plt.savefig(plot_name, dpi=PLOT_DPI, bbox_inches='tight') + plt.close('all') + print(f"DIAGNOSTIC PLOT SAVED: {plot_name}") + +def bin_to_nc_2d(inFile, vName, outPath, topo_fName): + # 1. Grid & Mask + try: + with xr.open_dataset(topo_fName, decode_times=False) as ds_topo: + nj, ni = ds_topo.depth.shape + ocean_mask = ~np.isnan(ds_topo.depth.values) + except Exception as e: + sys.exit(f"ERROR: Problem reading topography file {topo_fName}: {e}") + + # 2. Read Binary + val_date, vals = read_ncoda_2d(inFile, ni, nj) + + # 3. Clean Special Values + vals = np.where(vals < NCODA_SPVAL_THRESHOLD, FILL_VALUE, vals) + + # 4. Scaling + scale_factor, v_units, v_long_name, _ = var_atts(vName) + vals = np.where(vals != FILL_VALUE, vals * scale_factor, FILL_VALUE) + + # 5. Masking + vals = np.where(ocean_mask, vals, FILL_VALUE) + + # 6. Construct Dataset + # Handle the "sic" renaming logic here + out_var_name = "sic" if vName == "icecov" else vName + + ds_out = xr.Dataset( + data_vars={out_var_name: (("ny", "nx"), vals.astype(np.float32))}, + coords={"time": (("time"), val_date)} + ) + + # 7. Metadata & Encoding + ds_out[out_var_name].encoding = {"_FillValue": FILL_VALUE} + ds_out[out_var_name].attrs = { + "units": v_units, + "long_name": v_long_name, + "scale_factor_applied": scale_factor + } + + # 8. Stats + valid_data = vals[vals != FILL_VALUE] + if valid_data.size > 0: + v_min, v_max = valid_data.min(), valid_data.max() + print(f"Summary for {out_var_name}: Min={v_min:.4f}, Max={v_max:.4f}") + if vName == "icecov" and (v_max > 1.01 or v_min < -0.01): + print(f"CRITICAL WARNING: {out_var_name} out of bounds! [Min: {v_min:.4f}, Max: {v_max:.4f}]") + else: + print(f"WARNING: No valid ocean data found for {out_var_name}") + + # 9. Save + out_name = os.path.join(outPath, f"{os.path.basename(inFile)}.nc") + try: + ds_out.to_netcdf(out_name) + print(f"SUCCESS: Created NetCDF [{out_var_name}] at {out_name}") + except Exception as e: + sys.exit(f"ERROR: Failed to write NetCDF {out_name}: {e}") + + return ds_out