Update streams and namelist options for ismip7_run

compass/landice/tests/ismip7_run/ismip7_ais/albany_input.yaml

@@ -2,23 +2,22 @@
 ---
 ANONYMOUS:
   Problem:
+    Depth Integrated Model: true
     LandIce Field Norm:
       sliding_velocity_basalside:
         Regularization Type: Given Value
         Regularization Value: 1.0e-4
     LandIce BCs:
       BC 0:
+        Cubrature Degree: 8
         Basal Friction Coefficient:
           Type: Power Law
-          Power Exponent: 0.2
+          Power Exponent: 0.33333333333
           Mu Type: Field
-          Effective Pressure Type: Constant
-          Effective Pressure: 1.0
-          # Zero Effective Pressure On Floating Ice At Nodes: true
+          Effective Pressure Type: Hydrostatic Computed At Nodes
+          Zero Effective Pressure On Floating Ice At Nodes: true
           Zero Beta On Floating Ice: true
 
-    Cubature Degree: 8
-
 # Discretization Description
   Discretization:
     #Exodus Output File Name: albany_output.exo

compass/landice/tests/ismip7_run/ismip7_ais/ismip7_ais.cfg

@@ -59,18 +59,6 @@ calving_method = restore
 # Only required if calving_method is set to 'von_mises'.
 von_mises_parameter_path = NotAvailable
 
-# Whether to apply the ISMIP7 Path C ice shelf collapse mask (from
-# ismip7_forcing fracture) as hydrofracture-gated mask calving. Expected
-# layout: {forcing_basepath}/{model}_{scenario}/shelf_collapse/*.nc
-# Not used for historical, ctrl, or ocx experiments. If true and the mask
-# file is missing for an experiment that should have one, setup fails.
-use_hydrofracture_forcing = true
-
-# Ice fracture toughness (K_IC, Pa m^0.5) used in the hydrofracture
-# vulnerability criterion gating mask calving. Only used when
-# use_hydrofracture_forcing is true.
-calving_fracture_toughness = 2.0e5
-
 # Which glacial isostatic adjustment (GIA) coupling to use. Options:
 #   none         - no GIA coupling (fixed bed)
 #   1dSLM        - couple to the 1-D sea-level model (uses the slm_* and

compass/landice/tests/ismip7_run/ismip7_ais/ismip7_ais_test.cfg

@@ -66,6 +66,14 @@ von_mises_parameter_path = NotAvailable
 #   FastIsostasy - couple to FastIsostasy (uses the fastiso_* options below)
 gia_model = FastIsostasy
 
+# Whether to use damage-based calving with a spatially-varying threshold
+use_damage_calving = false
+
+# Path to the file containing the damage calving threshold field.
+# If provided, the damageCalvingThreshold stream will be added.
+# Actual use is controlled by namelist options in MALI.
+damage_calving_threshold_path = NotAvailable
+
 # Path to the directory containing globally defined ice thickness
 # field for the sea-level model
 slm_input_ice = NotAvailable

compass/landice/tests/ismip7_run/ismip7_ais/namelist.landice

@@ -5,8 +5,10 @@
     config_nonconvergence_error = .false.
     config_flowParamA_calculation = 'PB1982'
 
-    config_thickness_advection = 'fo'
-    config_tracer_advection = 'fo'
+    config_thickness_advection = 'fct'
+    config_tracer_advection = 'fct'
+    config_horiz_tracer_adv_order = 3
+    config_advection_coef_3rd_order = 1.0
     config_zero_sfcMassBalApplied_over_bare_land = .true.
     config_apply_smb_sat_lapse_rate = .true.
 
@@ -16,14 +18,17 @@
     config_SLM_to_MALI_weights_file = 'mapfile_slm_to_mali.nc'
     config_fastisostasy_coupling_interval = 1
 
-    config_calving = 'none'
-    config_apply_calving_mask = .false.
+    config_calving = 'crevasse_depth'
+    config_apply_calving_mask = .true.
+    config_require_extensional_stresses_for_mask_calving = .true.
+    config_calving_fracture_toughness = 2.0e5
     config_restore_calving_front_prevent_retreat = .false.
-    config_calculate_damage = .true.
+    config_calculate_damage = .false.
     config_damage_calving_threshold = 0.95
+    config_damage_calving_threshold_source = 'scalar'
     config_damage_calving_method = 'none'
     config_calving_speed_limit = 0.00063492063
-    config_restore_calving_front = .true.
+    config_restore_calving_front = .false.
     config_remove_icebergs = .true.
     config_remove_small_islands = .true.
     config_distribute_unablatedVolumeDynCell = .true.
