Adding script for computing exf_albedo

src/mitgcm_inputs/__main__.py

@@ -17,6 +17,9 @@
 from mitgcm_inputs.cosmetic_mask import COMMAND_NAME as CMASK_COMMAND_NAME
 from mitgcm_inputs.cosmetic_mask import main as cosmetic_mask_main
 from mitgcm_inputs.cosmetic_mask import sub_arguments as cmask_sub_arguments
+from mitgcm_inputs.exf_albedo import COMMAND_NAME as EXF_ALBEDO_COMMAND_NAME
+from mitgcm_inputs.exf_albedo import main as exf_albedo_main
+from mitgcm_inputs.exf_albedo import sub_arguments as exf_albedo_sub_arguments
 from mitgcm_inputs.k_extinction import COMMAND_NAME as KEXT_COMMAND_NAME
 from mitgcm_inputs.k_extinction import main as kext_main
 from mitgcm_inputs.k_extinction import sub_arguments as kext_sub_arguments
@@ -81,6 +84,7 @@ def argument(sys_argv=None):
     rbcs_sub_arguments(subparsers)
     ob_indices_sub_arguments(subparsers)
     cmask_sub_arguments(subparsers)
+    exf_albedo_sub_arguments(subparsers)
 
     subparser = subparsers.add_parser(
         "FLUXES",
@@ -193,6 +197,7 @@ def execute_flux_commands(args: argparse.Namespace):
     KEXT_COMMAND_NAME: kext_main,
     RBCS_COMMAND_NAME: rbcs_main,
     OB_INDICES_COMMAND_NAME: ob_indices_main,
+    EXF_ALBEDO_COMMAND_NAME: exf_albedo_main,
     "FLUXES": execute_flux_commands,
 }
 

src/mitgcm_inputs/cosmetic_mask/__init__.py

@@ -83,7 +83,7 @@ def main(args: argparse.Namespace) -> int:
     else:
         file_with_river = args.mask
 
-    with xr.open_dataset(file_with_river) as river_ds:
+    with xr.open_dataset(file_with_river, engine="netcdf4") as river_ds:
         river_map = river_ds["rivers"].load()
 
     river_mask = MaskWithRivers(

src/mitgcm_inputs/exf_albedo/__init__.py

@@ -0,0 +1,65 @@
+import argparse
+import logging
+
+import xarray as xr
+from bitsea.utilities.argparse_types import existing_file_path
+from bitsea.utilities.argparse_types import path_inside_an_existing_dir
+
+if __name__ == "__main__":
+    LOGGER = logging.getLogger()
+else:
+    LOGGER = logging.getLogger(__name__)
+
+
+COMMAND_NAME = "exf_albedo"
+
+
+def sub_arguments(subparser):
+    parser = subparser.add_parser(
+        COMMAND_NAME, help="Compute the EXF albedo for a specific domain"
+    )
+
+    parser.add_argument(
+        "-m",
+        "--mask",
+        required=True,
+        type=existing_file_path,
+        help="The meshmask file that defines the domain",
+    )
+
+    parser.add_argument(
+        "-o",
+        "--output",
+        required=True,
+        type=path_inside_an_existing_dir,
+        help="Path of the output file",
+    )
+
+
+def compute_albedo(average_lat: float) -> float:
+    return 0.06 + 0.0022 * (average_lat - 37)
+
+
+def main(args: argparse.Namespace) -> int:
+    if args.cmd != COMMAND_NAME:
+        LOGGER.error(
+            "EXF albedo command has been invoked with command: %s", args.cmd
+        )
+        return 1
+
+    LOGGER.debug("Opening file %s", args.mask)
+    with xr.open_dataset(args.mask, engine="netcdf4") as ds:
+        min_lat = float(ds.attrs["domain_minimum_latitude"])
+        max_lat = float(ds.attrs["domain_maximum_latitude"])
+    LOGGER.debug("Domain latitude range: %s - %s", min_lat, max_lat)
+
+    average_lat = (min_lat + max_lat) / 2
+    albedo = compute_albedo(average_lat)
+
+    LOGGER.debug("Writing output to file %s", args.output)
+    with open(args.output, "w") as f:
+        f.write(f"exf_albedo={albedo:.3f}\n")
+
+    LOGGER.info("Execution completed!")
+
+    return 0

src/mitgcm_inputs/k_extinction/k_extinction.py

@@ -34,7 +34,7 @@ def compute_k_extinction(
     LOGGER.debug("Meshmask shape: %s", meshmask.shape)
 
     LOGGER.debug("Reading data file: %s", data_file)
-    climatological_data = xr.load_dataset(data_file)
+    climatological_data = xr.load_dataset(data_file, engine="netcdf4")
 
     data_lon_min = np.min(climatological_data.longitude.values)
     data_lon_max = np.max(climatological_data.longitude.values)

src/mitgcm_inputs/ob_indices/__init__.py

@@ -114,7 +114,7 @@ def main(args: argparse.Namespace) -> int:
     # If there are rivers in this domain
     if args.rivers_positions is not None:
         with xr.open_dataset(
-            file_with_river, mask_and_scale=False
+            file_with_river, mask_and_scale=False, engine="netcdf4"
         ) as river_ds:
             river_map = river_ds["rivers"].load()
 

src/mitgcm_inputs/rbcs/rbcs_gen.py

@@ -107,7 +107,7 @@ def get_spatial_description_from_mit_static(input_file: Path) -> xr.Dataset:
     Returns:
         A spatial description of our domain
     """
-    with xr.open_dataset(input_file) as ds:
+    with xr.open_dataset(input_file, engine="netcdf4") as ds:
         spatial_description = xr.Dataset(
             data_vars={
                 "bathymetry": (("latitude", "longitude"), ds.Depth.values),

src/mitgcm_inputs/tools/read_mesh_mask.py

@@ -6,7 +6,7 @@
 
 
 def read_mesh_mask(mask_file_path: Path):
-    with xr.open_dataset(mask_file_path) as ds:
+    with xr.open_dataset(mask_file_path, engine="netcdf4") as ds:
         if "mer_mesh_mask_version" not in ds.attrs:
             return Mask.from_file(mask_file_path)
         latitude = ds.latitude.values