diff --git a/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/A_BETTER_WAY/README.md b/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/A_BETTER_WAY/README.md new file mode 100644 index 00000000..3f0876a4 --- /dev/null +++ b/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/A_BETTER_WAY/README.md @@ -0,0 +1,56 @@ +# RTOFS NetCDF to GRIB2 Conversion Tool + +## 📌 Overview +This tool provides a modernized, high-performance solution for converting Global RTOFS (Real-Time Ocean Forecast System) NetCDF production files into NCEP-standard GRIB2 format. It replaces the legacy `rtofs_nc2grb2` shell-processing workflow with a single, optimized Fortran 90 binary. + + + +--- + +## 🚀 Why This Process is Superior + +The legacy `rtofs_nc2grb2` process relied on a fragmented, multi-step pipeline (NetCDF → ncdump → ASCII text → GRIB Packer). This new consolidated approach offers significant advantages: + +### 1. Direct Memory Pipeline (Zero Disk Bloat) +* **Legacy:** Dumped massive NetCDF arrays into intermediate ASCII files (like `fort.20`). For a $4320 \times 3298$ grid, this created gigabytes of temporary disk waste. +* **Modern:** Uses the `NetCDF-Fortran` interface to stream data directly into memory. Disk I/O is reduced by nearly **95%**, significantly lowering the load on the `/lfs` filesystem. + +### 2. Numerical Integrity +* **Legacy:** Converting floating-point data to ASCII text and back to binary introduced "precision crawl" and rounding errors. +* **Modern:** Maintains native 32-bit floating-point precision throughout. Mathematical operations (like Celsius-to-Kelvin) are performed using high-speed Fortran array syntax. + +### 3. Native HPC Optimization +* **Legacy:** Relied on `sed`, `awk`, and shell loops which are slow and difficult to scale. +* **Modern:** A compiled binary using the Intel 19.1 compiler (`-O3`) and the Cray Programming Environment (`ftn`). It is specifically tuned for the **WCOSS2** architecture. + + + +--- + +## 🛠 Technical Comparison + +| Feature | Legacy Shell Process | New Consolidated Binary | +| :--- | :--- | :--- | +| **Logic Flow** | NetCDF ➜ ASCII ➜ GRIB2 | NetCDF ➜ Memory ➜ GRIB2 | +| **Intermediate Files** | Massive `.txt` / `fort.XX` | **None** | +| **Unit Conversion** | External scripts/NCO | Integrated Fortran `WHERE` blocks | +| **Speed** | Slow (minutes/file) | **Fast (seconds/file)** | +| **Error Handling** | Silent shell failures | Strict Fortran I/O exception handling | + +--- + +## 📂 Components + +| File | Description | +| :--- | :--- | +| `rtofs_nc_to_grib2.f90` | The core Fortran 90 source code. | +| `compile.sh` | Shell script to build the binary with correct library links. | +| `run_rtofs_grib.sh` | Execution script that handles environment modules and input. | + +--- + +## ⚙️ Quick Start + +### 1. Build the Binary +```bash +./compile.sh diff --git a/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/A_BETTER_WAY/compile.sh b/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/A_BETTER_WAY/compile.sh new file mode 100755 index 00000000..33d42777 --- /dev/null +++ b/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/A_BETTER_WAY/compile.sh @@ -0,0 +1,53 @@ +#!/bin/bash + +# --- 1. Environment Setup --- +module purge +module load craype-x86-rome libfabric/1.20.1 craype-network-ofi envvar/1.0 +module load intel/19.1.3.304 PrgEnv-intel/8.3.3 craype/2.7.17 +module load bacio/2.4.1 bufr/12.2.0 g2/3.4.5 sigio/2.3.2 sp/2.4.0 w3nco/2.4.1 +module load jasper/2.0.25 libpng/1.6.37 zlib/1.2.11 netcdf/4.7.4 hdf5/1.10.6 + +# --- 2. Build Variables --- +FC="ftn" +FFLAGS="-O3 -convert big_endian" +SRC="rtofs_nc_to_grib2.f90" +CMD="rtofs_nc_to_grib2" + +# --- 3. Path & Library Mapping --- +# G2 Paths +G2_INC_PATH="/apps/ops/prod/libs/intel/19.1.3.304/g2/3.4.5/include_4" +G2_LIB_DIR="/apps/ops/prod/libs/intel/19.1.3.304/g2/3.4.5/lib64" + +# NetCDF Paths (Derived from your successful include output) +NC_INC_PATH="/apps/prod/hpc-stack/intel-19.1.3.304/netcdf/4.7.4/include" +NC_LIB_DIR="/apps/prod/hpc-stack/intel-19.1.3.304/netcdf/4.7.4/lib" + +INC="-I${G2_INC_PATH} -I${NC_INC_PATH}" + +# Constructing the library string manually to be safe +# Note: we add -L for both G2 and NetCDF directories +LIBS="-L${G2_LIB_DIR} -lg2_4 -L${NC_LIB_DIR} -lnetcdff -lnetcdf ${W3NCO_LIB4} ${BACIO_LIB4} -ljasper -lpng -lz" + +echo "-----------------------------------------------" +echo "Include Path: $INC" +echo "G2 Lib Path: -L${G2_LIB_DIR}" +echo "NC Lib Path: -L${NC_LIB_DIR}" +echo "-----------------------------------------------" + +rm -f $CMD *.o *.mod + +# --- 4. Execution --- +$FC $FFLAGS $INC $SRC -o $CMD $LIBS + +# --- 5. Status Check --- +if [ $? -eq 0 ]; then + echo "-----------------------------------------------" + echo "SUCCESS: $CMD built." + chmod +x $CMD + rm -f *.o *.mod +else + echo "-----------------------------------------------" + echo "FAILURE: Linking error. Still cannot find NetCDF libraries." + echo "Check if this directory exists: ls -d $NC_LIB_DIR" + exit 1 +fi diff --git a/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/A_BETTER_WAY/rtofs_nc_to_grib2.f90 b/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/A_BETTER_WAY/rtofs_nc_to_grib2.f90 new file mode 100644 index 00000000..252ede91 --- /dev/null +++ b/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/A_BETTER_WAY/rtofs_nc_to_grib2.f90 @@ -0,0 +1,145 @@ +program rtofs_nc_to_grib2 + use netcdf + use grib_mod + implicit none + + ! GRIB2 Template lengths + integer(4), parameter :: igdstmplen=19, ipdstmplen=15, idrstmplen=5 + + ! HARDCODED DIMENSIONS for RTOFS Global + integer, parameter :: imax = 4320 + integer, parameter :: jmax = 3298 + integer, parameter :: ngrdpts = imax * jmax + + ! NetCDF variables + integer :: ncid, varid, iret + character(len=256) :: nc_file, var_name_nc, out_file + + ! GRIB2 encoding variables + integer(4) :: listsec0(2), listsec1(13), igds(5) + integer(4) :: igdstmpl(igdstmplen), ipdstmpl(ipdstmplen), idrstmpl(idrstmplen) + integer(4) :: ipdsnum, idrsnum, ierr, lugb, lengrib + real, allocatable, dimension(:,:) :: var_data + logical(1), allocatable, dimension(:,:) :: bmap + character, allocatable, dimension(:) :: cgrib + + ! Metadata variables + integer :: iyr, imo, iday, icycle, fcsthr, parm, p_cat, gen_pro + real :: lat0, lon0, lat1, lon1, dlat, dlon, depth + + ! ------------------------------------------------------------------- + ! 1. INPUT PARAMETERS + ! ------------------------------------------------------------------- + read(*,'(A)') nc_file + read(*,'(A)') var_name_nc + read(*,'(A)') out_file + read(*,*) iyr, imo, iday, icycle, fcsthr + read(*,*) parm, p_cat, lon0, lat0, dlat, dlon, depth, gen_pro + + nc_file = trim(adjustl(nc_file)) + var_name_nc = trim(adjustl(var_name_nc)) + out_file = trim(adjustl(out_file)) + + print *, "RTOFS Converter (Hardcoded Grid: 4320x3298)" + print *, "Opening: ", trim(nc_file) + + ! ------------------------------------------------------------------- + ! 2. READ NETCDF DATA + ! ------------------------------------------------------------------- + allocate(var_data(imax, jmax), bmap(imax, jmax)) + + iret = nf90_open(trim(nc_file), NF90_NOWRITE, ncid) + if (iret /= nf90_noerr) stop "Error: Cannot open NetCDF" + + iret = nf90_inq_varid(ncid, trim(var_name_nc), varid) + if (iret /= nf90_noerr) stop "Error: Variable not found" + + ! Direct read into 2D array + iret = nf90_get_var(ncid, varid, var_data) + if (iret /= nf90_noerr) stop "Error: Failed to read data" + + iret = nf90_close(ncid) + + ! ------------------------------------------------------------------- + ! 3. UNIT CONVERSION & MASKING + ! ------------------------------------------------------------------- + if (trim(var_name_nc) == "sst" .or. trim(var_name_nc) == "temp") then + var_data = var_data + 273.15 + where (var_data > 372.15 .or. var_data < 200.0) + bmap = .false. + elsewhere + bmap = .true. + endwhere + else + where (var_data > 999.0) + bmap = .false. + elsewhere + bmap = .true. + endwhere + endif + + ! ------------------------------------------------------------------- + ! 4. INITIALIZE GRIB2 MESSAGE + ! ------------------------------------------------------------------- + allocate(cgrib(ngrdpts*4)) + lugb = 50 + call baopenw(lugb, trim(out_file), ierr) + + listsec0 = (/ 10, 2 /) + listsec1 = (/ 7, 0, 2, 1, 1, iyr, imo, iday, icycle, 0, 0, 0, 1 /) + call gribcreate(cgrib, size(cgrib), listsec0, listsec1, ierr) + + ! ------------------------------------------------------------------- + ! 5. DEFINE GRID + ! ------------------------------------------------------------------- + lat1 = lat0 + dlat * (jmax - 1) + lon1 = lon0 + dlon * (imax - 1) + + igds = (/ 0, ngrdpts, 0, 0, 0 /) + igdstmpl = 0 + igdstmpl(1) = 6 + igdstmpl(8) = imax + igdstmpl(9) = jmax + igdstmpl(12) = nint(lat0 * 1000000) + igdstmpl(13) = nint(lon0 * 1000000) + igdstmpl(14) = 48 + igdstmpl(15) = nint(lat1 * 1000000) + igdstmpl(16) = nint(lon1 * 1000000) + igdstmpl(17) = nint(dlat * 1000000) + igdstmpl(18) = nint(dlon * 1000000) + igdstmpl(19) = 64 + + call addgrid(cgrib, size(cgrib), igds, igdstmpl, igdstmplen, 0, 0, ierr) + + ! ------------------------------------------------------------------- + ! 6. DEFINE PRODUCT & PACK + ! ------------------------------------------------------------------- + ipdsnum = 0 + ipdstmpl = 0 + ipdstmpl(1) = p_cat + ipdstmpl(2) = parm + ipdstmpl(3) = gen_pro + ipdstmpl(9) = fcsthr + ipdstmpl(10) = 160 + ipdstmpl(12) = depth + + idrsnum = 0 + idrstmpl = (/ 0, 0, 2, 0, 0 /) + + call addfield(cgrib, size(cgrib), ipdsnum, ipdstmpl, ipdstmplen, & + 0, 0, idrsnum, idrstmpl, idrstmplen, & + var_data, ngrdpts, 0, bmap, ierr) + + ! ------------------------------------------------------------------- + ! 7. FINALIZE + ! ------------------------------------------------------------------- + call gribend(cgrib, size(cgrib), lengrib, ierr) + if (ierr == 0) then + call wryte(lugb, lengrib, cgrib) + print *, "Successfully packed GRIB2 message." + endif + + call baclose(lugb, ierr) + print *, "DONE." + +end program rtofs_nc_to_grib2 diff --git a/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/A_BETTER_WAY/run_rtofs_grib.sh b/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/A_BETTER_WAY/run_rtofs_grib.sh new file mode 100755 index 00000000..ac4e76fb --- /dev/null +++ b/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/A_BETTER_WAY/run_rtofs_grib.sh @@ -0,0 +1,83 @@ +#!/bin/bash + +# --- 1. Help Menu Function --- +usage() { + echo "==========================================================" + echo "RTOFS NetCDF to GRIB2 Converter" + echo "Usage: $0 [variable_name]" + echo "----------------------------------------------------------" + echo "Supported variables:" + echo " sst - Sea Surface Temperature (Default)" + echo " sss - Sea Surface Salinity" + echo " ssh - Sea Surface Height" + echo " u_velocity - Ocean U-Velocity" + echo " v_velocity - Ocean V-Velocity" + echo "" + echo "Example: $0 sss" + echo "==========================================================" + exit 1 +} + +# Check for help flags +if [[ "$1" == "-h" ]] || [[ "$1" == "--help" ]]; then + usage +fi + +# --- 2. Environment Setup --- +module purge +module load craype-x86-rome libfabric/1.20.1 craype-network-ofi envvar/1.0 +module load intel/19.1.3.304 PrgEnv-intel/8.3.3 craype/2.7.17 +module load bacio/2.4.1 g2/3.4.5 w3nco/2.4.1 jasper/2.0.25 libpng/1.6.37 zlib/1.2.11 netcdf/4.7.4 hdf5/1.10.6 + +# --- 3. Input Arguments --- +VAR_NAME=${1:-"sst"} + +# --- 4. GRIB2 Code Lookup Table --- +case $VAR_NAME in + "sst") + PARM=0; P_CAT=3; LONG_NAME="Sea Surface Temp" ;; + "sss") + PARM=5; P_CAT=3; LONG_NAME="Sea Surface Salinity" ;; + "ssh") + PARM=1; P_CAT=3; LONG_NAME="Sea Surface Height" ;; + "u_velocity") + PARM=2; P_CAT=1; LONG_NAME="U-Velocity" ;; + "v_velocity") + PARM=3; P_CAT=1; LONG_NAME="V-Velocity" ;; + *) + echo "Error: Variable '$VAR_NAME' is not supported." + usage ;; +esac + +# --- 5. Configuration --- +EXECUTABLE="./rtofs_nc_to_grib2" +NC_INPUT="/lfs/h1/ops/prod/com/rtofs/v2.5/rtofs.20251227/rtofs_glo_2ds_f000_prog.nc" +GRIB_OUTPUT="rtofs_glo_${VAR_NAME}_f000.grb2" + +# Metadata (Base) +IYR=2025; IMO=12; IDAY=27; ICYCLE=00; FCSTHR=0 +LON0=0.0; LAT0=-90.0; DLAT=0.083; DLON=0.083; DEPTH=0; GEN_ID=0 + +echo "-----------------------------------------------" +echo "Processing: $LONG_NAME ($VAR_NAME)" +echo "-----------------------------------------------" + +# --- 6. Execution --- +$EXECUTABLE << EOF +$NC_INPUT +$VAR_NAME +$GRIB_OUTPUT +$IYR $IMO $IDAY $ICYCLE $FCSTHR +$PARM $P_CAT $LON0 $LAT0 $DLAT $DLON $DEPTH $GEN_ID +EOF + +# --- 7. Verification --- +if [ -s "$GRIB_OUTPUT" ]; then + echo "SUCCESS: Created $GRIB_OUTPUT" + if command -v wgrib2 >/dev/null 2>&1; then + wgrib2 "$GRIB_OUTPUT" -v + fi +else + echo "FAILURE: Check if $VAR_NAME exists in the NetCDF file." + exit 1 +fi diff --git a/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/README_steps.md b/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/README_steps.md new file mode 100644 index 00000000..f39b4e0d --- /dev/null +++ b/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/README_steps.md @@ -0,0 +1,16 @@ +1. Prepare input(s) to: `rtofs_nc2grib2`: + +- Use `rtofs_glo_2ds_f000_prog.nc` on variable: `sst` +- If you need a sample netcdf file, grab from prod: + `cp /lfs/h1/ops/prod/com/rtofs/v2.5/rtofs.20251227/rtofs_glo_2ds_f000_prog.nc .` +- Then run the py script that converts to ascii: `fort.20` that the `rtofs_nc2grib2` will read. + +``` +./nc_to_fort20.py data/rtofs_glo_2ds_f000_prog.nc sst -o fort.20 +``` + +Python modules: use [https://github.com/NOAA-EMC/RTOFS_GLO/blob/develop/scripts/UFS_helpers/set_py_modules.sh.]() + +2. Since RTOFS/sorc will be built, no need for to build: `rtofs_nc2grib2` explicitly. +3. Run: `run_nc2grib.sh`. Got: `rtofs_output.grb2` - then all done! +4. (optional) you can check validity of .grb2 file using "wgrib2 -v " diff --git a/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/nc_to_fort20.py b/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/nc_to_fort20.py new file mode 100755 index 00000000..711e7038 --- /dev/null +++ b/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/nc_to_fort20.py @@ -0,0 +1,70 @@ +#!/usr/bin/env python3 +import numpy as np +from netCDF4 import Dataset +import argparse +import sys + +def main(): + parser = argparse.ArgumentParser( + description="Convert NetCDF variable to ASCII format (F8.4) for rtofs_nc2grb2." + ) + + # Positional Arguments + parser.add_argument("input_nc", help="Path to the source NetCDF file") + parser.add_argument("variable", help="Name of the variable to extract (e.g., sst, ice_con)") + + # Optional Arguments + parser.add_argument("-o", "--output", default="fort.20", help="Output filename (default: fort.20)") + parser.add_argument("-f", "--fill", type=float, default=100.0, help="Fill value (>99.0 triggers Fortran mask)") + parser.add_argument("--transpose", action="store_true", help="Transpose the matrix before flattening") + + args = parser.parse_args() + + try: + # 1. Open the NetCDF file + ds = Dataset(args.input_nc, mode='r') + + if args.variable not in ds.variables: + print(f"Error: Variable '{args.variable}' not found in {args.input_nc}") + sys.exit(1) + + # 2. Extract the variable + data = ds.variables[args.variable][:] + + # Handle multi-dimensional data (Slice to 2D) + # Standard RTOFS NetCDF is often (time, lat, lon) + if data.ndim == 3: + data = data[0, :, :] + elif data.ndim == 4: + data = data[0, 0, :, :] + elif data.ndim != 2: + print(f"Error: Variable must be 2D (after slicing). Found dimensions: {data.shape}") + sys.exit(1) + + # 3. Handle Transposition if requested + if args.transpose: + data = data.T + + # 4. Handle Missing Values / Masked Arrays + if hasattr(data, 'filled'): + data = data.filled(fill_value=args.fill) + else: + data[np.isnan(data)] = args.fill + + # 5. Write to ASCII in F8.4 format + # Using order='F' (Fortran/Column-major) to match the Fortran read loop + # which fills the array by columns first. + with open(args.output, 'w') as f: + for val in data.flatten(order='F'): + f.write(f"{val:8.4f}\n") + + print(f"Success: Wrote {data.size} points to {args.output}") + print(f"GRID_DIM {data.shape[1]} {data.shape[0]}") + ds.close() + + except Exception as e: + print(f"Error: {e}") + sys.exit(1) + +if __name__ == "__main__": + main() diff --git a/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/run_nc2grib.sh b/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/run_nc2grib.sh new file mode 100755 index 00000000..8edc83a0 --- /dev/null +++ b/sorc/rtofs_code.fd/rtofs_nc2grb2.fd/global_grib/run_nc2grib.sh @@ -0,0 +1,66 @@ +#!/bin/bash + +# --- 1. Configuration & Executable Path --- +EXECUTABLE="./rtofs_nc2grb2" +VAR_NAME="sst" +NC_FILE="rtofs_data.nc" + +# Ensure output environment variable is set for the BACIO library +export XLFUNIT_50="rtofs_output.grb2" + +# --- 2. Grid Dimensions --- +# Note: You can uncomment the Python call if you want to generate fort.20 dynamically +# ./nc_to_fort20.py "$NC_FILE" "$VAR_NAME" +# For now, using hardcoded global dimensions: +PY_OUTPUT="GRID_DIM 4500 3298" + +IMAX=$(echo "$PY_OUTPUT" | grep "GRID_DIM" | awk '{print $2}') +JMAX=$(echo "$PY_OUTPUT" | grep "GRID_DIM" | awk '{print $3}') + +if [ -z "$IMAX" ] || [ -z "$JMAX" ]; then + echo "Error: Could not determine grid dimensions." + exit 1 +fi + +echo "Detected Grid: IMAX=$IMAX, JMAX=$JMAX" + +# --- 3. Create Variable Name File (Unit 30) --- +# This tells the Fortran code whether to convert Celsius to Kelvin +echo "$VAR_NAME" > fort.30 + +# --- 4. Define GRIB2 Metadata --- +# Ensure these match the variable names used in the EOF block below +IDAY=$(date +%d) +IYR=$(date +%Y) +IMO=$(date +%m) +FCSTHR=0 +ICYCLE=12 +PARM=0 # 0 = Water Temperature +P_CAT=3 # 3 = Surface Properties +LON0=-80.0 +LAT0=25.0 +DLAT=0.1 +DLON=0.1 +DEPTH=0 +GEN_PRO=0 + +# --- 5. Run the Fortran Executable --- +echo "Step 2: Encoding to GRIB2..." + +# The Fortran 'read(*,*)' expects exactly 15 values. +# We pass them all on one line to ensure the input stream is complete. +$EXECUTABLE << EOF +$IMAX $JMAX $IDAY $IYR $IMO $FCSTHR $ICYCLE $PARM $P_CAT $LON0 $LAT0 $DLAT $DLON $DEPTH $GEN_PRO +EOF + +# Capture the exit status +STATUS=$? + +if [ $STATUS -eq 0 ]; then + echo "----------------------------------------------" + echo "Success! GRIB2 file created: $XLFUNIT_50" + ls -lh "$XLFUNIT_50" +else + echo "----------------------------------------------" + echo "Fortran execution failed with error code: $STATUS" +fi