Boundary conditions - v13 based on MOM6/v2 ocean bathymetry (OM4)

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_cube.py

@@ -5,15 +5,37 @@
 import os
 from make_bcs_questionary import *
 from make_bcs_shared import * 
+from datetime import datetime
+import subprocess
 
 cube_template = """
 
 ln -s {MAKE_BCS_INPUT_DIR}/ocean/MOM5/360x200 data/MOM5/360x200
 ln -s {MAKE_BCS_INPUT_DIR}/ocean/MOM5/720x410 data/MOM5/720x410
 ln -s {MAKE_BCS_INPUT_DIR}/ocean/MOM5/1440x1080 data/MOM5/1440x1080
-ln -s {MAKE_BCS_INPUT_DIR}/ocean/MOM6/72x36 data/MOM6/72x36
-ln -s {MAKE_BCS_INPUT_DIR}/ocean/MOM6/540x458 data/MOM6/540x458
-ln -s {MAKE_BCS_INPUT_DIR}/ocean/MOM6/1440x1080 data/MOM6/1440x1080
+
+if ( {TRIPOL_OCEAN} == True ) then
+  set mom6v    = {mom6_bathy_version}
+  set mom6root = {MAKE_BCS_INPUT_DIR}/ocean/MOM6/$mom6v
+  set req      = {imo}x{jmo}
+
+  if ( ! -d $mom6root/$req ) then
+    echo "ERROR: MOM6/$mom6v/$req missing under {MAKE_BCS_INPUT_DIR}/ocean/MOM6"
+    echo "       Selected via questionnaire '{lbcsv}' -> MOM6_BATHY_VERSION=$mom6v"
+    exit 10
+  endif
+
+  if ( ! -d data/MOM6 ) mkdir -p data/MOM6
+  if ( -e data/MOM6/$req ) /bin/rm -f data/MOM6/$req
+  ln -s $mom6root/$req data/MOM6/$req
+endif
+
+ln -s {MAKE_BCS_INPUT_DIR}/ocean/MOM6/{mom6_bathy_version}/72x36 data/MOM6/72x36
+ln -s {MAKE_BCS_INPUT_DIR}/ocean/MOM6/{mom6_bathy_version}/540x458 data/MOM6/540x458
+ln -s {MAKE_BCS_INPUT_DIR}/ocean/MOM6/{mom6_bathy_version}/1440x1080 data/MOM6/1440x1080
+ln -s {MAKE_BCS_INPUT_DIR}/ocean/MOM6/{mom6_bathy_version}/720x576 data/MOM6/720x576
+ln -s {MAKE_BCS_INPUT_DIR}/ocean/MOM6/{mom6_bathy_version}/2880x2240 data/MOM6/2880x2240
+
 
 if( -e CF{NC}x6C{SGNAME}_{DATENAME}{IMO}x{POLENAME}{JMO}.stdout ) /bin/rm -f CF{NC}x6C{SGNAME}_{DATENAME}{IMO}x{POLENAME}{JMO}.stdout
 
@@ -137,6 +159,35 @@ def make_bcs_cube(config):
   if not os.path.exists(log_dir):
     os.makedirs(log_dir)
 
+  TOPO_VERSION       = topo_version_for_bcs(config['lbcsv'])
+  MOM6_BATHY_VERSION = mom6_bathy_version_for_bcs(config['lbcsv'])
+
+  # ---------- INPUT CHECKS (abort before sbatch) ----------
+  # 1) TOPO must exist for this grid
+  topo_dir = f"CF{NC}x6C{SGNAME}"              # CF0090x6C or CF0540x6C-SG001
+  topo_src = os.path.join(config['inputdir'], "atmosphere", "TOPO", TOPO_VERSION, topo_dir)
+  if not os.path.isdir(topo_src):
+      print(f"ABORT: Missing TOPO: {topo_src} "
+            f"(LBCSV={config['lbcsv']} TOPO_VERSION={TOPO_VERSION})")
+      return
+
+  # 2) MOM6 bathymetry: strict, check ONLY the size used by this run
+  if config["TRIPOL_OCEAN"]:
+      req = f"{config['imo']}x{config['jmo']}"  # 540x458 or 1440x1080
+      mom6_src = os.path.join(config['inputdir'], "ocean", "MOM6", MOM6_BATHY_VERSION, req)
+      if not os.path.isdir(mom6_src):
+          print(f"ABORT: Missing MOM6 bathymetry: {mom6_src} "
+                f"(LBCSV={config['lbcsv']} MOM6={MOM6_BATHY_VERSION})")
+          return
+
+  msg = f"[make_bcs_cube] LBCSV={config['lbcsv']}  TOPO={TOPO_VERSION}  GRID={GRIDNAME}"
+
+  if config["TRIPOL_OCEAN"]:
+      msg = msg + f"  MOM6={MOM6_BATHY_VERSION}  REQ={config['imo']}x{config['jmo']}"
+
+  print(msg)
+  # -------------------------------------------------------------------
+
   script_template = get_script_head() + cube_template + get_script_mv(config['grid_type'])
 
   script_string = script_template.format(\
@@ -175,6 +226,8 @@ def make_bcs_cube(config):
            STRETCH = STRETCH, \
            SGNAME  = SGNAME, \
            SGPARAM = SGPARAM, \
+           TOPO_VERSION = TOPO_VERSION,
+           mom6_bathy_version = MOM6_BATHY_VERSION,
            IS_STRETCHED = IS_STRETCHED, \
            NCPUS = config['NCPUS'])
 
@@ -211,6 +264,5 @@ def make_bcs_cube(config):
    answers = ask_questions()
    configs = get_configs_from_answers(answers)
    for config in configs:
-      if grid_type in ["Stretched_CS", "Cubed-Sphere"] :
-         make_bcs_cube(config)
-
+       if config['grid_type'] in ["Stretched_CS", "Cubed-Sphere"]:
+           make_bcs_cube(config)   

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py

@@ -65,7 +65,7 @@ def get_configs_from_answers(answers):
 
         maskfile = ''
 
-        if orslv in['O1','T2','T3','T4','T1MOM6','T3MOM6','T4MOM6']:
+        if orslv in['O1','T2','T3','T4','T1MOM6','T3MOM6','T4MOM6','T5MOM6','T8MOM6']:
            maskfile = 'GEOS5_10arcsec_mask_freshwater-lakes.nc'
            if lbcsv in ['F25', 'GM4', 'ICA']:
               maskfile = 'global.cat_id.catch.DL'
@@ -198,7 +198,7 @@ def ask_questions(default_grid="Cubed-Sphere"):
    "v10 : NL3 + PEATMAP + MODIS snow alb v2", \
    "v11 : NL3 + JPL veg height + PEATMAP + MODIS snow alb v2", \
    "v12 : NL3 + JPL veg height + PEATMAP + MODIS snow alb v2 + Argentina peatland fix", \
-   "v13 : NL3 + JPL veg height + PEATMAP + MODIS snow alb v2 + Argentina peatland fix + mean land elevation fix", \
+   "v13 : As in v13 + mean land elevation fix + MOM6 v2 (OM4) ocean-seaice bathymetry", \
    "ICA : Icarus        (archived*: /discover/nobackup/projects/gmao/bcs_shared/legacy_bcs/Icarus/)", \
    "GM4 : Ganymed-4_0   (archived*: /discover/nobackup/projects/gmao/bcs_shared/legacy_bcs/Ganymed-4_0/)", \
    "F25 : Fortuna-2_5   (archived*: n/a)"], 
@@ -320,6 +320,8 @@ def ask_questions(default_grid="Cubed-Sphere"):
                  "T1MOM6 --  Tripolar (MOM6-Tripolar-Ocean:    $72x36$  )", \
                  "T3MOM6 --  Tripolar (MOM6-Tripolar-Ocean:   $540x458$ )", \
                  "T4MOM6 --  Tripolar (MOM6-Tripolar-Ocean:  $1440x1080$)", \
+                 "T5MOM6 --  Tripolar (MOM6-Tripolar-Ocean:   $720x576$ )", \
+                 "T8MOM6 --  Tripolar (MOM6-Tripolar-Ocean:  $2880x2240$)", \
                  "CS     --  Cubed-Sphere Ocean  (Cubed-Sphere Data-Ocean)"],
             "when": lambda x:  "Stretched_CS" == x['grid_type'] or "Cubed-Sphere" == x['grid_type'],
         },

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_shared.py

@@ -4,6 +4,53 @@
 
 import os
 import glob
+import re
+
+# --- BEGIN VERSION MATRIX ---
+
+# Independent version mapping per questionnaire 'lbcsv'
+#
+# TOPO_VERSION refers to topography inputs for atm model
+
+_VERSION_MATRIX = {
+    "NL3": {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v1"},
+    "NL4": {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v1"},
+    "NL5": {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v1"},
+    "ICA": {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v1"},
+    "GM4": {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v1"},
+    "F25": {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v1"},
+    "v06": {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v1"},
+    "v07": {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v1"},
+    "v08": {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v1"},
+    "v09": {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v1"},
+    "v10": {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v1"},
+    "v11": {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v1"},
+    "v12": {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v1"},
+    "v13": {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v2"},
+}
+
+_DEFAULTS = {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v1"}
+
+def _normalize_lbcsv(label: str) -> str:
+    v = (label or "").strip()
+    m = re.match(r'[vV]?\d+|[A-Za-z0-9]+', v)  # accepts v13 / 13 / NL3 / etc.
+    if not m:
+        return v
+    key = m.group(0)
+    # normalize numeric like '13' -> 'v13'
+    if key.isdigit():
+        key = 'v' + key
+    return key
+
+def resolve_bcs_matrix(bcs_version: str):
+    key = _normalize_lbcsv(bcs_version)
+    return {**_DEFAULTS, **_VERSION_MATRIX.get(key, {})}
+
+def topo_version_for_bcs(bcs_version: str) -> str:
+    return resolve_bcs_matrix(bcs_version)["TOPO_VERSION"]
+
+def mom6_bathy_version_for_bcs(bcs_version: str) -> str:
+    return resolve_bcs_matrix(bcs_version)["MOM6_BATHY_VERSION"]
 
 def get_script_head() :
 
@@ -172,12 +219,36 @@ def get_script_mv(grid_type):
     echo "Successfully copied CO2_MonthlyMean_DiurnalCycle.nc4 to bcs dir."
 endif
 
-# adjust permissions
+"""
+
+   if grid_type in ("Cubed-Sphere", "Stretched_CS"):
+        mv_template = mv_template + """
 
-chmod +rX -R geometry land logs
+# Link (atm) TOPO into this BCS directory based on bcs_version (only needed if grid_type is [stretched] cube-sphere)
+set topo_version = {TOPO_VERSION}
+
+if ( ! -d TOPO ) mkdir -p TOPO
+set topo_dir  = CF{NC}x6C{SGNAME}     # e.g., CF0024x6C or CF0540x6C-SG001
+set topo_root = {MAKE_BCS_INPUT_DIR}/atmosphere/TOPO
+set topo_src  = $topo_root/$topo_version/$topo_dir
+
+if ( -e TOPO/$topo_dir ) then
+    echo "TOPO/$topo_dir already exists; not relinking."
+else if ( -d $topo_src ) then
+    /bin/ln -s $topo_src TOPO/$topo_dir
+    echo "Linked TOPO/$topo_dir -> $topo_src"
+else
+    echo "WARNING: TOPO source not found: $topo_src"
+endif
 
 """
 
