diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt index 895b59f9d..bd8458931 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt @@ -37,7 +37,6 @@ else () target_compile_definitions(${this} PRIVATE STDC) endif () -esma_add_subdirectories (@ncar_topo utils_topo) ecbuild_add_executable (TARGET CombineRasters.x SOURCES CombineRasters.F90 LIBS MAPL ${this}) ecbuild_add_executable (TARGET mkCatchParam.x SOURCES mkCatchParam.F90 LIBS MAPL ${this} OpenMP::OpenMP_Fortran) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/scrip_to_restart_topo.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/scrip_to_restart_topo.py deleted file mode 100755 index 351a62039..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/scrip_to_restart_topo.py +++ /dev/null @@ -1,99 +0,0 @@ -#!/usr/bin/env python3 - -#------------- -# Load modules -#------------- -from netCDF4 import Dataset -import numpy -import argparse -import math - -def parse_args(): - p = argparse.ArgumentParser(description='convert old style cube to new style cube input') - p.add_argument('-i','--input',type=str,help='input file',default=None) - p.add_argument('-o','--output',type=str,help='output file',default=None) - p.add_argument('-g','--grid',type=str,help='grid identifier (e.g., sg001, sg002)',default=None) - return vars(p.parse_args()) - -#------------------ -# Opening the file -#------------------ -comm_args = parse_args() -Input_file = comm_args['input'] -Output_file = comm_args['output'] - -ncFid = Dataset(Input_file, mode='r') -ncFidOut = Dataset(Output_file, mode='w', format='NETCDF4') - -grid_type = comm_args['grid'] - -if grid_type == 'sg001': - ncFidOut.STRETCH_FACTOR = 2.5 - ncFidOut.TARGET_LAT = 39.5 - ncFidOut.TARGET_LON = -98.35 -elif grid_type == 'sg002': - ncFidOut.STRETCH_FACTOR = 3.0 - ncFidOut.TARGET_LAT = 39.5 - ncFidOut.TARGET_LON = -98.35 - -#--------------------- -# Extracting variables -#--------------------- - -ntiles = len(ncFid.dimensions['ncol']) -haveRdg = False -for dim in ncFid.dimensions: - if dim == 'nrdg': - haveRdg = True - rdgSize = len(ncFid.dimensions['nrdg']) - - -cRes = ntiles/6 -cRes = int(math.sqrt(cRes)) - -Xdim = ncFidOut.createDimension('lon',cRes) -Ydim = ncFidOut.createDimension('lat',cRes*6) - -if haveRdg: - rdgOut = ncFidOut.createDimension('unknown_dim1',rdgSize) - -vXdim = ncFidOut.createVariable('lon','f8',('lon')) -vYdim = ncFidOut.createVariable('lat','f8',('lat')) -setattr(ncFidOut.variables['lon'],'units','degrees_east') -setattr(ncFidOut.variables['lat'],'units','degrees_north') -setattr(ncFidOut.variables['lon'],'long_name','Longitude') -setattr(ncFidOut.variables['lat'],'long_name','Latitude') -vXdim[:]=range(1,cRes+1) -vYdim[:]=range(1,(6*cRes)+1) - -temp1d = numpy.zeros([6*cRes,cRes]) -if haveRdg: - temp2d = numpy.zeros([rdgSize,6*cRes,cRes]) - -exclude = ['lon','lat'] -for var in ncFid.variables: - if var not in exclude: - temp = ncFid.variables[var][:] - dim_size =len(temp.shape) - - if dim_size == 2: - tout = ncFidOut.createVariable(var,'f8',('unknown_dim1','lat','lon'),fill_value=1.0e15) - for att in ncFid.variables[var].ncattrs(): - if att != "_FillValue": - setattr(ncFidOut.variables[var],att,getattr(ncFid.variables[var],att)) - temp2d = numpy.reshape(temp,[rdgSize,cRes*6,cRes]) - tout[:,:,:] = temp2d[:,:,:] - - elif dim_size == 1: - tout = ncFidOut.createVariable(var,'f8',('lat','lon'),fill_value=1.0e15) - for att in ncFid.variables[var].ncattrs(): - if att != "_FillValue": - setattr(ncFidOut.variables[var],att,getattr(ncFid.variables[var],att)) - temp1d = numpy.reshape(temp,[cRes*6,cRes]) - tout[:,:] = temp1d[:,:] -#----------------- -# Closing the file -#---------------- -ncFidOut.close() -ncFid.close() - diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/CMakeLists.txt index d4eca3cfc..cc745169e 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/CMakeLists.txt @@ -1 +1 @@ -esma_add_subdirectories (soil routing) +esma_add_subdirectories (soil routing topography) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/.gitignore b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/.gitignore similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/.gitignore rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/.gitignore diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/CMakeLists.txt new file mode 100644 index 000000000..a8fe29cff --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/CMakeLists.txt @@ -0,0 +1,12 @@ +# CMakeLists.txt for Utils/Raster/preproc/topography +if(COMMAND esma_add_subdirectories) + esma_add_subdirectories(@ncar_topo utils_topo) +else() + if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/@ncar_topo/CMakeLists.txt") + add_subdirectory(@ncar_topo) + endif() + if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/utils_topo/CMakeLists.txt") + add_subdirectory(utils_topo) + endif() +endif() + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/README_topo.md b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/README_topo.md new file mode 100755 index 000000000..12d9008bf --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/README_topo.md @@ -0,0 +1,306 @@ +# Topography Pipeline — Quickstart & Reference + +This repository builds **GWD-ready topography** on cubed‑sphere grids (uniform and Schmidt‑stretched). The entry point for most users is **`make_topo.py`**, which interactively generates a Slurm job script and runs the full pipeline end‑to‑end via compiled utilities in `bin/`. + +--- + +## Quickstart + +```bash +# 1) Ensure binaries and modules are available +ls bin/{bin_to_cube.x,generate_scrip_cube_topo.x,cube_to_target.x,convert_to_gmao_output_topo.x} # should exist + +source g5_modules + +# 2) Run the driver +./make_topo.py + +# 3) Answer prompts (bin dir, output dir, input GMTED path, and resolutions) +# The script writes a Slurm job like topo_{res}.j in your chosen out_dir. + +# 4) Submit +cd +sbatch topo_c{res}.j +``` + +Outputs (per‑resolution) in `out_dir/output_/` as NetCDF files and binary alongside a GMAO restart. + +--- + +## What the pipeline does + +1. **Prepare a high‑res intermediate cube (`c3000`)** + + * Converts the global GMTED lat‑lon input + * Produces: `c3000.gmted_fixedanarticasuperior.nc` + +2. **Generate a SCRIP descriptor** (`generate_scrip_cube_topo.x`) + + * Uniform or Schmidt‑stretched (SG001/SG002) per `GenScrip.yaml` built by the job + * Writes `PEx-CF.nc4` + `c_coords.nc4` + +3. **Select `rrfac_max` (stretched grids only)** using `cdo` + `awk` + + * For **uniform grids**, `rrfac ≡ 1` everywhere, so `rrfac_max=1` and the flag is effectively a no‑op. + +4. **Remap to target grid** (`cube_to_target.x`) + + * Uses the intermediate cube + `--smoothing_scale` + optional `--rrfac_max` and Laplacian iterations + * Produces PE* NetCDF in `output_/` + +5. **Convert to restart + GMAO outputs** + + * `scrip_to_restart_topo.py` (adds `-g sg001/sg002` as needed) + * `convert_to_gmao_output_topo.x -i PE* --im ` + +--- + +## Inputs & dependencies + +* **Topography source (lat‑lon):** `gmted_fixed_anartica_superior_caspian.nc4` +* **Binaries (in `bin/`):** + + * `bin_to_cube.x` + * `generate_scrip_cube_topo.x` + * `cube_to_target.x` + * `convert_to_gmao_output_topo.x` + * `scrip_to_restart_topo.py` +* **Environment** + + * `bin/g5_modules` (loads platform compilers/MPI) + * Modules: `nco`, `cdo` , `python` + * `mpirun` available in `$PATH` +* **Scheduler:** Slurm (`sbatch`) + +The driver writes a **csh** job file + +--- + +## Running `make_topo.py` + +`make_topo.py` prompts for: + +* **bin_dir**: Absolute path to the `bin/` directory containing the compiled tools. +* **out_dir**: Where to put the generated Slurm script and outputs. Default: `/discover/nobackup//BCS_TOPO/`. +* **path_latlon**: Directory containing `gmted_fixed_anartica_superior_caspian.nc4`. +* **resolutions**: Pick from uniform (`C12 … C5760`) and stretched families (`SG001`, `SG002`). + + * If you choose a stretched family, you will be asked for its concrete IM sizes (e.g., `C270`, `C540`). + +The script writes **`topo_.j`** in `out_dir` and sets sane defaults for time, nodes, and job name. + +--- + +## Smoothing & GWD knobs + +* **`smoothing_scale`** (per grid) lives in the Python map inside `make_topo.py`. + + * *Bigger* `smoothing_scale` ⇒ *stronger* smoothing ⇒ *smaller* GWD (less roughness). +* **`ALPHA`** (stretched only) scales the Laplacian smooth step. Injected into `GenScrip.yaml` when applicable. +* **`RRFAC`** (regional refinement factor) **exists only on stretched grids**. Uniform grids have `rrfac=1` everywhere; you can ignore any `rrfac_max` logic for uniform runs. +* **Extra Laplacian cycles**: For `IM=5760`, the job adds `-l 13` to `cube_to_target.x`. + +### Defaults table + +| Grid | smoothing_scale | alpha (stretched only) | +| ------------- | --------------: | ---------------------: | +| C12 | 512.0 | – | +| C24 | 305.0 | – | +| C48 | 166.0 | – | +| C90 | 96.2 | – | +| C180 | 51.21 | – | +| C360 | 28.95 | – | +| C720 | 19.5 | – | +| C1120 | 8.26 | – | +| C1440 | 12.0 | – | +| C2880 | 3.285 | – | +| C5760 | 3.0 | – | +| C270 (SG001) | 100.0 | 2.9 | +| C540 (SG001) | 53.3 | 2.73 | +| C1080 (SG001) | 17.0 | 7.0 | +| C1536 (SG002) | 26.8 | 12.1 | +| C2160 (SG001) | 2.98 | 14.25 | + +To add/retune a grid, edit the `smoothmap`/`alpha` dicts in `make_topo.py` and re‑run. + +--- + +## Runtime guidance + +Wall estimates assume Discover nodes; adjust for your system. + +* `c180 … c2880, c1536, c1120, c2160`: **1 node**, ~**1h** +* `c5760, c540, c270, c48`: **2 nodes**, ~**4h** +* `c90`: **1 node**, ~**3h** +* `c24`: **2 nodes**, ~**8h** +* `c12`: **2 nodes**, ~**19h**, use `qos=long` + +The driver currently sets 1 node in the Slurm header by default; edit the generated job for heavy cases. + +--- + +## Outputs & where to find them + +For each `IM`: + +``` +/ + topo_.j # the job script you submit + output_/ + PEx* .nc # remapped topography (latest picked by timestamp) + gwd_internal_rst.* # restart written by scrip_to_restart_topo.py + *_gmao_*.nc # GMAO‑formatted outputs from convert_to_gmao_output_topo.x +``` + +The job also leaves the descriptor in `/PEx-CF.nc4`. + +--- + +## Stretching presets + +* **SG001** ⇒ `TARGET_LON=-98.35`, `TARGET_LAT=39.5`, `STRETCH_FACTOR=2.5` and `IM ∈ {270,540,1080,2160}` +* **SG002** ⇒ `TARGET_LON=-98.35`, `TARGET_LAT=39.5`, `STRETCH_FACTOR=3.0` and `IM = 1536` + +These are injected into `GenScrip.yaml` only for matching resolutions. +Users may add new stretching options with different target centers and stretch factors. Note, however, that the stretch factor must be less than 3.0 (3.0 is the maximum supported). + +--- + +## How to add a new resolution + +1. Edit `smoothmap` (and `alpha` if stretched) in `make_topo.py`. +2. If it’s a *stretched* grid, extend the SG001/SG002 lists in the job template section. +3. Re‑run `make_topo.py` and select the new resolution. + +--- + +## Design notes + +* The pipeline **builds `c3000` once** and **reuses** it across resolutions for consistency and speed. +* For **stretched** grids, `rrfac_max` is computed via `cdo infon` + `awk` on the SCRIP file to inform `cube_to_target.x`. +* For **uniform** grids, `rrfac=1`, so `rrfac_max` is 1 and can be omitted. +* For `IM=5760`, the driver adds `-l 13` Laplacian cycles for additional smoothing. + +--- + +## Example for one session + +``` +$ ./make_topo.py +Enter the root path of the bin. Where code was compiled install/bin. +Enter the path of the output directory: /discover/nobackup/{USER}/BCS_TOPO +Enter the path contains gmted_fixed_anartica_superior_caspian.nc4 (confirm selection): + /discover/nobackup/projects/gmao/bcs_shared/preprocessing_bcs_inputs/land/topo/v1/ +Select resolutions: [x] C360 [x] SG001 +Select resolution of SG001 grid: [x] C540 [ ] C1080 [ ] C2160 [ ] C270 +# => writes topo_c360_c540.j under out_dir +$ cd /discover/nobackup/{USER}/BCS_TOPO +$ sbatch topo_c360_c540.j +# or submit one job at the time to have clean run log for each job: topo_c360.j. +Save it is res dependent dir to have whole run isolated: /discover/nobackup/{USER}/BCS_TOPO/c360/ +``` + +--- + +## **List of Files in directory `utils_topo`** + +Below are the main files in this directory, with one-line summaries. +Each file also has its own inline header with more detail. + +- **CMakeLists.txt** — builds the utilities. +- **convert_bin_to_netcdf.F90** — small helper to convert raw binary topography into NetCDF (intermediate or testing). +- **convert_to_gmao_output.F90** — final step: produces GMAO deliverables (`gmted_DYN_ave_*.nc4`, `gmted_GWD_var_*.nc4`, `gmted_TRB_var_*.nc4`). +- **geompack.F90** — bundled geometry library (Burkardt routines: convex hull, triangle quality, etc.), needed for SCRIP generation. +- **generate_scrip_cube.F90** — builds SCRIP descriptors (uniform or Schmidt-stretched). +- **make_topo.py** — interactive driver, generates the Slurm job script to run the full pipeline. +- **scrip_to_cube_topo.py** — converts flat `ncol` → cube layout, using example file geometry. +- **scrip_to_restart_topo.py** — converts PE outputs into GWD restart format. + +--- + +## **List of Files in directory `@ncar_topo`** + +# **@ncar_topo — GMAO Version of NCAR Topography Tools** + +This directory contains the **GMAO-adapted version** of the NCAR *Topography Generation Software* +(originally `NCAR_Topo_2_0_1`, February 2023). +It provides the core Fortran utilities used to build cubed-sphere topography for the **GEOS** model. + +For the original NCAR user guide, see: +➡️ [https://github.com/NCAR/Topo/wiki/User's-Guide](https://github.com/NCAR/Topo/wiki/User's-Guide) + +### **`bin_to_cube/`** +Utility to **bin high-resolution global elevation data** (e.g., GMTED) from a latitude–longitude grid +onto a **cubed-sphere intermediate grid**. + +**Purpose:** +Creates an intermediate cube (e.g., `c3000.gmted_fixedanarticasuperior.nc`) used by the next stage, +`cube_to_target`. + +**Key files:** +- **bin_to_cube.F90** — main Fortran program +- **bin_to_cube\*.nl** — sample namelists for various cube sizes (e.g., 540, 3000) + +- Normally run **once** unless the source dataset or cube resolution changes. +- Generates data used for all downstream remapping and smoothing. + +--- + +### **`cube_to_target/`** +Utility to **process, smooth, and remap** the cubed-sphere intermediate topography to a **target model grid**, +either **uniform** or **Schmidt-stretched** variable resolution. + +**Capabilities:** +- Iterative **Laplacian “no-leak” smoother** +- Ridge detection and sub-grid statistics (`SGH`, `SGH30`, etc.) +- Seamless **variable-resolution** support using `rrfac` and `STRETCH_FACTOR` + +**Key files:** +- **cube_to_target.F90** — main driver +- **neighbor_search_mod.F90**, **kdtree_mod.F90** — geometric utilities +- **f90getopt.F90** — command line argument parser + +--- + +### **Upstream NCAR Software** + +Peter H. Lauritzen, Julio T. Bacmeister, Patrick Callaghan, and Mark A. Taylor (2015). +*NCAR Global Model Topography Generation Software for Unstructured Grids.* +*Geosci. Model Dev.*, **8**, 3975–3986. +DOI: **10.5194/gmd-8-3975-2015** + +--- + +### **GMAO Modifications** + +- Discover-compatible (Slurm) version +- Schmidt stretching and `rrfac` support +- Neighbor-repair logic using KD-tree search +- NetCDF compliance +- GEOS restart generation cleanup + +--- +## FAQ + +**Q: Can I reuse an existing `c3000`?** +A: Yes. The job checks and reuses `c3000.gmted_fixedanarticasuperior.nc` if present. + +**Q: Where do I change the stretch center/factor?** +A: Inside the job template section (look for `TARGET_LON/LAT` and `STRETCH_FACTOR`). + +**Q: Do uniform grids use RRFAC?** +A: No. Uniform grids have `rrfac=1` everywhere. The `rrfac_max` step is only meaningful for stretched grids. + +**Q: Can I tune GWD magnitude or general smoothing?** +A: Yes. smoothing_scale (all grids) and alpha (stretched only) control smoothing strength and GWD amplitude. Be cautious: aggressive changes can destabilize the Laplacian smoother. In such cases you may need to recode the smoother or carefully retune these knobs. + +**Q: Why does this code differ from the NCAR_Topo repository?** +A: There are several reasons: +This repository started from an outdated fork of NCAR_Topo. +We introduced major improvements in job submission and processing, especially through dynamic segments. +We cannot push changes back upstream (no test platform there), so their repo may follow different design choices. + +**Q: Can I widen or narrow the stretched‑grid refinement “dome”?** +A: Yes. The dome’s footprint is set by the half‑power radius (half_power_radius_deg) inside the stretched‑grid rrfac logic. By default it’s ~40°/sqrt(max(1, stretch_factor)), so for SF≈2.5 the half‑power radius is ~25° (covers most of CONUS). Increase the constant (e.g., 45–50) to broaden the dome; decrease to tighten it. This only affects stretched grids. + +*Last updated: 2025-10-08* diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/generate_topo.sh b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/obsolete/generate_topo.sh similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/generate_topo.sh rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/obsolete/generate_topo.sh diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/CMakeLists.txt similarity index 97% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/CMakeLists.txt rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/CMakeLists.txt index d7fe66616..828e2ed7d 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/CMakeLists.txt @@ -27,5 +27,4 @@ endif () install(PROGRAMS scrip_to_cube_topo.py DESTINATION bin) install(PROGRAMS scrip_to_restart_topo.py DESTINATION bin) -install(PROGRAMS generate_topo.sh DESTINATION bin) install(PROGRAMS make_topo.py DESTINATION bin) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/convert_bin_to_netcdf.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/convert_bin_to_netcdf.F90 similarity index 71% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/convert_bin_to_netcdf.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/convert_bin_to_netcdf.F90 index 86244ad16..197af5568 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/convert_bin_to_netcdf.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/convert_bin_to_netcdf.F90 @@ -1,4 +1,34 @@ program create_example +!------------------------------------------------------------------------------ +! Program: create_example +! +! Purpose: +! Convert a raw, unformatted 2D topography field (binary) into: +! (a) GEOS-style NetCDF with Xdim/Ydim and a 2D variable "z" (--output) +! (b) NCAR/SE-style NetCDF with 1D "ncol" and variable "PHIS" (--ncar) +! +! Inputs: +! -i, --input : raw unformatted binary (REAL*4/REAL*8) of size (IM, JM) +! --im : IM (x) size of the field +! --jm : JM (y) size; if omitted, assumes a cubed-sphere layout with JM=IM*6 +! +! Outputs (optional; supply flags to enable): +! -o, --output : write GEOS-style NetCDF to this path +! --ncar : write NCAR/SE-style NetCDF to this path +! +! Notes: +! * If --jm is not provided, JM defaults to IM*6 (cubed-sphere “isCube” mode). +! * GEOS output: +! - Dimensions: Xdim=IM, Ydim=JM +! - Variable: z(Xdim,Ydim) [units: m], with coordinate vars Xdim/Ydim (1..IM, 1..JM) +! * NCAR output: +! - Dimension: ncol = IM*JM (row-major flattening) +! - Variable: PHIS(ncol) [units: m] +! * Basic error handling via subroutine `check` (aborts with NetCDF error message). +! +!------------------------------------------------------------------------------ + + use netcdf use, intrinsic :: iso_fortran_env, only: REAL64 implicit none diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/convert_to_gmao_output.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/convert_to_gmao_output.F90 similarity index 83% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/convert_to_gmao_output.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/convert_to_gmao_output.F90 index 9df4a7292..e1bf251d9 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/convert_to_gmao_output.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/convert_to_gmao_output.F90 @@ -1,4 +1,39 @@ program create_example +!------------------------------------------------------------------------------ +! Program: convert_to_gmao_output (entry point: program create_example) +! +! Purpose: +! Convert NCAR-format PE topography outputs (with PHIS, SGH, SGH30) into +! GMAO-standard NetCDF and binary files, split into three deliverables: +! +! * gmted_DYN_ave_IMxJM.{nc4,data} → mean elevation (DYN), in meters +! * gmted_GWD_var_IMxJM.{nc4,data} → variance for Gravity Wave Drag (GWD), m^2 +! * gmted_TRB_var_IMxJM.{nc4,data} → variance for Turbulent Orography (TRB), m^2 +! +! Inputs: +! -i, --input : NCAR-style NetCDF containing +! * PHIS (m^2 s^-2) → converted to meters by dividing g +! * SGH (std. dev. for GWD, m) +! * SGH30 (std. dev. for TRB, m) +! --im : IM (x) size of the grid +! --jm : JM (y) size of the grid; if omitted, JM = IM*6 (cubed-sphere) +! +! Outputs: +! * NetCDF4 + raw unformatted binary pairs (.nc4 + .data) for each of DYN/GWD/TRB. +! * Dimensions: lon=IM, lat=JM. +! * Units: +! - DYN: meters (mean elevation) +! - GWD: m^2 (variance of subgrid orography for GWD) +! - TRB: m^2 (variance of subgrid orography for TRB) +! +! Notes: +! * Converts PHIS → elevation by dividing by g=9.80616. +! * Squares SGH/SGH30 (std. dev. in m) to produce variance in m^2. +! * Guards against missing/sentinel values (≈1e36). +! * Optionally smooths the poles in lat-lon mode (--jm given). +! +!------------------------------------------------------------------------------ + use netcdf implicit none diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/generate_scrip_cube.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/generate_scrip_cube.F90 similarity index 98% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/generate_scrip_cube.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/generate_scrip_cube.F90 index 881bc0a9e..a3eea5233 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/generate_scrip_cube.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/generate_scrip_cube.F90 @@ -2,6 +2,37 @@ #define _VERIFY(A) if(A/=0) call local_abort(A,__LINE__) program ESMF_GenerateCSGridDescription +!------------------------------------------------------------------------------ +! Program: ESMF_GenerateCSGridDescription (binary: generate_scrip_cube_topo.x) +! +! Purpose: +! Build a SCRIP descriptor and GEOS coordinate file for a cubed-sphere grid. +! Supports both uniform and Schmidt-stretched grids (target_lon/lat + factor). +! Writes per-cell centers/corners (deg), spherical area (sr), imask, and: +! - grid_fallback_mask : 1 if a tiny safe polygon was used for geometry +! - rrfac (stretched only) and rrfac_max (global attribute) +! +! Key features: +! • Robust corner ordering & periodicity handling (0/360 wraps). +! • Degeneracy repair & tiny-polygon fallback for bad cells. +! • Global min/max edge-length diagnostics and Gaussian rrfac for stretched. +! • MPI layout requires 6 PETs (one per cube face). +! +! Inputs: +! GenScrip.yaml with keys: +! CUBE_DIM, output_scrip, output_geos, +! [optional] DO_SCHMIDT, TARGET_LON, TARGET_LAT, STRETCH_FACTOR +! +! Outputs: +! - : SCRIP NetCDF (unstructured 1D layout) +! - : GEOS-style per-face centers/corners (deg) +! +! Notes: +! • Uniform grids do NOT carry rrfac (it is implicitly 1 everywhere). +! • For stretched grids, rrfac is radial (Gaussian) about the chosen midpoint. +! • Fallback cells are flagged; their areas are kept small but positive. +! +!------------------------------------------------------------------------------- ! ESMF Framework module use ESMF diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/geompack.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/geompack.F90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/geompack.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/geompack.F90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/make_topo.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/make_topo.py similarity index 92% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/make_topo.py rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/make_topo.py index ef9f29dd5..446216289 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/make_topo.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/make_topo.py @@ -8,6 +8,35 @@ # c90 1 node : 3h # c24 2 nodes : 8h # c12 2 nodes : 19h - has to be run with qos=long +#!/usr/bin/env python3 +""" +make_topo.py + +Purpose +------- +Interactive driver for the topography build pipeline. Prompts the user +for paths and resolutions, then writes a Slurm job script +(topo_.j) that runs the full sequence: + + 1. Build/reuse c3000 intermediate from GMTED lat-lon. + 2. Generate a SCRIP descriptor (uniform or Schmidt-stretched). + 3. Compute rrfac_max (stretched grids only). + 4. Run cube_to_target.x with smoothing/alpha. + 5. Convert outputs to GMAO restart and deliverables. + +Key knobs +--------- +- smoothing_scale : controls roughness → GWD amplitude (all grids). +- alpha : Laplacian scaling (stretched only). +- stretch center/factor : hard-coded for SG001/SG002, can be modified. +- runtime : varies by resolution; edit nodes/time in job if needed. + +Usage +----- + ./make_topo.py # interactively answer prompts + sbatch topo_.j # submit generated job script +""" + # import os import subprocess diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/scrip_to_cube_topo.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/scrip_to_cube_topo.py similarity index 83% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/scrip_to_cube_topo.py rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/scrip_to_cube_topo.py index 9608a1456..9332969af 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/scrip_to_cube_topo.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/scrip_to_cube_topo.py @@ -1,4 +1,32 @@ #!/usr/bin/env python3 +""" +scrip_to_cube_topo.py + +Purpose +------- +Convert old-style cubed-sphere data stored on a flat 'ncol' dimension +(ncol = 6 * IM * IM) into a GEOS-style cubed-sphere file with explicit +faces and (Xdim,Ydim), using an example file to copy the grid geometry. + +Inputs +------ +-i, --input : source NetCDF with 'ncol' (and optional 'nrdg') +-e, --example : GEOS-style example NetCDF providing 'lons', 'lats', + 'corner_lons', 'corner_lats', and Xdim/Ydim sizes +-o, --output : destination NetCDF4 +-v, --vars ... : optional list of variable names to convert (defaults to all) + +What it does +-------- +• Infers IM from the example file's Xdim; sets dims: + nf = 6, Xdim = IM, Ydim = IM, XCdim = IM+1, YCdim = IM+1 +• Reshapes variables: + (ncol) -> (nf, Ydim, Xdim) + (nrdg, ncol) -> (nrdg, nf, Ydim, Xdim) +• Copies variable attributes (excluding _FillValue). +• Writes grid-mapping metadata and coordinates ('lons','lats','corner_*') + by copying them from the example file. +""" #------------- # Load modules diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/scrip_to_restart_topo.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/scrip_to_restart_topo.py new file mode 100755 index 000000000..f45d379b3 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/scrip_to_restart_topo.py @@ -0,0 +1,140 @@ +#!/usr/bin/env python3 +""" +scrip_to_restart_topo.py + +Purpose: + Convert old-style cubed-sphere NetCDF files with a flat 'ncol' dimension + into new-style cube files with explicit (lat, lon) dimensions. + Optionally annotate stretched grids (sg001, sg002) with STRETCH_FACTOR + and target center information. + +Key points: + - Input: NetCDF with dimension 'ncol' (and possibly 'nrdg') + - Output: NetCDF with dimensions ('lat','lon'[,'unknown_dim1']) + - Reshapes variables from (ncol) → (lat, lon) or (nrdg,ncol) → (nrdg,lat,lon) + - Adds lon/lat index coordinate variables (1..IM, 1..6*IM) + - Sets global attrs for sg001/sg002 grids +""" +#------------- +# Load modules +#------------- +from netCDF4 import Dataset +import numpy +import argparse +import math + +def parse_args(): + p = argparse.ArgumentParser(description='convert old style cube to new style cube input') + p.add_argument('-i','--input',type=str,help='input file',default=None) + p.add_argument('-o','--output',type=str,help='output file',default=None) + p.add_argument('-g','--grid',type=str,help='grid identifier (e.g., sg001, sg002)',default=None) + return vars(p.parse_args()) + +#------------------ +# Opening the file +#------------------ +comm_args = parse_args() +Input_file = comm_args['input'] +Output_file = comm_args['output'] + +ncFid = Dataset(Input_file, mode='r') +ncFidOut = Dataset(Output_file, mode='w', format='NETCDF4') + +grid_type = comm_args['grid'] + +if grid_type == 'sg001': + ncFidOut.STRETCH_FACTOR = 2.5 + ncFidOut.TARGET_LAT = 39.5 + ncFidOut.TARGET_LON = -98.35 +elif grid_type == 'sg002': + ncFidOut.STRETCH_FACTOR = 3.0 + ncFidOut.TARGET_LAT = 39.5 + ncFidOut.TARGET_LON = -98.35 + +#--------------------- +# Extracting variables +#--------------------- + +ntiles = len(ncFid.dimensions['ncol']) +haveRdg = False +for dim in ncFid.dimensions: + if dim == 'nrdg': + haveRdg = True + rdgSize = len(ncFid.dimensions['nrdg']) + + +cRes = ntiles/6 +cRes = int(math.sqrt(cRes)) + +Xdim = ncFidOut.createDimension('lon',cRes) +Ydim = ncFidOut.createDimension('lat',cRes*6) + +if haveRdg: + rdgOut = ncFidOut.createDimension('unknown_dim1',rdgSize) + +vXdim = ncFidOut.createVariable('lon','f8',('lon')) +vYdim = ncFidOut.createVariable('lat','f8',('lat')) +setattr(ncFidOut.variables['lon'],'units','degrees_east') +setattr(ncFidOut.variables['lat'],'units','degrees_north') +setattr(ncFidOut.variables['lon'],'long_name','Longitude') +setattr(ncFidOut.variables['lat'],'long_name','Latitude') +vXdim[:]=range(1,cRes+1) +vYdim[:]=range(1,(6*cRes)+1) + +temp1d = numpy.zeros([6*cRes,cRes]) +if haveRdg: + temp2d = numpy.zeros([rdgSize,6*cRes,cRes]) + +exclude = ['lon', 'lat'] +for var in ncFid.variables: + if var in exclude: + continue + + v = ncFid.variables[var] + temp = v[:] + nd = temp.ndim + + # Choose fill value at creation time + is_angle = var.upper() in ('ANGLX', 'ANGLL') + fv = -9999.0 if is_angle else 1.0e15 + + # Create destination variable with correct fill immediately + if nd == 2: # (nrdg, ncol) -> (unknown_dim1, lat, lon) + tout = ncFidOut.createVariable( + var, 'f8', ('unknown_dim1', 'lat', 'lon'), + fill_value=fv + ) + elif nd == 1: # (ncol) -> (lat, lon) + tout = ncFidOut.createVariable( + var, 'f8', ('lat', 'lon'), + fill_value=fv + ) + else: + # unexpected rank — skip safely + continue + + # Copy attributes verbatim EXCEPT _FillValue (already set) + for att in v.ncattrs(): + if att != '_FillValue': + setattr(tout, att, getattr(v, att)) + + # For angle variables, make the metadata sentinel explicit + if is_angle: + setattr(tout, 'missing_value', -9999.0) + + # Simple reshape write (no masking/NaN munging in this tool) + if nd == 2: + tout[:, :, :] = numpy.reshape(temp, (len(ncFid.dimensions['nrdg']), + int((len(ncFid.dimensions['ncol'])//6)**0.5)*6, + int((len(ncFid.dimensions['ncol'])//6)**0.5))) + else: + tout[:, :] = numpy.reshape(temp, + (int((len(ncFid.dimensions['ncol'])//6)**0.5)*6, + int((len(ncFid.dimensions['ncol'])//6)**0.5))) + +#----------------- +# Closing the file +#---------------- +ncFidOut.close() +ncFid.close() +