@@ -36,7 +41,7 @@
     config_surface_air_temperature_source = 'file'
     config_basal_heat_flux_source = 'file'
 
-    config_basal_mass_bal_float = 'ismip6'
+    config_basal_mass_bal_float = 'ismip7'
     config_front_mass_bal_grounded = 'ismip6'
     config_use_3d_thermal_forcing_for_face_melt = .true.
     config_add_ocean_thermal_forcing = 0.0
@@ -48,7 +53,9 @@
     config_dynamic_thickness = 10.0
 
     config_dt = '0000-01-00_00:00:00'
-    config_time_integration = 'forward_euler'
+    config_time_integration = 'runge_kutta'
+    config_rk_order = 3
+    config_rk3_stages = 3
     config_adaptive_timestep = .true.
     config_adaptive_timestep_calvingCFL_fraction = 0.8
     config_adaptive_timestep_include_calving = .true.

compass/landice/tests/ismip7_run/ismip7_ais/set_up_experiment.py

@@ -57,18 +57,24 @@ def setup(self):  # noqa: C901
         forcing_basepath = section.get('forcing_basepath')
         init_cond_path = section.get('init_cond_path')
         init_cond_fname = os.path.split(init_cond_path)[-1]
+        # Extract mesh name from init cond filename
+        # (e.g., AIS_2to20km_r03_20260922)
+        # Typical format: AIS_<resolution>_<version>_<date>_<something>.nc
+        # Strip .nc extension first, then take first 4 underscore-separated
+        # parts
+        base = init_cond_fname.replace('.nc', '')
+        mesh_name = '_'.join(base.split('_')[:4])
+        print(f"    Extracted mesh name: {mesh_name} from init cond: "
+              f"{init_cond_fname}")
         melt_params_path = section.get('melt_params_path')
         melt_params_fname = os.path.split(melt_params_path)[-1]
         region_mask_path = section.get('region_mask_path')
         region_mask_fname = os.path.split(region_mask_path)[-1]
         reference_surface_path = section.get('reference_surface_path')
         reference_surface_fname = os.path.split(reference_surface_path)[-1]
         calving_method = section.get('calving_method')
-        use_hydrofracture_forcing = section.getboolean(
-            'use_hydrofracture_forcing')
-        calving_fracture_toughness = section.get(
-            'calving_fracture_toughness')
         sea_level_model, fastisostasy = parse_gia_model(config)
+        sea_level_model = section.getboolean('sea_level_model')
 
         exp_info = self.exp_info
         scenario = exp_info['scenario']
@@ -122,27 +128,35 @@ def setup(self):  # noqa: C901
             # Find atmosphere climatology files
             smb_files = glob.glob(os.path.join(ctrl_atm_path, '*SMB*.nc'))
             smb_files = [f for f in smb_files if 'gradient' not in f]
+            # Filter by mesh name to handle directories with multiple meshes
+            smb_files = [f for f in smb_files if mesh_name in f]
             if len(smb_files) == 1:
                 smb_fname = os.path.split(smb_files[0])[-1]
                 os.symlink(smb_files[0],
                            os.path.join(self.work_dir, smb_fname))
             else:
-                sys.exit(f"ERROR: Expected 1 SMB climatology file in "
-                         f"{ctrl_atm_path}, found {len(smb_files)}")
+                sys.exit(f"ERROR: Expected 1 SMB climatology file matching "
+                         f"mesh {mesh_name} in {ctrl_atm_path}, found "
+                         f"{len(smb_files)}")
 
             temp_files = glob.glob(
                 os.path.join(ctrl_atm_path, '*temperature*.nc'))
             temp_files = [f for f in temp_files if 'gradient' not in f]
+            # Filter by mesh name to handle directories with multiple meshes
+            temp_files = [f for f in temp_files if mesh_name in f]
             if len(temp_files) == 1:
                 temp_fname = os.path.split(temp_files[0])[-1]
                 os.symlink(temp_files[0],
                            os.path.join(self.work_dir, temp_fname))
             else:
-                sys.exit(f"ERROR: Expected 1 temperature climatology file in "
-                         f"{ctrl_atm_path}, found {len(temp_files)}")
+                sys.exit(f"ERROR: Expected 1 temperature climatology file "
+                         f"matching mesh {mesh_name} in {ctrl_atm_path}, "
+                         f"found {len(temp_files)}")
 
             runoff_files = glob.glob(
                 os.path.join(ctrl_atm_path, '*runoff*.nc'))