+   mv_template = mv_template + """
+
+# adjust permissions (for all grid types)
+chmod +rX -R geometry land logs
+
+"""
    return mv_template
 
 def check_script(expdir, script):

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCatchParam.F90

@@ -12,7 +12,7 @@ PROGRAM mkCatchParam
 !     -y: Size of latitude dimension of input raster.             DEFAULT: 4320
 !     -b: Position of dateline w.r.t. first grid cell boundaries. DEFAULT: DC (dateline-on-center)
 !     -g: Gridname (name of the .til or .rst file without file extension)  
-!     -v: LBCSV : Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08, v09)
+!     -v: LBCSV : Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08, v09, v11, v12 and v13 )
 !     
 !
 ! This program is good to generate  

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rmTinyCatchParaMod.F90

@@ -56,7 +56,7 @@ module rmTinyCatchParaMod
   public Get_MidTime, Time_Interp_Fac
   public ascat_r0, jpl_canoph, NC_VarID, init_bcs_config  
 
-  ! The following variables define the details of the BCS version (data sources).
+  ! The following variables define the details of the BCs versions (data sources).
   ! Initialize to dummy values here and set to desired values in init_bcs_config().
 
   logical,      public, save :: use_PEATMAP = .false.
@@ -227,17 +227,17 @@ SUBROUTINE init_bcs_config(LBCSV)
        jpl_height  = .true.
 
     case ("v12","v13")  
