Ignore:
Timestamp:
May 21, 2025, 4:01:25 PM (4 weeks ago)
Author:
jbclement
Message:

PEM:
Following r3770, adaptation and improvements of Python scripts to visualize the layered deposits.
JBC

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/evolution/changelog.txt

    r3742 r3771  
    642642== 28/04/2025 == JBC
    643643Resulting modifications due to changes in the Mars PCM with r3739 and r3741.
     644
     645== 21/05/2025 == JBC
     646Revision of the layering structure and algorithm:
     647    - 'stratum' components are now expressed in height
     648    - deletion of the (redundant) thickness feature of 'stratum'
     649    - 'stratif' in the PEM main program is renamed 'deposits'
     650    - addition of several procedures to get useful information about a stratum (major component or thickness)
     651    - all subsurface layers are now represented in the layering structure by strata with negative top elevation
     652    - simplification of the different situations arising from the input tendencies
     653    - porosity is determined for the entire stratum (and not anymore for each component) based on dominant component
     654    - sublimation of CO2 and H2O ice is now handled simultaneously (more realistic) in a stratum
     655    - linking the layering algorithm with the PEM initilization/finalization regarding PCM data and with the PEM stopping criteria
     656    - making separate cases for glaciers vs layering management
     657    - H2O sublimation flux correction is now handled with the layering when a dust lag layer layer is growing
     658    - update of 'h2o_ice_depth' and 'zdqsdif' accordingly at the PEM end for the PCM
     659
     660== 21/05/2025 == JBC
     661Following r3770, adaptation and improvements of Python scripts to visualize the layered deposits.

  • trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py

    r3673 r3771  
    44
    55import os
     6import sys
    67import numpy as np
    78from netCDF4 import Dataset
    … …  
    910from scipy import interpolate
    1011
    11 #########################
    12 ### Parameters to fill in
    13 folder_path = "starts"
    14 base_name = "restartpem"
    15 dz = 0.1 # Discrization step of the reference grid for the elevation [m]
    16 infofile = 'info_PEM.txt'
    17 #########################
    18 
    19 
    20 ##############################################################################################
     12### Function to get inputs
     13def get_user_inputs():
     14    folder_path = input("Enter the folder path containing the NetCDF files (press the Enter key for default [starts]): ").strip()
     15    if not folder_path:
     16        folder_path = "starts"
     17    while not os.path.isdir(folder_path):
     18        print("Invalid folder path. Please try again.")
     19        folder_path = input("Enter the folder path containing the NetCDF files (press the Enter key for default [starts]): ").strip()
     20        if not folder_path:
     21            folder_path = "starts"
     22
     23    base_name = input("Enter the base name of the NetCDF files (press the Enter key for default [restartpem]): ").strip()
     24    if not base_name:
     25        base_name = "restartpem"
     26
     27    infofile = input("Enter the name of the PEM info file (press the Enter key for default [info_PEM.txt]): ").strip()
     28    if not infofile:
     29        infofile = "info_PEM.txt"
     30    while not os.path.isfile(infofile):
     31        print("Invalid file path. Please try again.")
     32        infofile = input("Enter the name of the PEM info file (press the Enter key for default [info_PEM.txt]): ").strip()
     33        if not infofile:
     34            infofile = "info_PEM.txt"
     35
     36    return folder_path, base_name, infofile
     37
    2138### Function to read the "startpem.nc" files and process their stratification data
    2239def process_files(folder_path,base_name):
    … …  
    3552    # Process each file and collect data
    3653    datasets = []
    37     max_top_elevation = 0
     54    min_base_elevation = -56943.759374550937 # Base elevation of the deepest subsurface layer
     55    max_top_elevation = 0.
    3856    max_nb_str = 0
    39     ngrid = Dataset(os.path.join(folder_path,base_name + "1.nc"),'r').dimensions['physical_points'].size # ngrid is the same for all files
    40     nslope = Dataset(os.path.join(folder_path,base_name + "1.nc"),'r').dimensions['nslope'].size # nslope is the same for all files
    41     longitude = Dataset(os.path.join(folder_path,base_name + "1.nc"),'r').variables['longitude'][:]
    42     latitude = Dataset(os.path.join(folder_path,base_name + "1.nc"),'r').variables['latitude'][:]
     57    with Dataset(os.path.join(folder_path, base_name + "1.nc"), 'r') as ds:
     58        ngrid = ds.dimensions['physical_points'].size # ngrid is the same for all files
     59        nslope = ds.dimensions['nslope'].size # nslope is the same for all files
     60        longitude = ds.variables['longitude'][:]
     61        latitude = ds.variables['latitude'][:]
     62
    4363    for i in range(1,nfile + 1):
    4464        file_path = os.path.join(folder_path,base_name + str(i) + ".nc")
    … …  
    5373        for k in range(1,nslope + 1):
    5474            slope_var_name = f"stratif_slope{k:02d}_top_elevation"
     75            min_base_elevation = min(min_base_elevation,np.min(ds.variables[slope_var_name][:]))
    5576            max_top_elevation = max(max_top_elevation,np.max(ds.variables[slope_var_name][:]))
    5677
    57     print(f"> number of files    = {nfile}")
    58     print(f"> ngrid              = {ngrid}")
    59     print(f"> nslope             = {nslope}")
    60     print(f"> max(nb_str_max)    = {max_nb_str}")
    61     print(f"> max(top_elevation) = {max_top_elevation}")
     78    print(f"> number of files     = {nfile}")
     79    print(f"> ngrid               = {ngrid}")
     80    print(f"> nslope              = {nslope}")
     81    print(f"> max(nb_str_max)     = {max_nb_str}")
     82    print(f"> min(base_elevation) = {min_base_elevation}")
     83    print(f"> max(top_elevation)  = {max_top_elevation}")
    6284
    6385    # Concatenate stratif variables with dimension 'nb_str_max' along the "Time" dimension
    … …  
    6890    stratif_dust = np.zeros((ngrid,nfile,nslope,max_nb_str))
    6991    stratif_pore = np.zeros((ngrid,nfile,nslope,max_nb_str))
     92    stratif_poreice = np.zeros((ngrid,nfile,nslope,max_nb_str))
    7093    for var_name in datasets[0].variables:
    7194        if 'top_elevation' in var_name:
    … …  
    7699                            stratif_heights[i,j,k,:datasets[j].variables[var_name].shape[1]] = datasets[j].variables[var_name][0,:,i]
    77100            print(f"Processed variable: {var_name}")
    78         elif 'co2ice_volfrac' in var_name:
     101        elif 'h_co2ice' in var_name:
    79102            for i in range(0,ngrid):
    80103                for j in range(0,nfile):
    … …  
    83106                            stratif_co2ice[i,j,k,:datasets[j].variables[var_name].shape[1]] = datasets[j].variables[var_name][0,:,i]
    84107            print(f"Processed variable: {var_name}")
    85         elif 'h2oice_volfrac' in var_name:
     108        elif 'h_h2oice' in var_name:
    86109            for i in range(0,ngrid):
    87110                for j in range(0,nfile):
    … …  
    90113                            stratif_h2oice[i,j,k,:datasets[j].variables[var_name].shape[1]] = datasets[j].variables[var_name][0,:,i]
    91114            print(f"Processed variable: {var_name}")
    92         elif 'dust_volfrac' in var_name:
     115        elif 'h_dust' in var_name:
    93116            for i in range(0,ngrid):
    94117                for j in range(0,nfile):
    … …  
    97120                            stratif_dust[i,j,k,:datasets[j].variables[var_name].shape[1]] = datasets[j].variables[var_name][0,:,i]
    98121            print(f"Processed variable: {var_name}")
    99         elif 'pore_volfrac' in var_name:
     122        elif 'h_pore' in var_name:
    100123            for i in range(0,ngrid):
    101124                for j in range(0,nfile):
    … …  
    103126                        if f'slope{k + 1:02d}' in var_name:
    104127                            stratif_pore[i,j,k,:datasets[j].variables[var_name].shape[1]] = datasets[j].variables[var_name][0,:,i]
     128            print(f"Processed variable: {var_name}")
     129        elif 'icepore_volfrac' in var_name:
     130            for i in range(0,ngrid):
     131                for j in range(0,nfile):
     132                    for k in range(0,nslope):
     133                        if f'slope{k + 1:02d}' in var_name:
     134                            stratif_poreice[i,j,k,:datasets[j].variables[var_name].shape[1]] = datasets[j].variables[var_name][0,:,i]
    105135            print(f"Processed variable: {var_name}")
    106136
    … …  
    110140
    111141    stratif_data = [stratif_heights,stratif_co2ice,stratif_h2oice,stratif_dust,stratif_pore]
    112     return stratif_data, max_top_elevation, longitude, latitude
     142
     143    while True:
     144        try:
     145            dz = float(input("Enter the discretization step of the reference grid for the elevation [m]: ").strip())
     146            if dz <= 0:
     147                print("Discretization step must be strictly positive!")
     148                continue
     149            if dz > max_top_elevation:
     150                print("Discretization step is higher than the maximum top elevation: please provide a correct value!")
     151                continue
     152            break
     153        except ValueError:
     154            print("Invalid value.")
     155    return stratif_data, min_base_elevation, max_top_elevation, longitude, latitude, dz
    113156
    114157### Function to interpolate the stratification data on a reference grid
    115 def interpolate_data(stratif_data,max_top_elevation,dz):
     158def interpolate_data(stratif_data,min_base_elevation,max_top_elevation,dz):
    116159    # Define the reference ref_grid
    117     ref_grid = np.arange(0,max_top_elevation,dz)
     160    ref_grid = np.arange(min_base_elevation,max_top_elevation,dz)
    118161    print(f"> Number of ref_grid points = {len(ref_grid)}")
    119162
    … …  
    143186        date_time = []
    144187        for line in file:
    145             data = list(map(int,line.split()))
     188            data = list(map(float,line.split()))
    146189            data_lines.append(data)
    147190            date_time.append(data[0])
    … …  
    150193
    151194### Processing
    152 stratif_data, max_top_elevation, longitude, latitude = process_files(folder_path,base_name)
    153 if dz > max_top_elevation:
    154     print('The discretization step is higher than the maximum top elevation: please provide a correct value!')
    155     exit()
    156 ref_grid, gridded_stratif_data = interpolate_data(stratif_data,max_top_elevation,dz)
     195folder_path, base_name, infofile = get_user_inputs()
     196stratif_data, min_base_elevation, max_top_elevation, longitude, latitude, dz = process_files(folder_path,base_name)
     197ref_grid, gridded_stratif_data = interpolate_data(stratif_data,min_base_elevation,max_top_elevation,dz)
    157198date_time = read_infofile(infofile)
    158199
    … …  
    168209        for ax in axes.flat:
    169210            time_mesh, elevation_mesh = np.meshgrid(date_time,ref_grid)
    170             #im = ax.imshow(np.transpose(gridded_stratif_data[iprop][0,:,0,:]),aspect = 'auto',cmap = 'viridis',origin = 'lower',extent = [date_time[0],date_time[-1],ref_grid[0],ref_grid[-1]],vmin = 0,vmax = 1)
     211            #im = ax.imshow(np.transpose(gridded_stratif_data[iprop][igr,:,isl,:]),aspect = 'auto',cmap = cmap,origin = 'lower',extent = [date_time[0],date_time[-1],ref_grid[0],ref_grid[-1]],vmin = 0,vmax = 1)
    171212            im = ax.pcolormesh(time_mesh,elevation_mesh,np.transpose(gridded_stratif_data[iprop][igr,:,isl,:]),cmap = cmap,shading = 'auto',vmin = 0,vmax = 1)
    172213            ax.set_title(subtitle[iprop])

  • trunk/LMDZ.COMMON/libf/evolution/deftank/visu_layering.py

    r3673 r3771  
    33####################################################################################
    44
     5import os
     6import sys
     7import numpy as np
    58import netCDF4 as nc
    6 import numpy as np
    79import matplotlib.pyplot as plt
    8 import sys
    9 import os.path
    10 
    11 ##############################
    12 ### Parameters to fill in
    13 filename = 'startpem.nc' #'starts/restartpem226.nc' # File name
    14 igrid = 1  # Grid point
    15 islope = 1 # Slope number
    16 istr = 4   # Stratum number
    17 ##############################
    18 
    19 
    20 ####################################################################################
    21 ### Open the NetCDF file
    22 if os.path.isfile(filename):
    23     nc_file = nc.Dataset(filename,'r')
    24 else:
    25     sys.exit('The file \"' + filename + '\" does not exist!')
    26 
    27 ### Get the dimensions
    28 Time = len(nc_file.dimensions['Time'])
    29 ngrid = len(nc_file.dimensions['physical_points'])
    30 nslope = len(nc_file.dimensions['nslope'])
    31 nb_str_max = len(nc_file.dimensions['nb_str_max'])
    32 if igrid > ngrid or igrid < 1:
    33     sys.exit('Asked grid point is not possible!')
    34 if islope > nslope or islope < 1:
    35     sys.exit('Asked slope number is not possible!')
    36 if istr > nb_str_max or istr < 1:
    37     sys.exit('Asked stratum number is not possible!')
    38 
    39 ### Get the stratification properties
    40 stratif_thickness = []
    41 stratif_top_elevation = []
    42 stratif_co2ice_volfrac = []
    43 stratif_h2oice_volfrac = []
    44 stratif_dust_volfrac = []
    45 stratif_pore_volfrac = []
    46 for i in range(1,nslope + 1):
    47     stratif_thickness.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_thickness'][:])
    48     stratif_top_elevation.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_top_elevation'][:])
    49     stratif_co2ice_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_co2ice_volfrac'][:])
    50     stratif_h2oice_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_h2oice_volfrac'][:])
    51     stratif_dust_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_dust_volfrac'][:])
    52     stratif_pore_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_pore_volfrac'][:])
    53 
    54 ### Display the data
    55 igrid = igrid - 1
    56 islope = islope - 1
    57 labels = 'CO2 ice', 'H2O ice', 'Dust', 'Pore'
    58 contents = np.zeros([4,len(stratif_top_elevation[islope][0,:,igrid]) + 1])
    59 height = np.zeros(len(stratif_top_elevation[islope][0,:,igrid]) + 1)
    60 height[0] = stratif_top_elevation[islope][0,0,:] - stratif_thickness[islope][0,0,:]
    61 height[1:] = stratif_top_elevation[islope][0,:,igrid]
    62 for i in range(len(stratif_top_elevation[islope][0,:,igrid])):
    63     contents[0,1 + i] = stratif_co2ice_volfrac[islope][0,i,igrid]
    64     contents[1,1 + i] = stratif_h2oice_volfrac[islope][0,i,igrid]
    65     contents[2,1 + i] = stratif_dust_volfrac[islope][0,i,igrid]
    66     contents[3,1 + i] = stratif_pore_volfrac[islope][0,i,igrid]
    67 contents[:,0] = contents[:,1]
    68 
    69 # Simple subplots for a layering
    70 fig, (ax1, ax2, ax3, ax4) = plt.subplots(1,4,layout = 'constrained',sharey = True)
    71 fig.suptitle('Simple content profiles for the layering')
    72 ax1.step(contents[0,:],height,where = 'post')
    73 ax2.step(contents[1,:],height,where = 'post')
    74 ax3.step(contents[2,:],height,where = 'post')
    75 ax4.step(contents[3,:],height,where = 'post')
    76 ax1.set_ylabel('Elevation [m]')
    77 ax1.set_xlabel('Volume fraction [m3/m3]')
    78 ax2.set_xlabel('Volume fraction [m3/m3]')
    79 ax3.set_xlabel('Volume fraction [m3/m3]')
    80 ax4.set_xlabel('Volume fraction [m3/m3]')
    81 ax1.set_title(labels[0])
    82 ax2.set_title(labels[1])
    83 ax3.set_title(labels[2])
    84 ax4.set_title(labels[3])
    85 plt.savefig('layering_simpleprofiles.png')
    86 
    87 # Content profiles for a layering
    88 plt.figure()
    89 plt.step(contents[0,:],height,where = 'post',color = 'r',label = labels[0])
    90 #plt.plot(contents[0,:],height,'o--',color = 'r',alpha = 0.3)
    91 plt.step(contents[1,:],height,where = 'post',color = 'b',label = labels[1])
    92 #plt.plot(contents[1,:],height,'o--',color = 'b',alpha = 0.3)
    93 plt.step(contents[2,:],height,where = 'post',color = 'y',label = labels[2])
    94 #plt.plot(contents[2,:],height,'o--',color = 'y',alpha = 0.3)
    95 plt.step(contents[3,:],height,where = 'post',color = 'g',label = labels[3])
    96 #plt.plot(contents[3,:],height,'o--',color = 'g',alpha = 0.3)
    97 plt.grid(axis='x', color='0.95')
    98 plt.grid(axis='y', color='0.95')
    99 plt.xlim(0,1)
    100 plt.xlabel('Volume fraction [m3/m3]')
    101 plt.ylabel('Elevation [m]')
    102 plt.title('Content profiles for the layering')
    103 plt.savefig('layering_simpleprofiles.png')
    104 
    105 # Stack content profiles for a layering
    106 fig, ax = plt.subplots()
    107 ax.fill_betweenx(height,0,contents[0,:],label = labels[0],color = 'r',step = 'pre')
    108 ax.fill_betweenx(height,contents[0,:],sum(contents[0:2,:]),label = labels[1],color = 'b',step = 'pre')
    109 ax.fill_betweenx(height,sum(contents[0:2,:]),sum(contents[0:3,:]),label = labels[2],color = 'y',step = 'pre')
    110 ax.fill_betweenx(height,sum(contents[0:3,:]),sum(contents),label = labels[3],color = 'g',step = 'pre')
    111 plt.vlines(x = 0.,ymin = height[0],ymax = height[len(height) - 1],color = 'k',linestyle = '-')
    112 plt.vlines(x = 1.,ymin = height[0],ymax = height[len(height) - 1],color = 'k',linestyle = '-')
    113 for i in range(len(height)):
    114     plt.hlines(y = height[i],xmin = 0.0,xmax = 1.0,color = 'k',linestyle = '--')
    115 ax.set_title('Stack content profiles for the layering')
    116 plt.xlabel('Volume fraction [m3/m3]')
    117 plt.ylabel('Elevation [m]')
    118 ax.legend(loc = 'center left',bbox_to_anchor = (1,0.5))
    119 plt.savefig('layering_stackprofiles.png')
    120 
    121 plt.show()
    122 
    123 ### Close the NetCDF file
    124 nc_file.close()
     10
     11
     12def get_int_input(prompt: str, min_val: int, max_val: int) -> int:
     13    """
     14    Prompt the user for an integer between min_val and max_val (inclusive).
     15    If min_val == max_val, return min_val without prompting.
     16    """
     17    if min_val == max_val:
     18        print(f"{prompt} (only possible value: {min_val})")
     19        return min_val
     20
     21    while True:
     22        try:
     23            value = int(input(f"{prompt} (integer between {min_val} and {max_val}): "))
     24            if min_val <= value <= max_val:
     25                return value
     26            print(f"Invalid value! Please enter a number between {min_val} and {max_val}.")
     27        except ValueError:
     28            print("Invalid input. Please enter a valid integer.")
     29
     30
     31def get_yes_no_input(prompt: str) -> bool:
     32    """
     33    Prompt the user with a yes/no question. Returns True for yes, False for no.
     34    """
     35    while True:
     36        choice = input(f"{prompt} (y/n): ").strip().lower()
     37        if choice in ['y', 'yes']:
     38            return True
     39        elif choice in ['n', 'no']:
     40            return False
     41        else:
     42            print("Please respond with y or n.")
     43
     44
     45def load_slope_variables(nc_dataset: nc.Dataset, slope_index: int) -> dict:
     46    """
     47    Load all relevant stratification variables for a given slope index (0-based).
     48    Returns a dict of NumPy arrays.
     49    """
     50    idx_str = str(slope_index + 1).zfill(2)
     51    vars_base = {
     52        'top_elev': f"stratif_slope{idx_str}_top_elevation",
     53        'h_co2': f"stratif_slope{idx_str}_h_co2ice",
     54        'h_h2o': f"stratif_slope{idx_str}_h_h2oice",
     55        'h_dust': f"stratif_slope{idx_str}_h_dust",
     56        'h_air': f"stratif_slope{idx_str}_h_pore",
     57        'volfrac': f"stratif_slope{idx_str}_icepore_volfrac",
     58    }
     59
     60    data = {}
     61    for key, var_name in vars_base.items():
     62        if var_name not in nc_dataset.variables:
     63            sys.exit(f"Error: Variable '{var_name}' not found in the NetCDF file.")
     64        data[key] = nc_dataset.variables[var_name][:]
     65    return data
     66
     67
     68def calculate_contents(data: dict, grid_index: int, exclude_subsurface: bool) -> (np.ndarray, np.ndarray, np.ndarray):
     69    """
     70    For a given slope's stratification data and a specified grid index (0-based),
     71    compute the height array, layer thicknesses, and layer-wise volume fractions.
     72    If exclude_subsurface is True, omit layers whose top elevation <= 0.
     73
     74    Returns:
     75      - height: 1D array of elevations (length = number_of_layers + 1)
     76      - contents: 2D array of shape (5, number_of_layers + 1), each row for a component
     77      - thicknesses: 1D array of layer thicknesses (length = number_of_layers)
     78    """
     79    top = data['top_elev']
     80    h_co2 = data['h_co2']
     81    h_h2o = data['h_h2o']
     82    h_dust = data['h_dust']
     83    h_air = data['h_air']
     84    volfrac = data['volfrac']
     85
     86    layers = top.shape[1]
     87    # Compute raw height and thickness arrays
     88    raw_height = np.zeros(layers + 1)
     89    raw_thickness = np.zeros(layers)
     90
     91    total_thickness0 = (
     92        h_co2[0, 0, grid_index]
     93        + h_h2o[0, 0, grid_index]
     94        + h_dust[0, 0, grid_index]
     95        + h_air[0, 0, grid_index]
     96    )
     97    raw_height[0] = top[0, 0, grid_index] - total_thickness0
     98    raw_height[1:] = top[0, :, grid_index]
     99
     100    for i in range(layers):
     101        raw_thickness[i] = (
     102            h_co2[0, i, grid_index]
     103            + h_h2o[0, i, grid_index]
     104            + h_dust[0, i, grid_index]
     105            + h_air[0, i, grid_index]
     106        )
     107
     108    include_mask = np.ones(layers, dtype=bool)
     109    if exclude_subsurface:
     110        include_mask = raw_height[1:] > 0
     111    if exclude_subsurface and not include_mask.any() and raw_height[0] <= 0:
     112        sys.exit("Error: No layers remain above the surface (elevation > 0).")
     113
     114    filt_layers = np.where(include_mask)[0]
     115    num_filt = len(filt_layers)
     116    height = np.zeros(num_filt + 1)
     117    thicknesses = np.zeros(num_filt)
     118
     119    if exclude_subsurface and raw_height[0] <= 0:
     120        height[0] = raw_height[filt_layers[0] + 1] - raw_thickness[filt_layers[0]]
     121    else:
     122        height[0] = raw_height[0]
     123
     124    for idx, layer_idx in enumerate(filt_layers):
     125        height[idx + 1] = raw_height[layer_idx + 1]
     126        thicknesses[idx] = raw_thickness[layer_idx]
     127
     128    contents = np.zeros((5, num_filt + 1))
     129    for idx, layer_idx in enumerate(filt_layers):
     130        thickness = raw_thickness[layer_idx]
     131        if thickness <= 0:
     132            continue
     133        co2 = h_co2[0, layer_idx, grid_index] / thickness
     134        h2o = h_h2o[0, layer_idx, grid_index] / thickness
     135        dust = h_dust[0, layer_idx, grid_index] / thickness
     136        air = h_air[0, layer_idx, grid_index] * (1.0 - volfrac[0, layer_idx, grid_index]) / thickness
     137        pore = h_air[0, layer_idx, grid_index] * volfrac[0, layer_idx, grid_index] / thickness
     138        contents[:, idx + 1] = [co2, h2o, dust, air, pore]
     139
     140    if num_filt > 0:
     141        contents[:, 0] = contents[:, 1]
     142
     143    return height, contents, thicknesses
     144
     145
     146def plot_profiles(height: np.ndarray, contents: np.ndarray, thicknesses: np.ndarray, labels: tuple[str, ...]) -> None:
     147    """
     148    Create and save plots:
     149      1. Simple step profiles (each component in its own subplot)
     150      2. Overlaid step profiles (all components on one plot)
     151      3. Stacked fill-between profile with relative fractions
     152    """
     153    n_components = contents.shape[0]
     154    colors = ['r', 'b', 'y', 'violet', 'c']
     155
     156    # 1. Simple subplots
     157    fig, axes = plt.subplots(1, n_components, sharey=True, constrained_layout=True)
     158    fig.suptitle('Simple Content Profiles for Stratification')
     159
     160    for idx, ax in enumerate(axes):
     161        ax.step(contents[idx, :], height, where='post', color=colors[idx])
     162        ax.set_xlabel('Volume Fraction [m^3/m^3]')
     163        ax.set_title(labels[idx])
     164        if idx == 0:
     165            ax.set_ylabel('Elevation [m]')
     166        # Add major and minor grids on y-axis
     167        ax.grid(which='major', axis='y', linestyle='--', linewidth=0.5, color='0.8')
     168        ax.minorticks_on()
     169        ax.grid(which='minor', axis='y', linestyle=':', linewidth=0.3, color='0.9')
     170    plt.savefig('layering_simple_profiles.png')
     171
     172    # 2. Overlaid step profiles
     173    plt.figure()
     174    for idx, label in enumerate(labels):
     175        plt.step(contents[idx, :], height, where='post', color=colors[idx], label=label)
     176    plt.grid(axis='x', color='0.95')
     177    plt.grid(axis='y', color='0.95')
     178    plt.minorticks_on()
     179    plt.grid(which='minor', axis='y', linestyle=':', linewidth=0.3, color='0.9')
     180    plt.xlim(0, 1)
     181    plt.xlabel('Volume Fraction [m^3/m^3]')
     182    plt.ylabel('Elevation [m]')
     183    plt.title('Content Profiles for Stratification')
     184    plt.legend()
     185    plt.savefig('layering_overlaid_profiles.png')
     186
     187    # 3. Stacked fill-between profile with relative fractions
     188    plt.figure()
     189    cumulative = np.zeros_like(contents[0, :])
     190    for idx, label in enumerate(labels):
     191        plt.fill_betweenx(
     192            height,
     193            cumulative,
     194            cumulative + contents[idx, :],
     195            step='pre',
     196            label=label,
     197            color=colors[idx]
     198        )
     199        cumulative += contents[idx, :]
     200    plt.vlines(x=0., ymin=height[0], ymax=height[-1], color='k', linestyle='-')
     201    plt.vlines(x=1., ymin=height[0], ymax=height[-1], color='k', linestyle='-')
     202    for h in height:
     203        plt.hlines(y=h, xmin=0.0, xmax=1.0, color='k', linestyle='--', linewidth=0.5)
     204    plt.xlabel('Volume Fraction [m^3/m^3]')
     205    plt.ylabel('Elevation [m]')
     206    plt.title('Stacked Content Profiles for Stratification')
     207    plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
     208    plt.tight_layout()
     209    plt.savefig('layering_stacked_profiles.png')
     210
     211
     212def plot_stratum_layer(contents: np.ndarray, istr_index: int, labels: tuple[str, ...]) -> None:
     213    """
     214    Plot detailed composition of the specified stratum index (0-based) as a bar chart.
     215    """
     216    layer_fractions = contents[:, istr_index + 1]
     217    plt.figure()
     218    x_positions = np.arange(len(labels))
     219    plt.bar(x_positions, layer_fractions, color=['r', 'b', 'y', 'violet', 'c'])
     220    plt.xticks(x_positions, labels, rotation=45, ha='right')
     221    plt.ylim(0, 1)
     222    plt.ylabel('Volume Fraction [m^3/m^3]')
     223    plt.title(f'Composition of Stratum {istr_index + 1}')
     224    # Add major and minor grids on y-axis for the histogram
     225    plt.grid(which='major', axis='y', linestyle='--', linewidth=0.5, color='0.8')
     226    plt.minorticks_on()
     227    plt.grid(which='minor', axis='y', linestyle=':', linewidth=0.3, color='0.9')
     228    plt.tight_layout()
     229    plt.savefig(f'stratum_{istr_index+1}_composition.png')
     230
     231
     232def main():
     233    if len(sys.argv) > 1:
     234        filename = sys.argv[1]
     235    else:
     236        filename = input("Enter the NetCDF file name: ").strip()
     237
     238    if not os.path.isfile(filename):
     239        sys.exit(f"Error: File '{filename}' does not exist.")
     240
     241    with nc.Dataset(filename, 'r') as dataset:
     242        required_dims = ['Time', 'physical_points', 'nslope', 'nb_str_max']
     243        for dim in required_dims:
     244            if dim not in dataset.dimensions:
     245                sys.exit(f"Error: Missing dimension '{dim}' in file.")
     246
     247        ngrid = len(dataset.dimensions['physical_points'])
     248        nslope = len(dataset.dimensions['nslope'])
     249        nb_str_max = len(dataset.dimensions['nb_str_max'])
     250
     251        print(f"File '{filename}' opened successfully.")
     252        print(f"Number of grid points (physical_points): {ngrid}")
     253        print(f"Number of slopes: {nslope}")
     254        print(f"Maximum number of strata per slope: {nb_str_max}")
     255
     256        igrid_input = get_int_input(
     257            "Enter grid point number", 1, ngrid
     258        ) - 1
     259        islope_input = get_int_input(
     260            "Enter slope number", 1, nslope
     261        ) - 1
     262
     263        show_subsurface = get_yes_no_input("Show subsurface layers?")
     264        exclude_sub = not show_subsurface
     265
     266        # Load data for the chosen slope to determine number of surface strata
     267        slope_data = load_slope_variables(dataset, islope_input)
     268
     269        # Compute raw heights to count strata above surface
     270        top = slope_data['top_elev']
     271        layers = top.shape[1]
     272        raw_height = np.zeros(layers + 1)
     273        # Compute initial subsurface bottom
     274        h_co2 = slope_data['h_co2']
     275        h_h2o = slope_data['h_h2o']
     276        h_dust = slope_data['h_dust']
     277        h_air = slope_data['h_air']
     278        total_thickness0 = (
     279            h_co2[0, 0, igrid_input]
     280            + h_h2o[0, 0, igrid_input]
     281            + h_dust[0, 0, igrid_input]
     282            + h_air[0, 0, igrid_input]
     283        )
     284        raw_height[0] = top[0, 0, igrid_input] - total_thickness0
     285        raw_height[1:] = top[0, :, igrid_input]
     286
     287        include_mask = np.ones(layers, dtype=bool)
     288        if exclude_sub:
     289            include_mask = raw_height[1:] > 0
     290        # Count number of strata above surface
     291        filt_layers = np.where(include_mask)[0]
     292        nb_str_surf_max = len(filt_layers)
     293        if not show_subsurface:
     294            if nb_str_surf_max == 0:
     295                print("No stratum layers remain after filtering. Cannot proceed.")
     296                return
     297
     298        # Prompt for stratum index based on surface strata count
     299        istr_input = get_int_input(
     300            "Enter stratum index for detailed plot", 1, nb_str_surf_max if exclude_sub else nb_str_max
     301        ) - 1
     302
     303        height_arr, contents_arr, thicknesses_arr = calculate_contents(
     304            slope_data, igrid_input, exclude_sub
     305        )
     306
     307        component_labels = ("CO2 Ice", "H2O Ice", "Dust", "Air", "Pore ice")
     308
     309        plot_profiles(height_arr, contents_arr, thicknesses_arr, component_labels)
     310        plot_stratum_layer(contents_arr, istr_input, component_labels)
     311
     312        # Show all figures at once
     313        plt.show()
     314
     315
     316if __name__ == '__main__':
     317    main()
Note: See TracChangeset for help on using the changeset viewer.