+            # Filter by mesh name to handle directories with multiple meshes
+            runoff_files = [f for f in runoff_files if mesh_name in f]
             if len(runoff_files) == 1:
                 runoff_fname = os.path.split(runoff_files[0])[-1]
                 os.symlink(runoff_files[0],
@@ -152,6 +166,8 @@ def setup(self):  # noqa: C901
 
             smb_grad_files = glob.glob(
                 os.path.join(ctrl_atm_path, '*SMB_gradient*.nc'))
+            # Filter by mesh name to handle directories with multiple meshes
+            smb_grad_files = [f for f in smb_grad_files if mesh_name in f]
             smb_grad_fname = ''
             if len(smb_grad_files) == 1:
                 smb_grad_fname = os.path.split(smb_grad_files[0])[-1]
@@ -160,6 +176,8 @@ def setup(self):  # noqa: C901
 
             temp_grad_files = glob.glob(
                 os.path.join(ctrl_atm_path, '*temperature_gradient*.nc'))
+            # Filter by mesh name to handle directories with multiple meshes
+            temp_grad_files = [f for f in temp_grad_files if mesh_name in f]
             temp_grad_fname = ''
             if len(temp_grad_files) == 1:
                 temp_grad_fname = os.path.split(temp_grad_files[0])[-1]
@@ -174,30 +192,43 @@ def setup(self):  # noqa: C901
             # SMB forcing
             smb_search = os.path.join(atm_dir, '*SMB_*.nc')
             smb_list = glob.glob(smb_search)
+            print(f"    Found {len(smb_list)} SMB files before filtering")
             smb_list = [f for f in smb_list if 'gradient' not in f]
+            print(f"    Found {len(smb_list)} SMB files after removing "
+                  f"gradients")
+            # Filter by mesh name to handle directories with multiple meshes
+            smb_list = [f for f in smb_list if mesh_name in f]
+            print(f"    Found {len(smb_list)} SMB files after mesh filter: "
+                  f"{smb_list}")
             if len(smb_list) == 1:
                 smb_fname = os.path.split(smb_list[0])[-1]
                 os.symlink(smb_list[0],
                            os.path.join(self.work_dir, smb_fname))
             else:
-                sys.exit(f"ERROR: Expected 1 SMB file at {smb_search}, "
-                         f"found {len(smb_list)}: {smb_list}")
+                sys.exit(f"ERROR: Expected 1 SMB file matching mesh "
+                         f"{mesh_name} at {smb_search}, found "
+                         f"{len(smb_list)}: {smb_list}")
 
             # Temperature forcing
             temp_search = os.path.join(atm_dir, '*temperature_*.nc')
             temp_list = glob.glob(temp_search)
             temp_list = [f for f in temp_list if 'gradient' not in f]
+            # Filter by mesh name to handle directories with multiple meshes
+            temp_list = [f for f in temp_list if mesh_name in f]
             if len(temp_list) == 1:
                 temp_fname = os.path.split(temp_list[0])[-1]
                 os.symlink(temp_list[0],
                            os.path.join(self.work_dir, temp_fname))
             else:
-                sys.exit(f"ERROR: Expected 1 temperature file at "
-                         f"{temp_search}, found {len(temp_list)}")
+                sys.exit(f"ERROR: Expected 1 temperature file matching mesh "
+                         f"{mesh_name} at {temp_search}, found "
+                         f"{len(temp_list)}")
 
             # Runoff forcing (optional — may not exist for all experiments)
             runoff_search = os.path.join(atm_dir, '*runoff_*.nc')
             runoff_list = glob.glob(runoff_search)
+            # Filter by mesh name to handle directories with multiple meshes
+            runoff_list = [f for f in runoff_list if mesh_name in f]
             runoff_fname = ''
             if len(runoff_list) == 1:
                 runoff_fname = os.path.split(runoff_list[0])[-1]
@@ -207,6 +238,8 @@ def setup(self):  # noqa: C901
             # SMB gradient (lapse rate)
             smb_grad_search = os.path.join(atm_dir, '*SMB_gradient_*.nc')
             smb_grad_list = glob.glob(smb_grad_search)
+            # Filter by mesh name to handle directories with multiple meshes
+            smb_grad_list = [f for f in smb_grad_list if mesh_name in f]
             smb_grad_fname = ''
             if len(smb_grad_list) == 1:
                 smb_grad_fname = os.path.split(smb_grad_list[0])[-1]
@@ -217,6 +250,8 @@ def setup(self):  # noqa: C901
             temp_grad_search = os.path.join(atm_dir,
                                             '*temperature_gradient_*.nc')
             temp_grad_list = glob.glob(temp_grad_search)
+            # Filter by mesh name to handle directories with multiple meshes
+            temp_grad_list = [f for f in temp_grad_list if mesh_name in f]
             temp_grad_fname = ''
             if len(temp_grad_list) == 1:
                 temp_grad_fname = os.path.split(temp_grad_list[0])[-1]
@@ -226,34 +261,38 @@ def setup(self):  # noqa: C901
             # Thermal forcing
             tf_search = os.path.join(ocean_dir, '*thermal_forcing_*.nc')
             tf_list = glob.glob(tf_search)
+            # Filter by mesh name to handle directories with multiple meshes
+            tf_list = [f for f in tf_list if mesh_name in f]
             if len(tf_list) == 1:
                 tf_fname = os.path.split(tf_list[0])[-1]
                 os.symlink(tf_list[0],
                            os.path.join(self.work_dir, tf_fname))
             else:
-                sys.exit(f"ERROR: Expected 1 TF file at {tf_search}, "
-                         f"found {len(tf_list)}: {tf_list}")
+                sys.exit(f"ERROR: Expected 1 TF file matching mesh "
+                         f"{mesh_name} at {tf_search}, found "
+                         f"{len(tf_list)}: {tf_list}")
 
         # --- Find shelf collapse (calving) mask from ismip7_forcing
-        # fracture Path C, if requested ---
-        # historical, ctrl, and ocx experiments never use hydrofracture
-        # forcing, regardless of use_hydrofracture_forcing.
-        useCalvingMask = False
-        if (use_hydrofracture_forcing and
-                scenario not in ('historical', 'ctrl', 'ocx')):
+        # fracture Path C ---
+        # historical, ctrl, and ocx experiments do not use hydrofracture
+        # forcing.
+        if scenario not in ('historical', 'ctrl', 'ocx'):
             mask_search = os.path.join(forcing_dir, 'shelf_collapse',
                                        '*ice_shelf_collapse_mask_*.nc')
             mask_list = glob.glob(mask_search)
+            # Filter by mesh name to handle directories with multiple meshes
+            mask_list = [f for f in mask_list if mesh_name in f]
             if len(mask_list) == 1:
                 mask_fname = os.path.split(mask_list[0])[-1]
                 os.symlink(mask_list[0],
                            os.path.join(self.work_dir, mask_fname))
-                useCalvingMask = True
             else:
                 sys.exit(
-                    f"ERROR: use_hydrofracture_forcing is True but did not "
-                    f"find exactly 1 shelf collapse mask file at "
-                    f"{mask_search}: {mask_list}")
+                    f"ERROR: Did not find exactly 1 shelf collapse mask file "
+                    f"matching mesh {mesh_name} at {mask_search}: {mask_list}")
+        else:
+            # For historical/ctrl/ocx, set to empty - stream won't be used
+            mask_fname = ''
 
         # --- Set up streams ---
         # Determine forcing interval
