Merge remote-tracking branch 'origin/devel' into devel

Author: rurata

bin/MarsInterp.py

@@ -141,8 +141,8 @@ def wrapper(*args, **kwargs):
         f"(zstd), or altitude above ground level (zagl).\nWorks on "
         f"'daily', 'average', and 'diurn' files.\n"
         f"{Green}Example:\n"
-        f"> MarsInterp 01336.atmos_average.nc\n"
         f"> MarsInterp 01336.atmos_average.nc -t pstd\n"
+        f"> MarsInterp 01336.atmos_average.nc -t pstd -v pstd_default\n"
         f"{Nclr}\n\n"
     )
 )
@@ -300,7 +300,7 @@ def main():
     # Load all of the netcdf files
     file_list    = file_list = [f.name for f in args.input_file]
     interp_type  = args.interp_type  # e.g. pstd
-    custom_level = args.vertical_grid # e.g. p44
+    custom_level = args.vertical_grid # e.g. pstd_default
     grid_out     = args.print_grid
 
     # Create a namespace with numpy available
@@ -322,7 +322,7 @@ def main():
         if custom_level:
             lev_in = eval(f"np.array({custom_level})", namespace)
         else:
-            lev_in = np.array(namespace['pstd_default'])
+            lev_in = eval("np.array(pstd_default)", namespace)
 
     # =========================== zstd ===========================
     elif interp_type == "zstd":

bin/MarsVars.py

@@ -254,6 +254,18 @@ def wrapper(*args, **kwargs):
         ['Vg_sed', 'w'],
         ['pfull', 'pstd', 'zstd', 'zagl']
     ],
+    'dustref_per_pa': [
+        "Dust opacity per Pascal [dustref/delp]",
+        'op/Pa',
+        ['dustref', 'delp'],
+        ['pfull']
+    ],
+    'dustref_per_km': [
+        "Dust opacity per kilometer [dustref/delz]",
+        'op/km',
+        ['dustref', 'delz'],
+        ['pfull']
+    ],
     'wdir': [
         "Wind direction",
         'degree',
@@ -1173,6 +1185,59 @@ def compute_w_net(Vg, wvar):
     return w_net
 
 
+
+# =====================================================================
+def compute_dustref_per_pa(dustref, delp):
+    """
+    Computes visible dust opacity per Pascal from dustref and delp. 
+    dustref is visible dust opacity per level (model layer).
+
+    opacity per Pa = opacity per layer / layer thickness in Pa::
+
+        dustref_per_pa = dustref/delp
+
+    [Courtney Batterson, 2025]
+
+    :param dustref: Visible dust opacity [op/model layer]
+    :type  dustref: array [time, lev, lat, lon]
+    :param delp: Layer thickness [Pa]
+    :type  delp: array [time, lev, lat, lon]
+    :return: `dustref_per_pa` Visible dust opacity [op/Pa]
+    :rtype:  array [time, lev, lat, lon]
+    :raises ValueError: If the input dimensions are not compatible
+    :raises TypeError: If the input types are not compatible
+    :raises Exception: If any other error occurs
+    """
+
+    dustref_per_pa = dustref/delp
+    return dustref_per_pa
+
+# =====================================================================
+def compute_dustref_per_z(dustref, delz):
+    """
+    Computes visible dust opacity per kilometer from dustref and delz. 
+    dustref is visible dust opacity per level (model layer).
+
+    opacity per km = opacity per layer / layer thickness in m * 1000::
+
+        dustref_per_z = dustref/delz*1000
+
+    [Courtney Batterson, 2025]
+
+    :param dustref: Visible dust opacity [op/model layer]
+    :type  dustref: array [time, lev, lat, lon]
+    :param delz: Layer thickness [m]
+    :type  delz: array [time, lev, lat, lon]
+    :return: `dustref_per_z` Visible dust opacity [op/km]
+    :rtype:  array [time, lev, lat, lon]
+    :raises ValueError: If the input dimensions are not compatible
+    :raises TypeError: If the input types are not compatible
+    :raises Exception: If any other error occurs
+    """
+
+    dustref_per_z = dustref/delz*1000
+    return dustref_per_z
+
 # =====================================================================
 def compute_theta(p_3D, ps, T, f_type):
     """
@@ -1906,9 +1971,9 @@ def process_add_variables(file_name, add_list, master_list, debug=False):
                       f"file:{Nclr}")
                 for var, actual_var in already_in_file:
                     if var == actual_var:
-                        print(f"{Yellow}  - {var}{Nclr}")
+                        print(f"{Yellow} - {var}{Nclr}")
                     else:
-                        print(f"{Yellow}  - {var} (as '{actual_var}'){Nclr}")
+                        print(f"{Yellow} - {var} (as '{actual_var}'){Nclr}")
             return
 
 
@@ -2101,6 +2166,16 @@ def add_with_dependencies(var):
                 wvar = f.variables["w"][:]
                 OUT = compute_w_net(Vg, wvar)
 
+            if var == "dustref_per_pa":
+                dustref = f.variables["dustref"][:]
+                delp = f.variables["delp"][:]
+                OUT = compute_dustref_per_pa(dustref, delp)
+            
+            if var == "dustref_per_z":
+                dustref = f.variables["dustref"][:]
+                delz = f.variables["delz"][:]
+                OUT = compute_dustref_per_z(dustref, delz)
+            
             if var == "pfull3D":
                 OUT = p_3D
 

docs/source/description.rst

@@ -300,8 +300,6 @@ Open ``~/.amesgcm_profile`` with any text editor to see customizable grid defini
          1.0e-02, 5.0e-03, 3.0e-03, 5.0e-04, 3.0e-04, 1.0e-04, 5.0e-05,
          3.0e-05, 1.0e-05]
 
-    phalf_mb=[50]
-
 Use your custom grid with the ``[-v --vertical_grid]`` argument:
 
 .. code-block:: bash