-
-       ! "v12" and "v13" are identical except for:
-       ! - topography used for the atm (processed outside of make_bcs)
-       ! - bug fix for land elevation in catchment.def file
-       ! - generation of nc4-formatted tile file
+       !   v12 are v13 BCs are identical except for the following updates and bug fixes in v13:
+       !   - Fix for land elevation in catchment.def file.
+       !   - Generation of netCDF-4 (NC4) tile file (supplemental to ASCII file).
+       !   - Fix for inconsistency in land tile properties of coupled-model BCs vs. AGCM BCs.
+       !   - Coupled-model BCs use MOM6/v2 (OM4) ocean bathymetry.
 
        LAIBCS  = 'MODGEO'
        SOILBCS = 'HWSD_b'
        MODALB  = 'MODIS2'
        SNOWALB = 'MODC061v2'
-       OUTLETV = "v2"       
+       OUTLETV = "v2"
        GNU     = 1.0
        use_PEATMAP = .true.
        jpl_height  = .true.
@@ -5514,10 +5514,10 @@ SUBROUTINE svbksb(u,w,v,m,n,b,x)
     REAL*8 b(m),u(m,n),v(n,n),w(n),x(n) 
     PARAMETER (NMAX=500)  !Maximum anticipated value of n
     !------------------------------------------------------------------------------------------- 