@@ -282,6 +321,7 @@ def setup(self):  # noqa: C901
             'input_file_runoff_forcing': runoff_fname,
             'input_file_smb_gradient_forcing': smb_grad_fname,
             'input_file_temperature_gradient_forcing': temp_grad_fname,
+            'input_file_calving_mask_forcing': mask_fname,
             'forcing_interval_monthly': forcing_interval_monthly,
             'forcing_interval_annual': forcing_interval_annual,
         }
@@ -292,14 +332,6 @@ def setup(self):  # noqa: C901
             out_name='streams.landice',
             template_replacements=stream_replacements)
 
-        if useCalvingMask:
-            mask_stream_replacements = {
-                'input_file_calving_mask_forcing_name': mask_fname}
-            self.add_streams_file(
-                resource_location, 'streams.mask_calving',
-                out_name='streams.landice',
-                template_replacements=mask_stream_replacements)
-
         # --- Set up namelist ---
         self.add_namelist_file(
             resource_location, 'namelist.landice',
@@ -341,19 +373,6 @@ def setup(self):  # noqa: C901
                 out_name='streams.landice',
                 template_replacements=vM_stream_replacements)
 
-        # Mask calving options (ismip7_forcing fracture Path C), gated by
-        # hydrofracture vulnerability (MALI PR #187)
-        if useCalvingMask:
-            options = {
-                'config_apply_calving_mask': ".true.",
-                'config_restore_calving_front': ".false.",
-                'config_require_extensional_stresses_for_mask_calving':
-                ".true.",
-                'config_calving_fracture_toughness':
-                f'{calving_fracture_toughness}'}
-            self.add_namelist_options(options=options,
-                                      out_name='namelist.landice')
-
         # Sea-level model options
         if sea_level_model:
             slm_input_ice = section.get('slm_input_ice')