-    ! Solves A · X = B for a vector X, where A is specified by the arrays u, w, v as returned by 
-    ! svdcmp. m and n are the dimensions of a, and will be equal for square matrices. b(1:m) is 
-    ! the input right-hand side. x(1:n) is the output solution vector. No input quantities are 
-    ! destroyed, so the routine may be called sequentially with different b’s. 
+    ! Solves A . X = B for a vector X, where A is specified by the arrays u, w, v as returned by 
+    ! svdcmp.  m and n are the dimensions of a, and will be equal for square matrices. b(1:m) is 
+    ! the input right-hand side. x(1:n) is the output solution vector.  No input quantities are 
+    ! destroyed, so the routine may be called sequentially with different b values. 
     !-------------------------------------------------------------------------------------------
 
     INTEGER i,j,jj 
@@ -5552,7 +5552,7 @@ SUBROUTINE svdcmp(a,m,n,w,v)
     PARAMETER (NMAX=500)  !Maximum anticipated value of n. 
     !-------------------------------------------------------------------------------------- 
     ! Given a matrix A(1:m,1:n), this routine computes its singular value decomposition, 
-    ! A = U · W · Vt. The matrix U replaces A on output. The diagonal matrix of singular 
+    ! A = U . W . Vt. The matrix U replaces A on output.  The diagonal matrix of singular 
     ! values W is output as a vector W(1:n). The matrix V (not the transpose Vt) is output 
     ! as V(1:n,1:n). 
     !--------------------------------------------------------------------------------------