Changeset 4192

Timestamp:

Apr 20, 2026, 3:58:44 PM (3 weeks ago)

Author:

jbclement

Message:

PEM:
Codebase for pore ice representation and evolution in the layered deposits:

• replace single pore ice volume fraction with coefficients of polymial expansion to describe the stratum pore ice profile;
• add the formation date in the stratum information;
• make default initialization for 'stratum';
• centralize the definition of stratum types;
• merge of similar strata;
• add interpolation and polynomial fitting tools;
• remove subsurface strata.

JBC

Location:

trunk/LMDZ.COMMON

Files:

• libf/evolution/backup.F90 (modified) (1 diff)
• libf/evolution/changelog.txt (modified) (1 diff)
• libf/evolution/climate_init.F90 (modified) (7 diffs)
• libf/evolution/climate_rec.F90 (modified) (4 diffs)
• libf/evolution/deftank/visu_layering.py (modified) (5 diffs)
• libf/evolution/deftank/visu_layering_evo.py (modified) (8 diffs)
• libf/evolution/display.F90 (modified) (3 diffs)
• libf/evolution/layered_deposits.F90 (modified) (44 diffs)
• libf/evolution/maths.F90 (modified) (3 diffs)
• libf/evolution/pem.F90 (modified) (2 diffs)
• libf/evolution/soil.F90 (modified) (2 diffs)
• libf/evolution/subsurface_ice.F90 (modified) (18 diffs)
• libf/misc/job.F90 (modified) (3 diffs)
• makelmdz_fcm (modified) (2 diffs)

trunk/LMDZ.COMMON/libf/evolution/backup.F90

@@ -96 +96 @@
 use tracers,          only: build4PCM_tracers, nq
 use ice_table,        only: build4PCM_ssice
-use climate_rec,   only: write_restart, write_restartfi, write_restartevo
+use climate_rec,      only: write_restart, write_restartfi, write_restartevo
 use layered_deposits, only: layering
 

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

@@ -977 +977 @@
 - Rename "clim_state_init"/"clim_state_rec" into "climate_init"/"climate_rec".
 - Add routines to handle conservative remapping in 1d with slope limiters, positivity and boundedness enforcements.
+
+== 20/04/2026 == JBC
+Codebase for pore ice representation and evolution in the layered deposits:
+    - replace single pore ice volume fraction with coefficients of polymial expansion to describe the stratum pore ice profile;
+    - add the formation date in the stratum information;
+    - make default initialization for 'stratum';
+    - centralize the definition of stratum types;
+    - merge of similar strata;
+    - add interpolation and polynomial fitting tools;
+    - remove subsurface strata.

trunk/LMDZ.COMMON/libf/evolution/climate_init.F90

@@ -321 +321 @@
 use stoppage,           only: stop_clean
 use geometry,           only: ngrid, nslope, nsoil, nsoil_PCM
-use evolution,          only: dt
+use evolution,          only: dt, pem_ini_date, n_yr_sim
 use physics,            only: mugaz, r, alpha_clap_h2o, beta_clap_h2o
 use frost,              only: h2ofrost_PCM, h2o_frost4PCM, co2frost_PCM, co2_frost4PCM
@@ -331 +331 @@
 use sorption,           only: do_sorption, evolve_regolith_adsorption
 use tracers,            only: mmol, iPCM_qh2o
-use layered_deposits,   only: layering, do_layering, array2map, ini_layering, add_stratum
+use layered_deposits,   only: layering, do_layering, array2map, ini_layering, add_stratum, ncoef_poreice
 use surf_ice,           only: rho_co2ice, rho_h2oice
 use slopes,             only: subslope_dist, def_slope_mean
@@ -361 +361 @@
 ! LOCAL VARIABLES
 ! ---------------
-logical(k4)                               :: here
-integer(di)                               :: i, islope, k, nb_str_max, nsoil_startevo
-real(dp)                                  :: delta           ! Depth of the interface regolith/breccia, breccia/bedrock [m]
-real(dp), dimension(ngrid,nsoil,nslope)   :: TI_startevo     ! Soil thermal inertia saved in the startevo [SI]
-real(dp), dimension(ngrid,nsoil,nslope)   :: tsoil_startevo  ! Soil temperature saved in the startevo [K]
-real(dp), dimension(:,:,:,:), allocatable :: layerings_array ! Array for layerings
+logical(k4)                                   :: here
+integer(di)                                   :: i, islope, k, nb_str_max, nsoil_startevo
+real(dp)                                      :: delta           ! Depth of the interface regolith/breccia, breccia/bedrock [m]
+real(dp), dimension(ngrid,nsoil,nslope)       :: TI_startevo     ! Soil thermal inertia saved in the startevo [SI]
+real(dp), dimension(ngrid,nsoil,nslope)       :: tsoil_startevo  ! Soil temperature saved in the startevo [K]
+real(dp), dimension(:,:,:,:), allocatable     :: layerings_array ! Array for layerings
+real(dp), dimension(ncoef_poreice), parameter :: coef_poreice_0 = 0._dp
 
 ! CODE
@@ -489 +490 @@
     ! Layering
     if (do_layering) then
-        call print_msg('layerings_map is initialized with sub-surface strata.',LVL_WRN)
+        call print_msg('layerings_map is initialized with a surface reference stratum.',LVL_WRN)
         call print_msg("Ice is added with 'h2oice_huge_ini' where 'watercaptag' is true and otherwise with 'perennial_co2ice' found in the PCM.",LVL_WRN)
         do i = 1,ngrid
@@ -495 +496 @@
                 do islope = 1,nslope
                     call ini_layering(layerings_map(i,islope))
-                    call add_stratum(layerings_map(i,islope),1.05_dp*h2oice_huge_ini/rho_h2oice,0._dp,h2oice_huge_ini/rho_h2oice,0.05_dp*h2oice_huge_ini/rho_h2oice,0._dp,0._dp)
+                    call add_stratum(layerings_map(i,islope),1.05_dp*h2oice_huge_ini/rho_h2oice,0._dp,h2oice_huge_ini/rho_h2oice, &
+                                     0.05_dp*h2oice_huge_ini/rho_h2oice,0._dp,coef_poreice_0,pem_ini_date + n_yr_sim)
                 end do
             else
                 do islope = 1,nslope
                     call ini_layering(layerings_map(i,islope))
-                    if (co2_perice_PCM(i,islope) > 0.) call add_stratum(layerings_map(i,islope),1.05_dp*co2_perice_PCM(i,islope)/rho_co2ice,co2_perice_PCM(i,islope)/rho_co2ice,0._dp,0.05_dp*co2_perice_PCM(i,islope)/rho_co2ice,0._dp,0._dp)
+                    if (co2_perice_PCM(i,islope) > 0.) then
+                        call add_stratum(layerings_map(i,islope),1.05_dp*co2_perice_PCM(i,islope)/rho_co2ice,co2_perice_PCM(i,islope)/rho_co2ice,0._dp, &
+                                         0.05_dp*co2_perice_PCM(i,islope)/rho_co2ice,0._dp,coef_poreice_0,pem_ini_date + n_yr_sim)
+                    end if
                 end do
             end if
@@ -567 +572 @@
 if (do_layering) then
     call get_dim_nc('nb_str_max',nb_str_max)
-    allocate(layerings_array(ngrid,nslope,nb_str_max,6))
+    allocate(layerings_array(ngrid,nslope,nb_str_max,6 + ncoef_poreice))
     layerings_array = 0.
     call get_var_nc('stratif_top_elevation',layerings_array(:,:,:,1))
@@ -574 +579 @@
     call get_var_nc('stratif_h_dust',layerings_array(:,:,:,4))
     call get_var_nc('stratif_h_pore',layerings_array(:,:,:,5))
-    call get_var_nc('stratif_poreice_volfrac',layerings_array(:,:,:,6))
+    call get_var_nc('stratif_poreice_coef1',layerings_array(:,:,:,6))
+    call get_var_nc('stratif_poreice_coef2',layerings_array(:,:,:,7))
+    call get_var_nc('stratif_poreice_coef3',layerings_array(:,:,:,8))
+    call get_var_nc('stratif_poreice_coef4',layerings_array(:,:,:,9))
+    call get_var_nc('stratif_date',layerings_array(:,:,:,10))
     call array2map(layerings_array,layerings_map)
     deallocate(layerings_array)

trunk/LMDZ.COMMON/libf/evolution/climate_rec.F90

@@ -395 +395 @@
     call check_nc(nf90_put_att(ncid,varid,'title','Layering pore height'),'putting title attribute for stratif_h_pore')
 
-    call check_nc(nf90_def_var(ncid,'stratif_poreice_volfrac',nf90_double,(/dim_ngrid,dim_nb_str,dim_nslope,dim_time/),varid),'defining variable stratif_poreice_volfrac')
-    call check_nc(nf90_put_att(ncid,varid,'title','Layering ice pore volume fraction'),'putting title attribute for stratif_poreice_volfrac')
+    call check_nc(nf90_def_var(ncid,'stratif_poreice_coef1',nf90_double,(/dim_ngrid,dim_nb_str,dim_nslope,dim_time/),varid),'defining variable stratif_poreice_coef1')
+    call check_nc(nf90_put_att(ncid,varid,'title','Layering pore-ice coefficient #1'),'putting title attribute for stratif_poreice_coef1')
+
+    call check_nc(nf90_def_var(ncid,'stratif_poreice_coef2',nf90_double,(/dim_ngrid,dim_nb_str,dim_nslope,dim_time/),varid),'defining variable stratif_poreice_coef2')
+    call check_nc(nf90_put_att(ncid,varid,'title','Layering pore-ice coefficient #2'),'putting title attribute for stratif_poreice_coef2')
+
+    call check_nc(nf90_def_var(ncid,'stratif_poreice_coef3',nf90_double,(/dim_ngrid,dim_nb_str,dim_nslope,dim_time/),varid),'defining variable stratif_poreice_coef3')
+    call check_nc(nf90_put_att(ncid,varid,'title','Layering pore-ice coefficient #3'),'putting title attribute for stratif_poreice_coef3')
+
+    call check_nc(nf90_def_var(ncid,'stratif_poreice_coef4',nf90_double,(/dim_ngrid,dim_nb_str,dim_nslope,dim_time/),varid),'defining variable stratif_poreice_coef4')
+    call check_nc(nf90_put_att(ncid,varid,'title','Layering pore-ice coefficient #4'),'putting title attribute for stratif_poreice_coef4')
+
+    call check_nc(nf90_def_var(ncid,'stratif_date',nf90_double,(/dim_ngrid,dim_nb_str,dim_nslope,dim_time/),varid),'defining variable stratif_date')
+    call check_nc(nf90_put_att(ncid,varid,'title','Layering stratum formation date'),'putting title attribute for stratif_date')
 end if
 
@@ -436 +448 @@
 use geometry,         only: ngrid, nslope
 use evolution,        only: pem_ini_date, n_yr_sim
-use layered_deposits, only: layering, do_layering, map2array, get_nb_str_max
+use layered_deposits, only: layering, do_layering, map2array, get_nb_str_max, ncoef_poreice
 use soil,             only: do_soil, inertiedat
 use sorption,         only: do_sorption
@@ -493 +505 @@
 
 if (do_layering) then
-    allocate(layerings_array(ngrid,nslope,nb_str_max,6))
+    allocate(layerings_array(ngrid,nslope,nb_str_max,6 + ncoef_poreice))
     call map2array(layerings_map,layerings_array)
     call put_var_nc('stratif_top_elevation',layerings_array(:,:,:,1),itime)
@@ -500 +512 @@
     call put_var_nc('stratif_h_dust',layerings_array(:,:,:,4),itime)
     call put_var_nc('stratif_h_pore',layerings_array(:,:,:,5),itime)
-    call put_var_nc('stratif_poreice_volfrac',layerings_array(:,:,:,6),itime)
+    call put_var_nc('stratif_poreice_coef1',layerings_array(:,:,:,6),itime)
+    call put_var_nc('stratif_poreice_coef2',layerings_array(:,:,:,7),itime)
+    call put_var_nc('stratif_poreice_coef3',layerings_array(:,:,:,8),itime)
+    call put_var_nc('stratif_poreice_coef4',layerings_array(:,:,:,9),itime)
+    call put_var_nc('stratif_date',layerings_array(:,:,:,10),itime)
     deallocate(layerings_array)
 end if

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

@@ -6 +6 @@
 import os
 import sys
+from typing import Optional
 import numpy as np
 import netCDF4 as nc
@@ -44 +45 @@
 
 
+def _get_axis(dim_names: list[str], aliases: tuple[str, ...]) -> Optional[int]:
+    """Return axis index matching one alias in a list of NetCDF dimension names."""
+    for i, name in enumerate(dim_names):
+        lname = name.lower()
+        for alias in aliases:
+            if lname == alias or alias in lname:
+                return i
+    return None
+
+
+def _read_stratif_var_tsg(nc_dataset: nc.Dataset, var_name: str, slope_index: int) -> Optional[np.ndarray]:
+    """
+    Read one stratification variable and reshape it to (time, strata, grid).
+    If a slope dimension exists, keep only slope_index.
+    """
+    if var_name not in nc_dataset.variables:
+        return None
+
+    var = nc_dataset.variables[var_name]
+    arr = np.asarray(var[:], dtype=np.float64)
+    dims = list(var.dimensions)
+
+    ax_slope = _get_axis(dims, ('nslope', 'slope'))
+    if ax_slope is not None:
+        arr = np.take(arr, slope_index, axis=ax_slope)
+        dims.pop(ax_slope)
+
+    if arr.ndim == 2:
+        return arr[np.newaxis, :, :]
+    if arr.ndim != 3:
+        sys.exit(f"Error: Variable '{var_name}' has unsupported shape {arr.shape}.")
+
+    ax_time = _get_axis(dims, ('time',))
+    ax_str = _get_axis(dims, ('nb_str_max', 'nb_str'))
+    ax_grid = _get_axis(dims, ('physical_points', 'ngrid', 'grid'))
+
+    if ax_time is None or ax_str is None or ax_grid is None:
+        # Fallback to historical layout already used in this script.
+        return arr
+
+    return np.transpose(arr, (ax_time, ax_str, ax_grid))
+
+
+def _poreice_profile_from_coeff(coeff: np.ndarray, z_norm: np.ndarray) -> np.ndarray:
+    """Evaluate clamped cubic pore-ice profile on normalized depth z in [0,1]."""
+    profile = coeff[0] + coeff[1] * z_norm + coeff[2] * z_norm**2 + coeff[3] * z_norm**3
+    return np.clip(profile, 0.0, 1.0)
+
+
+def _mean_poreice_fraction(coeff: np.ndarray, n_profile: int = 33) -> float:
+    """Compute mean pore-ice volumetric fraction by integrating the cubic profile."""
+    z_norm = np.linspace(0.0, 1.0, n_profile)
+    profile = _poreice_profile_from_coeff(coeff, z_norm)
+    return float(np.clip(np.trapz(profile, z_norm), 0.0, 1.0))
+
+
 def load_slope_variables(nc_dataset: nc.Dataset, slope_index: int) -> dict:
     """
@@ -51 +108 @@
     idx_str = str(slope_index + 1).zfill(2)
     vars_base = {
-        'top_elev': f"stratif_slope{idx_str}_top_elevation",
-        'h_co2': f"stratif_slope{idx_str}_h_co2ice",
-        'h_h2o': f"stratif_slope{idx_str}_h_h2oice",
-        'h_dust': f"stratif_slope{idx_str}_h_dust",
-        'h_pore': f"stratif_slope{idx_str}_h_pore",
-        'poreice_volfrac': f"stratif_slope{idx_str}_poreice_volfrac",
+        'top_elev': (f"stratif_slope{idx_str}_top_elevation", "stratif_top_elevation"),
+        'h_co2': (f"stratif_slope{idx_str}_h_co2ice", "stratif_h_co2ice"),
+        'h_h2o': (f"stratif_slope{idx_str}_h_h2oice", "stratif_h_h2oice"),
+        'h_dust': (f"stratif_slope{idx_str}_h_dust", "stratif_h_dust"),
+        'h_pore': (f"stratif_slope{idx_str}_h_pore", "stratif_h_pore"),
     }
 
     data = {}
-    for key, var_name in vars_base.items():
-        if var_name not in nc_dataset.variables:
-            sys.exit(f"Error: Variable '{var_name}' not found in the NetCDF file.")
-        data[key] = nc_dataset.variables[var_name][:]
-    return data
+    for key, var_names in vars_base.items():
+        arr = None
+        for var_name in var_names:
+            arr = _read_stratif_var_tsg(nc_dataset, var_name, slope_index)
+            if arr is not None:
+                break
+        if arr is None:
+            sys.exit(f"Error: Variables {var_names} not found in the NetCDF file.")
+        data[key] = arr
+
+    coeff_arrays = []
+    for icoef in range(1, 5):
+        arr = _read_stratif_var_tsg(nc_dataset, f"stratif_slope{idx_str}_poreice_coef{icoef}", slope_index)
+        if arr is None:
+            arr = _read_stratif_var_tsg(nc_dataset, f"stratif_poreice_coef{icoef}", slope_index)
+        if arr is None:
+            coeff_arrays = []
+            break
+        coeff_arrays.append(arr)
+
+    if coeff_arrays:
+        data['poreice_coeff'] = np.stack(coeff_arrays, axis=-1)
+        return data
+
+    sys.exit("Error: Missing pore-ice coefficient variables (expected poreice_coef1..4).")
 
 
@@ -83 +159 @@
     h_dust = data['h_dust']
     h_pore = data['h_pore']
-    poreice_volfrac = data['poreice_volfrac']
+    poreice_coeff = data['poreice_coeff']
 
     layers = top.shape[1]
@@ -135 +211 @@
         h2o = h_h2o[0, layer_idx, grid_index] / thickness
         dust = h_dust[0, layer_idx, grid_index] / thickness
-        air = h_pore[0, layer_idx, grid_index] * (1.0 - poreice_volfrac[0, layer_idx, grid_index]) / thickness
-        poreice = h_pore[0, layer_idx, grid_index] * poreice_volfrac[0, layer_idx, grid_index] / thickness
+        poreice_frac = _mean_poreice_fraction(poreice_coeff[0, layer_idx, grid_index, :])
+        air = h_pore[0, layer_idx, grid_index] * (1.0 - poreice_frac) / thickness
+        poreice = h_pore[0, layer_idx, grid_index] * poreice_frac / thickness
         contents[:, idx + 1] = [co2, h2o, dust, air, poreice]
 

trunk/LMDZ.COMMON/libf/evolution/deftank/visu_layering_evo.py

@@ -97 +97 @@
 
 
+def _get_axis(dim_names, aliases):
+    """Return axis index matching aliases in a dimension-name list."""
+    for i, name in enumerate(dim_names):
+        lname = name.lower()
+        for alias in aliases:
+            if lname == alias or alias in lname:
+                return i
+    return None
+
+
+def _var_to_tsgi(nc_var):
+    """
+    Convert a NetCDF stratification variable to shape (time, slope, strata, grid).
+    Handles both slope-specific variables and generic variables with an nslope dimension.
+    """
+    arr = np.asarray(nc_var[:], dtype=np.float32)
+    dims = [d.lower() for d in nc_var.dimensions]
+
+    ax_time = _get_axis(dims, ('time',))
+    ax_slope = _get_axis(dims, ('nslope', 'slope'))
+    ax_str = _get_axis(dims, ('nb_str_max', 'nb_str'))
+    ax_grid = _get_axis(dims, ('physical_points', 'ngrid', 'grid'))
+
+    if arr.ndim == 4 and None not in (ax_time, ax_slope, ax_str, ax_grid):
+        return np.transpose(arr, (ax_time, ax_slope, ax_str, ax_grid))
+
+    if arr.ndim == 3 and None not in (ax_time, ax_str, ax_grid):
+        arr_tsg = np.transpose(arr, (ax_time, ax_str, ax_grid))
+        return arr_tsg[:, np.newaxis, :, :]
+
+    if arr.ndim == 3 and None not in (ax_slope, ax_str, ax_grid):
+        arr_ssg = np.transpose(arr, (ax_slope, ax_str, ax_grid))
+        return arr_ssg[np.newaxis, :, :, :]
+
+    if arr.ndim == 3:
+        # Historical fallback used by this script: (time, strata, grid)
+        return arr[:, np.newaxis, :, :]
+
+    if arr.ndim == 2:
+        return arr[np.newaxis, np.newaxis, :, :]
+
+    raise ValueError(f"Unsupported variable shape {arr.shape} for '{nc_var.name}'")
+
+
+def _slope_index_from_var(vname):
+    """Return 0-based slope index from a slope-specific variable name, else None."""
+    if 'slope' not in vname:
+        return None
+    try:
+        return int(vname.split('slope')[1].split('_')[0]) - 1
+    except (ValueError, IndexError):
+        return None
+
+
 def collect_stratification_variables(files, base_name):
     """
@@ -114 +168 @@
 
     property_markers = {
-        'heights':   'stratif_slope',    # "..._top_elevation"
+        'heights':   'top_elevation',
         'co2_ice':   'h_co2ice',
         'h2o_ice':   'h_h2oice',
         'dust':      'h_dust',
         'pore':      'h_pore',
-        'pore_ice':  'poreice_volfrac'
+        'poreice_coef1': 'poreice_coef1',
+        'poreice_coef2': 'poreice_coef2',
+        'poreice_coef3': 'poreice_coef3',
+        'poreice_coef4': 'poreice_coef4'
     }
     var_info = {prop: set() for prop in property_markers}
@@ -128 +185 @@
                 max_nb_str = max(max_nb_str, ds.dimensions['nb_str_max'].size)
 
-            nslope = ds.dimensions['nslope'].size
-            for k in range(1, nslope + 1):
-                var_name = f"stratif_slope{k:02d}_top_elevation"
-                if var_name in ds.variables:
-                    arr = ds.variables[var_name][:]
-                    min_base_elev = min(min_base_elev, np.min(arr))
-                    max_top_elev = max(max_top_elev, np.max(arr))
-                    var_info['heights'].add(var_name)
-
             for full_var in ds.variables:
+                if 'stratif' not in full_var:
+                    continue
+
+                if property_markers['heights'] in full_var:
+                    arr = np.asarray(ds.variables[full_var][:])
+                    if arr.size > 0:
+                        min_base_elev = min(min_base_elev, float(np.nanmin(arr)))
+                        max_top_elev = max(max_top_elev, float(np.nanmax(arr)))
+                    var_info['heights'].add(full_var)
+
                 for prop, marker in property_markers.items():
-                    if (marker in full_var) and prop != 'heights':
+                    if prop != 'heights' and marker in full_var:
                         var_info[prop].add(full_var)
 
@@ -174 +232 @@
     for t_idx, ds in enumerate(datasets):
         for var_name in var_info['heights']:
-            slope_idx = int(var_name.split("slope")[1].split("_")[0]) - 1
-            if 0 <= slope_idx < nslope:
-                raw = ds.variables[var_name][0, :, :]  # (n_strata, ngrid)
-                heights_data[t_idx][slope_idx] = raw.T  # (ngrid, n_strata)
+            if var_name not in ds.variables:
+                continue
+            raw_tsgi = _var_to_tsgi(ds.variables[var_name])
+            raw_sgi = raw_tsgi[0]  # (nslope_var, nstrata, ngrid)
+
+            slope_idx = _slope_index_from_var(var_name)
+            if slope_idx is not None and raw_sgi.shape[0] == 1:
+                if 0 <= slope_idx < nslope:
+                    heights_data[t_idx][slope_idx] = raw_sgi[0].T
+                continue
+
+            nsl = min(nslope, raw_sgi.shape[0])
+            for isl in range(nsl):
+                heights_data[t_idx][isl] = raw_sgi[isl].T
 
     raw_prop_arrays = {}
@@ -185 +253 @@
         raw_prop_arrays[prop] = np.zeros((ngrid, ntime, nslope, max_nb_str), dtype=np.float32)
 
-    def slope_index_from_var(vname):
-        return int(vname.split("slope")[1].split("_")[0]) - 1
-
     for prop in raw_prop_arrays:
-        slope_map = {}
-        for vname in var_info[prop]:
-            isl = slope_index_from_var(vname)
-            if 0 <= isl < nslope:
-                slope_map[isl] = vname
-
         arr = raw_prop_arrays[prop]
         for t_idx, ds in enumerate(datasets):
-            for isl, var_name in slope_map.items():
-                raw = ds.variables[var_name][0, :, :]  # (n_strata, ngrid)
-                n_strata_current = raw.shape[0]
-                arr[:, t_idx, isl, :n_strata_current] = raw.T
+            for var_name in var_info[prop]:
+                if var_name not in ds.variables:
+                    continue
+
+                raw_tsgi = _var_to_tsgi(ds.variables[var_name])
+                raw_sgi = raw_tsgi[0]  # (nslope_var, nstrata, ngrid)
+                slope_idx = _slope_index_from_var(var_name)
+
+                if slope_idx is not None and raw_sgi.shape[0] == 1:
+                    if 0 <= slope_idx < nslope:
+                        n_strata_current = min(raw_sgi.shape[1], max_nb_str)
+                        arr[:, t_idx, slope_idx, :n_strata_current] = raw_sgi[0, :n_strata_current, :].T
+                    continue
+
+                nsl = min(nslope, raw_sgi.shape[0])
+                for isl in range(nsl):
+                    n_strata_current = min(raw_sgi.shape[1], max_nb_str)
+                    arr[:, t_idx, isl, :n_strata_current] = raw_sgi[isl, :n_strata_current, :].T
 
     return heights_data, raw_prop_arrays, ntime
@@ -209 +282 @@
     Given raw_prop_arrays for 'co2_ice', 'h2o_ice', 'dust', 'pore' (in meters),
     normalize each set of strata so that the sum of those four = 1 per cell.
-    Returns:
-      - frac_arrays: dict mapping same keys -> 4D arrays of fractions (0..1).
+        Returns:
+            - frac_arrays: dict mapping same keys -> 4D arrays of fractions (0..1).
+            - includes an additional 'pore_ice' fraction derived from pore-ice coefficients.
     """
     co2 = raw_prop_arrays['co2_ice']
@@ -224 +298 @@
     frac_dust = np.zeros_like(dust, dtype=np.float32)
     frac_pore = np.zeros_like(pore, dtype=np.float32)
+    frac_pore_ice = np.zeros_like(pore, dtype=np.float32)
 
     frac_co2[mask] = co2[mask] / total[mask]
@@ -230 +305 @@
     frac_pore[mask] = pore[mask] / total[mask]
 
+    if all(key in raw_prop_arrays for key in ('poreice_coef1', 'poreice_coef2', 'poreice_coef3', 'poreice_coef4')):
+        z_norm = np.linspace(0.0, 1.0, 33, dtype=np.float32)
+        c1 = raw_prop_arrays['poreice_coef1']
+        c2 = raw_prop_arrays['poreice_coef2']
+        c3 = raw_prop_arrays['poreice_coef3']
+        c4 = raw_prop_arrays['poreice_coef4']
+        profile = (c1[..., None] + c2[..., None] * z_norm + c3[..., None] * z_norm**2 + c4[..., None] * z_norm**3)
+        profile = np.clip(profile, 0.0, 1.0)
+        poreice_mean = np.clip(np.trapz(profile, z_norm, axis=-1), 0.0, 1.0).astype(np.float32)
+    else:
+        poreice_mean = np.zeros_like(pore, dtype=np.float32)
+
+    frac_pore_ice[mask] = pore[mask] * poreice_mean[mask] / total[mask]
+
     return {
         'co2_ice': frac_co2,
         'h2o_ice': frac_h2o,
         'dust':     frac_dust,
-        'pore':     frac_pore
+        'pore':     frac_pore,
+        'pore_ice': frac_pore_ice
     }
 

trunk/LMDZ.COMMON/libf/evolution/display.F90

@@ -35 +35 @@
 character(32),  protected, private :: username = ' '    ! User name
 character(32),  protected, private :: hostname = ' '    ! Machine/station name
+character(128), protected, private :: cmd_pgrm = ' '    ! Command used to run the program
+character(8),   protected, private :: date = ' '        ! Current date (YYYYMMDD)
+character(10),  protected, private :: time = ' '        ! Current time (hhmmss.sss)
+character(5),   protected, private :: zone = ' '        ! UTC offset (+/-hhmm)
 integer(di),    protected, private :: verbosity_lvl = LVL_NFO
 logical(k4),    protected, private :: out2term = .true. ! Flag to output to terminal
@@ -125 +129 @@
 ! LOCAL VARIABLES
 ! ---------------
-character(8)              :: date
-character(10)             :: time
-character(5)              :: zone
 character(100)            :: msg
-integer(di), dimension(8) :: values
+integer(di), dimension(8) :: values ! Date/time components
 
 ! CODE
@@ -147 +148 @@
 call getlog(username)
 call hostnm(hostname)
+call get_command(cmd_pgrm)
 call print_msg('',LVL_NFO)
 call print_msg('********* PEM information *********',LVL_NFO)
-call print_msg('+ User     : '//trim(username),LVL_NFO)
-call print_msg('+ Machine  : '//trim(hostname),LVL_NFO)
-call print_msg('+ Directory: '//trim(curr_dir),LVL_NFO)
-write(msg,'(a,i2,a,i2,a,i4)') '+ Date     : ',values(3),'/',values(2),'/',values(1)
-call print_msg(msg,LVL_NFO)
-write(msg,'(a,i2,a,i2,a,i2,a)') '+ Time     : ',values(5),':',values(6),':',values(7)
-call print_msg(msg,LVL_NFO)
+call print_msg('+ User     : '//trim(adjustl(username)),LVL_NFO)
+call print_msg('+ Machine  : '//trim(adjustl(hostname)),LVL_NFO)
+call print_msg('+ Directory: '//trim(adjustl(curr_dir)),LVL_NFO)
+call print_msg('+ Command  : '//trim(adjustl(cmd_pgrm)),LVL_NFO)
+write(msg,'(a,a4,a,a2,a,a2)') '+ Date     : ',date(1:4),'-',date(5:6),'-',date(7:8)
+call print_msg(trim(msg),LVL_NFO)
+write(msg,'(a,a2,a,a2,a,a2,a,a3)') '+ Time     : ',time(1:2),':',time(3:4),':',time(5:6),'.',time(8:10)
+call print_msg(trim(msg),LVL_NFO)
+if (zone(1:1) == '+' .or. zone(1:1) == '-') then
+    write(msg,'(a,a1,a2,a,a2)') '+ Timezone : UTC',zone(1:1),zone(2:3),':',zone(4:5)
+else
+    write(msg,'(a,i0,a)') '+ Timezone : unknown (offset=',values(4),' min)'
+end if
+call print_msg(trim(msg),LVL_NFO)
 call print_msg('',LVL_NFO)
 call print_msg('********* Initialization *********',LVL_NFO)

trunk/LMDZ.COMMON/libf/evolution/layered_deposits.F90

@@ -5 +5 @@
 !
 ! DESCRIPTION
-!     Manage layered deposits in the PEM with a linked list data structure.
+!     Manages layered deposits thanks to a linked list data structure.
 !     Provides data types and subroutines to manipulate layers of ice and dust.
 !
 ! AUTHORS & DATE
-!     JB Clement, 2024
+!     JB Clement, 03/2024
 !
 ! NOTES
@@ -18 +18 @@
 ! ------------
 use numerics, only: dp, di, k4, tol, eps
-use geometry, only: ngrid, nslope
 use surf_ice, only: rho_co2ice, rho_h2oice
 
@@ -31 +30 @@
 
 ! Physical parameters
-logical(k4), protected, private :: impose_dust_ratio            ! Impose dust-to-ice ratio instead of dust tendency (see 'dust2ice_ratio')
-real(dp),    protected, private :: dust2ice_ratio               ! Dust-to-ice ratio when ice condenses (10% by default)
-real(dp),    protected, private :: d_dust                       ! Tendency of dust [kg/m2/y]
-real(dp),    parameter, private :: dry_lag_porosity  = 0.2_dp   ! Porosity of dust lag
-real(dp),    parameter, private :: wet_lag_porosity  = 0.15_dp  ! Porosity of water ice lag
-real(dp),    parameter, private :: co2ice_porosity   = 0.05_dp  ! Porosity of CO2 ice
-real(dp),    parameter, private :: h2oice_porosity   = 0.1_dp   ! Porosity of H2O ice
-real(dp),    parameter, private :: rho_dust          = 2500._dp ! Density of dust [kg.m-3]
-real(dp),    parameter, private :: h_patchy_dust     = 5.e-4_dp ! Height under which a dust lag is considered patchy [m]
-real(dp),    parameter, private :: h_patchy_ice      = 5.e-4_dp ! Height under which a H2O ice lag is considered patchy [m]
+logical(k4), protected, private :: impose_dust_ratio           ! Impose dust-to-ice ratio instead of dust tendency (see 'dust2ice_ratio')
+real(dp),    protected, private :: dust2ice_ratio              ! Dust-to-ice ratio when ice condenses (10% by default)
+real(dp),    protected, private :: d_dust                      ! Tendency of dust [kg/m2/y]
+real(dp),    parameter, private :: dry_lag_porosity = 0.2_dp   ! Porosity of dust lag
+real(dp),    parameter, private :: wet_lag_porosity = 0.15_dp  ! Porosity of water ice lag
+real(dp),    parameter, private :: co2ice_porosity  = 0.05_dp  ! Porosity of CO2 ice
+real(dp),    parameter, private :: h2oice_porosity  = 0.1_dp   ! Porosity of H2O ice
+real(dp),    parameter, private :: rho_dust         = 2500._dp ! Density of dust [kg.m-3]
+real(dp),    parameter, private :: h_patchy_dust    = 5.e-4_dp ! Height under which a dust lag is considered patchy [m]
+real(dp),    parameter, private :: h_patchy_ice     = 5.e-4_dp ! Height under which a H2O ice lag is considered patchy [m]
+
+! Stratum type codes
+integer(di), parameter, private :: STR_TYPE_EMPTY  = 0_di ! Empty stratum
+integer(di), parameter, private :: STR_TYPE_CO2ICE = 1_di ! CO2 ice stratum
+integer(di), parameter, private :: STR_TYPE_H2OICE = 2_di ! H2O ice stratum
+integer(di), parameter, private :: STR_TYPE_DUST   = 3_di ! Dust stratum
+integer(di), parameter, private :: STR_TYPE_MIXED  = 4_di ! Mixed stratum (ice and dust)
+
+! Parameters for pore ice profile reconstruction
+integer(di), parameter          :: ncoef_poreice = 4_di    ! Number of cubic coefficients for pore ice profile
+real(dp),    parameter, private :: dz_poreice    = 0.01_dp ! Target depth step for pore-ice profile reconstruction [m]
+integer(di), parameter, private :: nmin_poreice  = 4_di    ! Minimal number of points for pore-ice profile reconstruction
+integer(di), parameter, private :: nmax_poreice  = 64_di   ! Maximal number of points for pore-ice profile reconstruction
 
 ! Parameters for CO2 flux correction (inspired by Levrard et al. 2007)
@@ -50 +62 @@
 ! Stratum = node of the linked list
 type :: stratum
-    real(dp)               :: top_elevation   ! Layer top_elevation (top height from the surface) [m]
-    real(dp)               :: h_co2ice        ! Height of CO2 ice [m]
-    real(dp)               :: h_h2oice        ! Height of H2O ice [m]
-    real(dp)               :: h_dust          ! Height of dust [m]
-    real(dp)               :: h_pore          ! Height of pores [m]
-    real(dp)               :: poreice_volfrac ! Pore-filled ice volumetric fraction
-    type(stratum), pointer :: up => null()    ! Upper stratum (next node)
-    type(stratum), pointer :: down => null()  ! Lower stratum (previous node)
+    real(dp)                           :: top_elevation = 1._dp ! Layer top_elevation (top height from the surface) [m]
+    real(dp)                           :: h_co2ice      = 0._dp ! Height of CO2 ice [m]
+    real(dp)                           :: h_h2oice      = 0._dp ! Height of H2O ice [m]
+    real(dp)                           :: h_dust        = 1._dp ! Height of dust [m]
+    real(dp)                           :: h_pore        = 0._dp ! Height of pores [m]
+    real(dp), dimension(ncoef_poreice) :: coef_poreice  = 0._dp ! Coefficients for pore ice profile
+    real(dp)                           :: form_date     = 0._dp ! Stratum formation date [planetary years]
+    type(stratum), pointer             :: up   => null()        ! Upper stratum (next node)
+    type(stratum), pointer             :: down => null()        ! Lower stratum (previous node)
 end type stratum
 
 ! Layering = linked list
-type layering
+type :: layering
     integer(di)            :: nb_str = 0       ! Number of strata
-    type(stratum), pointer :: bottom => null() ! First stratum (head of the list)
-    type(stratum), pointer :: top => null()    ! Last stratum (tail of the list)
-    !real(dp)               :: h_toplag         ! Height of dust layer at the top [m]
-    !type(stratum), pointer :: subice => null() ! Shallower stratum containing ice
-    !type(stratum), pointer :: surf => null()   ! Surface stratum
+    type(stratum), pointer :: bottom => null() ! First stratum (head of the list) at the very bottom
+    type(stratum), pointer :: top    => null() ! Last stratum (tail of the list) at the very top
+    !type(stratum), pointer :: subice => null() ! Shallowest stratum containing compact ice
 end type layering
 
 ! Array of pointers to "stratum"
-type ptrarray
+type :: ptrarray
     type(stratum), pointer :: p
 end type ptrarray
@@ -125 +136 @@
 
 !=======================================================================
+FUNCTION get_str_type(str) RESULT(type)
+!-----------------------------------------------------------------------
+! NAME
+!     get_str_type
+!
+! DESCRIPTION
+!     Determine stratum type from composition and pore-ice state.
+!
+! AUTHORS & DATE
+!     JB Clement, 04/2026
+!
+! NOTES
+!
+!-----------------------------------------------------------------------
+
+! DECLARATION
+! -----------
+implicit none
+
+! ARGUMENTS
+! ---------
+type(stratum), intent(in) :: str
+
+! LOCAL VARIABLES
+! ---------------
+integer(di) :: type
+
+! CODE
+! ----
+if (str%h_co2ice <= tol .and. str%h_h2oice <= tol .and. str%h_dust <= tol .and. is_poreice_empty(str)) then
+    type = STR_TYPE_EMPTY
+else if (str%h_co2ice > str%h_h2oice .and. str%h_co2ice > str%h_dust) then
+    type = STR_TYPE_CO2ICE
+else if (str%h_h2oice > str%h_co2ice .and. str%h_h2oice > str%h_dust) then
+    type = STR_TYPE_H2OICE
+else if (str%h_dust > str%h_co2ice .and. str%h_dust > str%h_h2oice) then
+    type = STR_TYPE_DUST
+else
+    type = STR_TYPE_MIXED
+end if
+
+END FUNCTION get_str_type
+!=======================================================================
+
+!=======================================================================
+FUNCTION strtype2txt(type) RESULT(txt_type)
+!-----------------------------------------------------------------------
+! NAME
+!     strtype2txt
+!
+! DESCRIPTION
+!     Convert stratum type integer code into a readable string.
+!
+! AUTHORS & DATE
+!     JB Clement, 04/2026
+!
+! NOTES
+!
+!-----------------------------------------------------------------------
+
+! DECLARATION
+! -----------
+implicit none
+
+! ARGUMENTS
+! ---------
+integer(di), intent(in) :: type
+
+! LOCAL VARIABLES
+! ---------------
+character(16) :: txt_type
+
+! CODE
+! ----
+select case (type)
+    case (STR_TYPE_EMPTY)
+        txt_type = 'EMPTY'
+    case (STR_TYPE_CO2ICE)
+        txt_type = 'CO2_ICE'
+    case (STR_TYPE_H2OICE)
+        txt_type = 'H2O_ICE'
+    case (STR_TYPE_DUST)
+        txt_type = 'DUST'
+    case (STR_TYPE_MIXED)
+        txt_type = 'MIXED'
+    case default
+        txt_type = 'UNKNOWN'
+end select
+
+END FUNCTION strtype2txt
+!=======================================================================
+
+!=======================================================================
 SUBROUTINE print_stratum(this)
 !-----------------------------------------------------------------------
@@ -153 +257 @@
 type(stratum), intent(in) :: this
 
-! CODE
-! ----
-call print_msg('Top elevation       [m]: '//real2str(this%top_elevation),LVL_DBG)
-call print_msg('Height of CO2 ice   [m]: '//real2str(this%h_co2ice),LVL_DBG)
-call print_msg('Height of H2O ice   [m]: '//real2str(this%h_h2oice),LVL_DBG)
-call print_msg('Height of dust      [m]: '//real2str(this%h_dust),LVL_DBG)
-call print_msg('Height of pores     [m]: '//real2str(this%h_pore),LVL_DBG)
-call print_msg('Pore-filled ice [m3/m3]: '//real2str(this%poreice_volfrac),LVL_DBG)
+! LOCAL VARIABLES
+! ---------------
+integer(di)    :: i
+character(128) :: msg
+
+! CODE
+! ----
+call print_msg('Nature                     [-]: '//trim(strtype2txt(get_str_type(this))),LVL_DBG)
+call print_msg('Formation date     [plnt year]: '//real2str(this%form_date),LVL_DBG)
+call print_msg('Top elevation              [m]: '//real2str(this%top_elevation),LVL_DBG)
+call print_msg('Height of CO2 ice          [m]: '//real2str(this%h_co2ice),LVL_DBG)
+call print_msg('Height of H2O ice          [m]: '//real2str(this%h_h2oice),LVL_DBG)
+call print_msg('Height of dust             [m]: '//real2str(this%h_dust),LVL_DBG)
+call print_msg('Height of pores            [m]: '//real2str(this%h_pore),LVL_DBG)
+msg = 'Pore-ice profile coefs [m3/m3]: {'
+do i = 1, size(this%coef_poreice)
+    msg = trim(msg)//real2str(this%coef_poreice(i))//', '
+end do
+msg = trim(msg)//'}'
+call print_msg(trim(msg),LVL_DBG)
 
 END SUBROUTINE print_stratum
@@ -166 +282 @@
 
 !=======================================================================
-SUBROUTINE add_stratum(this,top_elevation,co2ice,h2oice,dust,pore,poreice)
+SUBROUTINE add_stratum(this,top_elevation,co2ice,h2oice,dust,pore,coef_poreice,form_date)
 !-----------------------------------------------------------------------
 ! NAME
@@ -187 +303 @@
 ! ARGUMENTS
 ! ---------
-type(layering), intent(inout)        :: this
-real(dp),       intent(in), optional :: top_elevation, co2ice, h2oice, dust, pore, poreice
+type(layering),               intent(inout)        :: this
+real(dp),                     intent(in), optional :: top_elevation, co2ice, h2oice, dust, pore
+real(dp),       dimension(:), intent(in), optional :: coef_poreice
+real(dp),                     intent(in), optional :: form_date
 
 ! LOCAL VARIABLES
@@ -198 +316 @@
 ! Creation of the stratum
 allocate(str)
-nullify(str%up,str%down)
-str%top_elevation   = 1._dp
-str%h_co2ice        = 0._dp
-str%h_h2oice        = 0._dp
-str%h_dust          = 1._dp
-str%h_pore          = 0._dp
-str%poreice_volfrac = 0._dp
 if (present(top_elevation)) str%top_elevation   = top_elevation
 if (present(co2ice))        str%h_co2ice        = co2ice
@@ -210 +321 @@
 if (present(dust))          str%h_dust          = dust
 if (present(pore))          str%h_pore          = pore
-if (present(poreice))       str%poreice_volfrac = poreice
+if (present(coef_poreice))  str%coef_poreice(:) = coef_poreice(:)
+if (present(form_date))     str%form_date       = form_date
 
 ! Increment the number of strata
@@ -229 +341 @@
 
 !=======================================================================
-SUBROUTINE insert_stratum(this,str_prev,top_elevation,co2ice,h2oice,dust,pore,poreice)
+SUBROUTINE insert_stratum(this,str_prev,top_elevation,co2ice,h2oice,dust,pore,coef_poreice,form_date)
 !-----------------------------------------------------------------------
 ! NAME
@@ -252 +364 @@
 type(layering),         intent(inout)        ::  this
 type(stratum), pointer, intent(inout)        ::  str_prev
-real(dp),               intent(in), optional :: top_elevation, co2ice, h2oice, dust, pore, poreice
+real(dp),               intent(in), optional :: top_elevation, co2ice, h2oice, dust, pore
+real(dp), dimension(:), intent(in), optional :: coef_poreice
+real(dp),               intent(in), optional :: form_date
 
 ! LOCAL VARIABLES
@@ -261 +375 @@
 ! ----
 if (.not. associated(str_prev%up)) then ! If str_prev is the top stratum of the layering
-    call add_stratum(this,top_elevation,co2ice,h2oice,dust,pore,poreice)
+    call add_stratum(this,top_elevation,co2ice,h2oice,dust,pore,coef_poreice,form_date)
 else
     ! Creation of the stratum
     allocate(str)
-    nullify(str%up,str%down)
-    str%top_elevation   = 1._dp
-    str%h_co2ice        = 0._dp
-    str%h_h2oice        = 0._dp
-    str%h_dust          = 1._dp
-    str%h_pore          = 0._dp
-    str%poreice_volfrac = 0._dp
     if (present(top_elevation)) str%top_elevation   = top_elevation
     if (present(co2ice))        str%h_co2ice        = co2ice
@@ -277 +384 @@
     if (present(dust))          str%h_dust          = dust
     if (present(pore))          str%h_pore          = pore
-    if (present(poreice))       str%poreice_volfrac = poreice
+    if (present(coef_poreice))  str%coef_poreice(:) = coef_poreice(:)
+    if (present(form_date))     str%form_date       = form_date
 
     ! Increment the number of strata
@@ -359 +467 @@
 
 !=======================================================================
+FUNCTION are_strata_similar(str_lo,str_hi,frac_tol,coef_tol,date_tol) RESULT(similar)
+!-----------------------------------------------------------------------
+! NAME
+!     are_strata_similar
+!
+! DESCRIPTION
+!     Check if two adjacent strata can be merged based on composition,
+!     pore-ice profile, and optional formation-date proximity.
+!
+! AUTHORS & DATE
+!     JB Clement, 04/2026
+!
+! NOTES
+!     Negative value for date_tol means that formation date is not checked
+!     for similarity.
+!-----------------------------------------------------------------------
+
+! DECLARATION
+! -----------
+implicit none
+
+! ARGUMENTS
+! ---------
+type(stratum), intent(in) :: str_lo, str_hi
+real(dp),      intent(in) :: frac_tol, coef_tol, date_tol
+
+! LOCAL VARIABLES
+! ---------------
+logical(k4)         :: similar
+real(dp)            :: h_lo, h_hi
+real(dp)            :: f_lo_co2, f_lo_h2o, f_lo_dust, f_lo_pore
+real(dp)            :: f_hi_co2, f_hi_h2o, f_hi_dust, f_hi_pore
+
+! CODE
+! ----
+similar = .false.
+
+! Keep the neutral base layer untouched.
+if (str_lo%top_elevation <= tol) return
+
+! Strata must be of the same type to be merged
+if (get_str_type(str_lo) /= get_str_type(str_hi)) return
+
+! Check that the composition fractions are similar within the tolerance
+h_lo = thickness(str_lo)
+f_lo_co2  = str_lo%h_co2ice/h_lo
+f_lo_h2o  = str_lo%h_h2oice/h_lo
+f_lo_dust = str_lo%h_dust/h_lo
+f_lo_pore = str_lo%h_pore/h_lo
+h_hi = thickness(str_hi)
+f_hi_co2   = str_hi%h_co2ice/h_hi
+f_hi_h2o   = str_hi%h_h2oice/h_hi
+f_hi_dust  = str_hi%h_dust/h_hi
+f_hi_pore  = str_hi%h_pore/h_hi
+if (abs(f_lo_co2 - f_hi_co2) > frac_tol) return
+if (abs(f_lo_h2o - f_hi_h2o) > frac_tol) return
+if (abs(f_lo_dust - f_hi_dust) > frac_tol) return
+if (abs(f_lo_pore - f_hi_pore) > frac_tol) return
+if (maxval(abs(str_lo%coef_poreice - str_hi%coef_poreice)) > coef_tol) return
+
+! Optionally check that the formation dates are close enough
+if (date_tol >= 0._dp) then
+    if (abs(str_lo%form_date - str_hi%form_date) > date_tol) return
+end if
+
+similar = .true.
+
+END FUNCTION are_strata_similar
+!=======================================================================
+
+!=======================================================================
+SUBROUTINE merge_adjacent_strata(this,frac_tol,coef_tol,date_tol,max_pass)
+!-----------------------------------------------------------------------
+! NAME
+!     merge_adjacent_strata
+!
+! DESCRIPTION
+!     Merge adjacent strata with similar properties to limit unnecessary
+!     stratification growth.
+!
+! AUTHORS & DATE
+!     JB Clement, 04/2026
+!
+! NOTES
+!     Negative value for date_tol means that formation date is not checked
+!     for similarity.
+!-----------------------------------------------------------------------
+
+! DECLARATION
+! -----------
+implicit none
+
+! ARGUMENTS
+! ---------
+type(layering), intent(inout)        :: this
+real(dp),       intent(in), optional :: frac_tol, coef_tol, date_tol
+integer(di),    intent(in), optional :: max_pass
+
+! LOCAL VARIABLES
+! ---------------
+type(stratum), pointer :: str_lo, str_hi
+real(dp)               :: h_lo, h_hi, h_tot, w_lo, w_hi
+real(dp)               :: frac_tol_loc, coef_tol_loc, date_tol_loc
+logical(k4)            :: merged_in_pass
+integer(di)            :: ipass, max_pass_loc
+
+! CODE
+! ----
+if (this%nb_str <= 2) return
+
+frac_tol_loc = 1.e-2_dp
+coef_tol_loc = 1.e-2_dp
+date_tol_loc = -1._dp
+if (present(frac_tol)) frac_tol_loc = max(0._dp,frac_tol)
+if (present(coef_tol)) coef_tol_loc = max(0._dp,coef_tol)
+if (present(date_tol)) date_tol_loc = date_tol
+
+max_pass_loc = this%nb_str
+if (present(max_pass)) max_pass_loc = max(1_di,max_pass)
+
+do ipass = 1,max_pass_loc
+    merged_in_pass = .false.
+    str_lo => this%bottom
+
+    do while (associated(str_lo) .and. associated(str_lo%up))
+        str_hi => str_lo%up
+
+        if (are_strata_similar(str_lo,str_hi,frac_tol_loc,coef_tol_loc,date_tol_loc)) then
+            h_lo = thickness(str_lo)
+            h_hi = thickness(str_hi)
+            h_tot = h_lo + h_hi
+            w_lo = h_lo/h_tot
+            w_hi = h_hi/h_tot
+
+            str_lo%top_elevation = str_hi%top_elevation
+            str_lo%h_co2ice = str_lo%h_co2ice + str_hi%h_co2ice
+            str_lo%h_h2oice = str_lo%h_h2oice + str_hi%h_h2oice
+            str_lo%h_dust = str_lo%h_dust + str_hi%h_dust
+            str_lo%h_pore = str_lo%h_pore + str_hi%h_pore
+            str_lo%coef_poreice(:) = w_lo*str_lo%coef_poreice(:) + w_hi*str_hi%coef_poreice(:)
+            str_lo%form_date = w_lo*str_lo%form_date + w_hi*str_hi%form_date
+
+            call del_stratum(this,str_hi)
+            merged_in_pass = .true.
+            cycle
+        end if
+
+        str_lo => str_lo%up
+    end do
+
+    if (.not. merged_in_pass) exit
+end do
+
+END SUBROUTINE merge_adjacent_strata
+!=======================================================================
+
+!=======================================================================
 FUNCTION get_stratum(lay,i) RESULT(str)
 !-----------------------------------------------------------------------
@@ -409 +674 @@
 !-----------------------------------------------------------------------
 ! NAME
-!     del_layering
+!     ini_layering
 !
 ! DESCRIPTION
@@ -421 +686 @@
 !-----------------------------------------------------------------------
 
-! DEPENDENCIES
-! ------------
-use soil,     only: do_soil, layer, index_breccia, index_bedrock, regolith_porosity, breccia_porosity
-use geometry, only: nsoil
+!~ ! DEPENDENCIES
+!~ ! ------------
+!~ use soil,     only: do_soil, layer, index_breccia, index_bedrock, regolith_porosity, breccia_porosity
+!~ use geometry, only: nsoil
 
 ! DECLARATION
@@ -436 +701 @@
 ! LOCAL VARIABLES
 ! ---------------
-integer(di) :: i
-real(dp)    :: h_soil, h_pore, h_tot
-
-! CODE
-! ----
-! Creation of strata at the bottom of the layering to describe the sub-surface
-if (do_soil) then
-    do i = nsoil,index_bedrock,-1
-        h_soil = layer(i) - layer(i - 1) ! No porosity
-        call add_stratum(this,-layer(i - 1),0._dp,0._dp,h_soil,0._dp,0._dp)
-    end do
-    do i = index_bedrock - 1,index_breccia,-1
-        h_tot = layer(i) - layer(i - 1)
-        h_pore = h_tot*breccia_porosity
-        h_soil = h_tot - h_pore
-        call add_stratum(this,-layer(i - 1),0._dp,0._dp,h_soil,h_pore,0._dp)
-    end do
-    do i = index_breccia - 1,2,-1
-        h_tot = layer(i) - layer(i - 1)
-        h_pore = h_tot*regolith_porosity
-        h_soil = h_tot - h_pore
-        call add_stratum(this,-layer(i - 1),0._dp,0._dp,h_soil,h_pore,0._dp)
-    end do
-    h_pore = layer(1)*regolith_porosity
-    h_soil = layer(1) - h_pore
-    call add_stratum(this,0._dp,0._dp,0._dp,h_soil,h_pore,0._dp)
-end if
+real(dp), dimension(ncoef_poreice), parameter :: coef_poreice_0 = 0._dp
+!~ integer(di) :: i
+!~ real(dp)    :: h_soil, h_pore, h_tot
+
+! CODE
+! ----
+!~ ! Creation of strata at the bottom of the layering to describe the sub-surface
+!~ if (do_soil) then
+!~     do i = nsoil,index_bedrock,-1
+!~         h_soil = layer(i) - layer(i - 1) ! No porosity
+!~         call add_stratum(this,-layer(i - 1),0._dp,0._dp,h_soil,0._dp,coef_poreice_0)
+!~     end do
+!~     do i = index_bedrock - 1,index_breccia,-1
+!~         h_tot = layer(i) - layer(i - 1)
+!~         h_pore = h_tot*breccia_porosity
+!~         h_soil = h_tot - h_pore
+!~         call add_stratum(this,-layer(i - 1),0._dp,0._dp,h_soil,h_pore,coef_poreice_0)
+!~     end do
+!~     do i = index_breccia - 1,2,-1
+!~         h_tot = layer(i) - layer(i - 1)
+!~         h_pore = h_tot*regolith_porosity
+!~         h_soil = h_tot - h_pore
+!~         call add_stratum(this,-layer(i - 1),0._dp,0._dp,h_soil,h_pore,coef_poreice_0)
+!~     end do
+!~     h_pore = layer(1)*regolith_porosity
+!~     h_soil = layer(1) - h_pore
+!~     call add_stratum(this,0._dp,0._dp,0._dp,h_soil,h_pore,coef_poreice_0)
+!~ end if
+
+! Creation of a neutral base stratum at the surface
+call add_stratum(this,0._dp,0._dp,0._dp,0._dp,0._dp,coef_poreice_0)
 
 END SUBROUTINE ini_layering
@@ -648 +917 @@
         k = 1
         do while (associated(current))
-            layerings_array(ig,islope,k,1) = current%top_elevation
-            layerings_array(ig,islope,k,2) = current%h_co2ice
-            layerings_array(ig,islope,k,3) = current%h_h2oice
-            layerings_array(ig,islope,k,4) = current%h_dust
-            layerings_array(ig,islope,k,5) = current%h_pore
-            layerings_array(ig,islope,k,6) = current%poreice_volfrac
+            layerings_array(ig,islope,k,1)   = current%top_elevation
+            layerings_array(ig,islope,k,2)   = current%h_co2ice
+            layerings_array(ig,islope,k,3)   = current%h_h2oice
+            layerings_array(ig,islope,k,4)   = current%h_dust
+            layerings_array(ig,islope,k,5)   = current%h_pore
+            layerings_array(ig,islope,k,6:9) = current%coef_poreice(:)
+            layerings_array(ig,islope,k,10)  = current%form_date
             current => current%up
             k = k + 1
@@ -704 +974 @@
                              layerings_array(ig,islope,k,2),layerings_array(ig,islope,k,3), &
                              layerings_array(ig,islope,k,4),layerings_array(ig,islope,k,5), &
-                             layerings_array(ig,islope,k,6))
+                             layerings_array(ig,islope,k,6:9),layerings_array(ig,islope,k,10))
         end do
     end do
@@ -752 +1022 @@
 do ig = 1,ngrid
     do islope = 1,nslope
-        if (is_co2ice_str(layerings_map(ig,islope)%top)) then
+        if (get_str_type(layerings_map(ig,islope)%top) == STR_TYPE_CO2ICE) then
             co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice*rho_co2ice
-        else if (is_h2oice_str(layerings_map(ig,islope)%top)) then
+        else if (get_str_type(layerings_map(ig,islope)%top) == STR_TYPE_H2OICE) then
             h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice*rho_h2oice
         else
-            call get_ssice_layering(layerings_map(ig,islope),h2oice_depth(ig,islope),h2o_ice(ig,islope),co2_ice(ig,islope))
+            call get_subice_layering(layerings_map(ig,islope),h2oice_depth(ig,islope),h2o_ice(ig,islope),co2_ice(ig,islope))
         end if
     end do
@@ -828 +1098 @@
 ! DEPENDENCIES
 ! ------------
-use utility, only: int2str
-use display, only: print_msg, LVL_DBG
-
-! DEPENDENCIES
-! ------------
 use geometry, only: ngrid, nslope
+use display,  only: print_msg, LVL_DBG
+use utility,  only: int2str
 
 ! DECLARATION
@@ -863 +1130 @@
 
 !=======================================================================
+SUBROUTINE define_str_grid(h,dz,nmin,nmax,z)
+!-----------------------------------------------------------------------
+! NAME
+!     define_str_grid
+!
+! DESCRIPTION
+!     Build a normalized-depth grid with adaptive point count.
+!
+! AUTHORS & DATE
+!     JB Clement, 04/2026
+!
+! NOTES
+!
+!-----------------------------------------------------------------------
+
+! DECLARATION
+! -----------
+implicit none
+
+! ARGUMENTS
+! ---------
+real(dp),                            intent(in)  :: h, dz
+integer(di),                         intent(in)  :: nmin, nmax
+real(dp), dimension(:), allocatable, intent(out) :: z
+
+! LOCAL VARIABLES
+! ---------------
+integer(di) :: n, i
+
+! CODE
+! ----
+n = nmin
+if (dz > eps) n = ceiling(max(0._dp,h)/dz)
+n = max(n,nmin)
+n = min(n,nmax)
+n = max(n,2_di)
+
+if (allocated(z)) deallocate(z)
+allocate(z(n))
+do i = 1,n
+    z(i) = real(i - 1,dp)/real(n - 1,dp)
+end do
+
+END SUBROUTINE define_str_grid
+!=======================================================================
+
+!=======================================================================
+FUNCTION is_poreice_empty(str) RESULT(profile_empty)
+!-----------------------------------------------------------------------
+! NAME
+!     is_poreice_empty
+!
+! DESCRIPTION
+!     Determine if a pore ice profile is empty.
+!
+! AUTHORS & DATE
+!     JB Clement, 04/2026
+!
+! NOTES
+!
+!-----------------------------------------------------------------------
+
+! DECLARATION
+! -----------
+implicit none
+
+! ARGUMENTS
+! ---------
+type(stratum), intent(in) :: str
+
+! LOCAL VARIABLES
+! ---------------
+logical(k4) :: profile_empty
+integer(di) :: i
+
+! CODE
+! ----
+profile_empty = .true.
+do i = 1,ncoef_poreice
+    if (abs(str%coef_poreice(i)) > tol) then
+        profile_empty = .false.
+        exit
+    end if
+end do
+
+END FUNCTION is_poreice_empty
+!=======================================================================
+
+!=======================================================================
+FUNCTION has_compact_ice(str) RESULT(has_ice)
+!-----------------------------------------------------------------------
+! NAME
+!     has_compact_ice
+!
+! DESCRIPTION
+!     Check if a stratum contains compact ice.
+!
+! AUTHORS & DATE
+!     JB Clement, 04/2026
+!
+! NOTES
+!
+!-----------------------------------------------------------------------
+
+! DECLARATION
+! -----------
+implicit none
+
+! ARGUMENTS
+! ---------
+type(stratum), intent(in) :: str
+
+! LOCAL VARIABLES
+! ---------------
+logical(k4) :: has_ice
+
+! CODE
+! ----
+has_ice = .false.
+if (str%h_co2ice > tol .or. str%h_h2oice > tol) has_ice = .true.
+
+END FUNCTION has_compact_ice
+!=======================================================================
+
+!=======================================================================
+FUNCTION is_pure_dust(str) RESULT(pure_dust)
+!-----------------------------------------------------------------------
+! NAME
+!     is_pure_dust
+!
+! DESCRIPTION
+!     Check if a stratum contains only dust (not cemented by ice).
+!
+! AUTHORS & DATE
+!     JB Clement, 04/2026
+!
+! NOTES
+!
+!-----------------------------------------------------------------------
+
+! DECLARATION
+! -----------
+implicit none
+
+! ARGUMENTS
+! ---------
+type(stratum), intent(in) :: str
+
+! LOCAL VARIABLES
+! ---------------
+logical(k4) :: pure_dust
+
+! CODE
+! ----
+pure_dust = .false.
+if (str%h_dust > tol .and. .not. has_compact_ice(str) .and. is_poreice_empty(str)) pure_dust = .true.
+
+END FUNCTION is_pure_dust
+!=======================================================================
+
+!=======================================================================
+FUNCTION is_dust_lag(str) RESULT(dust_lag)
+!-----------------------------------------------------------------------
+! NAME
+!     is_dust_lag
+!
+! DESCRIPTION
+!     Check if a stratum is a dust lag.
+!
+! AUTHORS & DATE
+!     JB Clement, 04/2026
+!
+! NOTES
+!
+!-----------------------------------------------------------------------
+
+! DECLARATION
+! -----------
+implicit none
+
+! ARGUMENTS
+! ---------
+type(stratum), intent(in) :: str
+
+! LOCAL VARIABLES
+! ---------------
+logical(k4) :: dust_lag
+
+! CODE
+! ----
+dust_lag = .false.
+if (str%h_dust > tol .and. .not. has_compact_ice(str)) dust_lag = .true.
+
+END FUNCTION is_dust_lag
+!=======================================================================
+
+!=======================================================================
+SUBROUTINE compute_poreice_dustlag(this)
+!-----------------------------------------------------------------------
+! NAME
+!     compute_poreice_dustlag
+!
+! DESCRIPTION
+!     Update top dust-lag pore-ice coefficients from dynamic pore filling.
+!
+! AUTHORS & DATE
+!     JB Clement, 04/2026
+!
+! NOTES
+!
+!-----------------------------------------------------------------------
+
+! DEPENDENCIES
+! ------------
+use subsurface_ice, only: dyn_ss_ice_m
+
+! DECLARATION
+! -----------
+implicit none
+
+! ARGUMENTS
+! ---------
+type(layering), intent(inout) :: this
+
+! LOCAL VARIABLES
+! ---------------
+
+! CODE
+! ----
+
+! Identify the contiguous top dust-lag column
+
+! Ensure the top dust-lag column overlays H2O ice
+
+! Build the grid for each dust-lag stratum
+
+! If the dust lag is too thin to be resolved by the grid, skip the update
+
+! Build the full grid and pore-ice profile for the dust-lag column
+
+! Compute the dynamic grid for the pore filling model, which is finer near the surface to better resolve the expected gradients
+
+! Interpolate the pore-ice profile onto the dynamic grid, assuming that all pores are filled below the dust lag
+
+! Get the updated pore-ice profile
+
+! Interpolate the updated pore-ice profile back onto the stratum-specific grids
+
+! Update the pore-ice coefficients for each dust-lag stratum
+
+
+END SUBROUTINE compute_poreice_dustlag
+!=======================================================================
+
+!=======================================================================
 FUNCTION thickness(this) RESULT(h_str)
 !-----------------------------------------------------------------------
@@ -898 +1420 @@
 
 !=======================================================================
-FUNCTION is_co2ice_str(str) RESULT(co2ice_str)
-!-----------------------------------------------------------------------
-! NAME
-!     is_co2ice_str
-!
-! DESCRIPTION
-!     Check if a stratum can be considered as a CO2 ice layer.
-!
-! AUTHORS & DATE
-!     JB Clement, 2024
-!
-! NOTES
-!
-!-----------------------------------------------------------------------
-
-! DECLARATION
-! -----------
-implicit none
-
-! ARGUMENTS
-! ---------
-type(stratum), intent(in) :: str
-
-! LOCAL VARIABLES
-! ---------------
-logical(k4) :: co2ice_str
-
-! CODE
-! ----
-co2ice_str = .false.
-if (str%h_co2ice > str%h_h2oice .and. str%h_co2ice > str%h_dust) co2ice_str = .true.
-
-END FUNCTION is_co2ice_str
-!=======================================================================
-
-!=======================================================================
-FUNCTION is_h2oice_str(str) RESULT(h2oice_str)
-!-----------------------------------------------------------------------
-! NAME
-!     is_h2oice_str
-!
-! DESCRIPTION
-!     Check if a stratum can be considered as a H2O ice layer.
-!
-! AUTHORS & DATE
-!     JB Clement, 2024
-!
-! NOTES
-!
-!-----------------------------------------------------------------------
-
-! DECLARATION
-! -----------
-implicit none
-
-! ARGUMENTS
-! ---------
-type(stratum), intent(in) :: str
-
-! LOCAL VARIABLES
-! ---------------
-logical(k4) :: h2oice_str
-
-! CODE
-! ----
-h2oice_str = .false.
-if (str%h_h2oice > str%h_co2ice .and. str%h_h2oice > str%h_dust) h2oice_str = .true.
-
-END FUNCTION is_h2oice_str
-!=======================================================================
-
-!=======================================================================
-FUNCTION is_dust_str(str) RESULT(dust_str)
-!-----------------------------------------------------------------------
-! NAME
-!     is_dust_str
-!
-! DESCRIPTION
-!     Check if a stratum can be considered as a dust layer.
-!
-! AUTHORS & DATE
-!     JB Clement, 2024
-!
-! NOTES
-!
-!-----------------------------------------------------------------------
-
-! DECLARATION
-! -----------
-implicit none
-
-! ARGUMENTS
-! ---------
-type(stratum), intent(in) :: str
-
-! LOCAL VARIABLES
-! ---------------
-logical(k4) :: dust_str
-
-! CODE
-! ----
-dust_str = .false.
-if (str%h_dust > str%h_co2ice .and. str%h_dust > str%h_h2oice) dust_str = .true.
-
-END FUNCTION is_dust_str
-!=======================================================================
-
-!=======================================================================
-FUNCTION is_dust_lag(str) RESULT(dust_lag)
-!-----------------------------------------------------------------------
-! NAME
-!     is_dust_lag
-!
-! DESCRIPTION
-!     Check if a stratum is a dust lag (no ice content).
-!
-! AUTHORS & DATE
-!     JB Clement, 2024
-!
-! NOTES
-!
-!-----------------------------------------------------------------------
-
-! DECLARATION
-! -----------
-implicit none
-
-! ARGUMENTS
-! ---------
-type(stratum), intent(in) :: str
-
-! LOCAL VARIABLES
-! ---------------
-logical(k4) :: dust_lag
-
-! CODE
-! ----
-dust_lag = .false.
-if (str%h_co2ice < tol .and. str%h_h2oice < tol .and. str%poreice_volfrac < tol) dust_lag = .true.
-
-END FUNCTION is_dust_lag
-!=======================================================================
-
-!=======================================================================
-SUBROUTINE get_ssice_layering(this,h_toplag,h2o_ice,co2_ice)
-!-----------------------------------------------------------------------
-! NAME
-!     get_ssice_layering
+SUBROUTINE get_subice_layering(this,h_toplag,h2o_ice,co2_ice)
+!-----------------------------------------------------------------------
+! NAME
+!     get_subice_layering
 !
 ! DESCRIPTION
@@ -1080 +1458 @@
 do
     h_toplag = h_toplag + thickness(current)
-    if (.not. associated(current%down)) then ! Bottom of the layering and only dust lag strata
+    if (.not. associated(current%down)) then ! Bottom of layering and only dust lag strata
         h_toplag = 0._dp
         return
@@ -1089 +1467 @@
 
 ! Is the stratum below the dust lag made of ice?
-if (current%h_h2oice > 0._dp .and. current%h_h2oice > current%h_co2ice) then ! Yes, there is more H2O ice than CO2 ice
+if (current%h_h2oice > tol .and. current%h_h2oice > current%h_co2ice) then ! Yes, there is more H2O ice than CO2 ice
     if (h_toplag < h_patchy_dust) then ! But the dust lag is too thin to considered ice as subsurface ice
         h_toplag = 0._dp
         h2o_ice = current%h_h2oice*rho_h2oice
     end if
-else if (current%h_co2ice > 0._dp .and. current%h_co2ice > current%h_h2oice) then ! No, there is more CO2 ice than H2O ice
+else if (current%h_co2ice > tol .and. current%h_co2ice > current%h_h2oice) then ! No, there is more CO2 ice than H2O ice
     if (h_toplag < h_patchy_ice) co2_ice = current%h_co2ice*rho_co2ice ! The dust lag is too thin to considered ice as subsurface ice
     h_toplag = 0._dp ! Because it matters only for H2O ice
 end if
 
-END SUBROUTINE get_ssice_layering
+END SUBROUTINE get_subice_layering
 !=======================================================================
 
@@ -1142 +1520 @@
 
 !=======================================================================
-SUBROUTINE sublimate_ice(this,str,h2subl_co2ice,h2subl_h2oice,new_lag,zlag)
+SUBROUTINE sublimate_ice(this,str,h2subl_co2ice,h2subl_h2oice,new_lag,zlag,current_date)
 !-----------------------------------------------------------------------
 ! NAME
@@ -1168 +1546 @@
 real(dp),               intent(inout) :: zlag
 real(dp),               intent(in)    :: h2subl_co2ice, h2subl_h2oice
-
-! LOCAL VARIABLES
-! ---------------
-real(dp)               :: h2subl, h_ice, h_pore, h_pore_new, h_tot
-real(dp)               :: hlag_dust, hlag_h2oice
-type(stratum), pointer :: current
+real(dp),               intent(in)    :: current_date
+
+! LOCAL VARIABLES
+! ---------------
+real(dp), dimension(ncoef_poreice), parameter :: coef_poreice_0 = 0._dp
+real(dp)                                      :: h2subl, h_ice, h_pore, h_pore_new, h_tot
+real(dp)                                      :: hlag_dust, hlag_h2oice
+type(stratum), pointer                        :: current
 
 ! CODE
@@ -1217 +1597 @@
 ! New stratum above the current stratum as a lag layer
 if (new_lag) then
-    call insert_stratum(this,str,str%top_elevation + h_tot,0._dp,hlag_h2oice,hlag_dust,h_pore_new,0._dp)
+    call insert_stratum(this,str,str%top_elevation + h_tot,0._dp,hlag_h2oice,hlag_dust,h_pore_new,coef_poreice_0,current_date)
     new_lag = .false.
 else
@@ -1268 +1648 @@
 h_toplag = 0._dp
 currentloc => current%up
-do while (is_dust_lag(currentloc) .and. associated(currentloc%up))
+do while (associated(currentloc%up) .and. is_dust_lag(currentloc))
     h_toplag = h_toplag + thickness(currentloc)
     currentloc => currentloc%up
@@ -1310 +1690 @@
 ! ARGUMENTS
 ! ---------
-real(dp),               intent(in)           :: h2oice_depth
-real(dp),               intent(inout)        :: flux_ssice, d_h2oice
-integer(di),            intent(in), optional :: icell, jcell
+real(dp),    intent(in)           :: h2oice_depth
+real(dp),    intent(inout)        :: flux_ssice, d_h2oice
+integer(di), intent(in), optional :: icell, jcell
 
 ! LOCAL VARIABLES
@@ -1367 +1747 @@
 
 !=======================================================================
-SUBROUTINE evolve_layering(this,d_co2ice,d_h2oice,h2oice_depth,flux_ssice,new_str,zshift_surf,new_lag,zlag,current,icell,jcell)
+SUBROUTINE evolve_layering(this,d_co2ice,d_h2oice,h2oice_depth,flux_ssice,new_str,zshift_surf,new_lag,zlag,current,current_date,icell,jcell)
 !-----------------------------------------------------------------------
 ! NAME
@@ -1398 +1778 @@
 ! ARGUMENTS
 ! ---------
+real(dp),               intent(in)           :: h2oice_depth
+real(dp),               intent(in)           :: current_date
+integer(di),            intent(in), optional :: icell, jcell
 real(dp),               intent(inout)        :: d_co2ice, d_h2oice, flux_ssice
-real(dp),               intent(in)           :: h2oice_depth
 type(layering),         intent(inout)        :: this
 type(stratum), pointer, intent(inout)        :: current
 logical(k4),            intent(inout)        :: new_str, new_lag
 real(dp),               intent(out)          :: zshift_surf, zlag
-integer(di),            intent(in), optional :: icell, jcell
-
-! LOCAL VARIABLES
-! ---------------
-real(dp)    :: dh_co2ice, dh_h2oice, dh_dust, dh_dust_co2, dh_dust_h2o
-real(dp)    :: h_co2ice_old, h_h2oice_old
-real(dp)    :: h2subl_co2ice, h2subl_h2oice, h2subl_co2ice_tot, h2subl_h2oice_tot
-real(dp)    :: h2lift, h2lift_tot, h_pore, h_tot!, R_subl
-logical(k4) :: unexpected
+
+! LOCAL VARIABLES
+! ---------------
+real(dp), dimension(ncoef_poreice), parameter :: coef_poreice_0 = 0._dp
+real(dp)                                      :: dh_co2ice, dh_h2oice, dh_dust, dh_dust_co2, dh_dust_h2o
+real(dp)                                      :: h_co2ice_old, h_h2oice_old
+real(dp)                                      :: h2subl_co2ice, h2subl_h2oice, h2subl_co2ice_tot, h2subl_h2oice_tot
+real(dp)                                      :: h2lift, h2lift_tot, h_pore, h_tot!, R_subl
+logical(k4)                                   :: unexpected
 
 ! CODE
@@ -1456 +1838 @@
     ! New stratum at the top of the layering
     if (new_str) then
-        call add_stratum(this,this%top%top_elevation + h_tot,dh_co2ice,dh_h2oice,dh_dust,h_pore,0._dp)
+        call add_stratum(this,this%top%top_elevation + h_tot,dh_co2ice,dh_h2oice,dh_dust,h_pore,coef_poreice_0,current_date)
         new_str = .false.
     else
@@ -1492 +1874 @@
             end if
             ! Ice sublimation
-            if (h2subl_co2ice > 0._dp .or. h2subl_h2oice > 0._dp) call sublimate_ice(this,current,h2subl_co2ice,h2subl_h2oice,new_lag,zlag)
+            if (h2subl_co2ice > 0._dp .or. h2subl_h2oice > 0._dp) call sublimate_ice(this,current,h2subl_co2ice,h2subl_h2oice,new_lag,zlag,current_date)
 
             ! Is there dust to lift?
             h2lift = 0._dp
-            if (is_dust_lag(current) .and. h2lift_tot > 0._dp) then
+            if (is_pure_dust(current) .and. h2lift_tot > 0._dp) then
                 h2lift = min(h2lift_tot,current%h_dust)
                 h2lift_tot = h2lift_tot - h2lift
@@ -1502 +1884 @@
                 if (h2lift > 0._dp) call lift_dust(this,current,h2lift)
             else if (associated(current%up)) then
-                if (is_dust_lag(current%up) .and. h2lift_tot > 0._dp) then
+                if (is_pure_dust(current%up) .and. h2lift_tot > 0._dp) then
                     h2lift = min(h2lift_tot,current%up%h_dust)
                     h2lift_tot = h2lift_tot - h2lift
@@ -1511 +1893 @@
 
             ! Still ice to be sublimated or dust to be lifted and no more ice in the current stratum so we move to the underlying stratum
-            if ((h2subl_co2ice_tot > 0. .or. h2subl_h2oice_tot > 0._dp .or. h2lift_tot > 0._dp) .and. current%h_co2ice < tol .and. current%h_h2oice < tol .and. current%poreice_volfrac < tol) then
+            if ((h2subl_co2ice_tot > 0. .or. h2subl_h2oice_tot > 0._dp .or. h2lift_tot > 0._dp) .and. &
+                current%h_co2ice < tol .and. current%h_h2oice < tol .and. is_poreice_empty(current)) then
                 current => current%down
                 new_lag = .true.
@@ -1550 +1933 @@
             end if
             ! Ice sublimation
-            if (h2subl_co2ice > 0._dp) call sublimate_ice(this,current,h2subl_co2ice,0._dp,new_lag,zlag)
+            if (h2subl_co2ice > 0._dp) call sublimate_ice(this,current,h2subl_co2ice,0._dp,new_lag,zlag,current_date)
 
             ! Is there dust to lift?
             h2lift = 0._dp
-            if (is_dust_lag(current) .and. h2lift_tot > 0._dp) then
+            if (is_pure_dust(current) .and. h2lift_tot > 0._dp) then
                 h2lift = min(h2lift_tot,current%h_dust)
                 h2lift_tot = h2lift_tot - h2lift
@@ -1560 +1943 @@
                 if (h2lift > 0._dp) call lift_dust(this,current,h2lift)
             else if (associated(current%up)) then
-                if (is_dust_lag(current%up) .and. h2lift_tot > 0._dp) then
+                if (is_pure_dust(current%up) .and. h2lift_tot > 0._dp) then
                     h2lift = min(h2lift_tot,current%up%h_dust)
                     h2lift_tot = h2lift_tot - h2lift
@@ -1585 +1968 @@
         ! New stratum at the top of the layering
         if (new_str) then
-            call add_stratum(this,this%top%top_elevation + h_tot,0._dp,dh_h2oice,dh_dust_h2o,h_pore,0._dp)
+            call add_stratum(this,this%top%top_elevation + h_tot,0._dp,dh_h2oice,dh_dust_h2o,h_pore,coef_poreice_0,current_date)
             new_str = .false.
         else
@@ -1599 +1982 @@
         ! New stratum at the top of the layering
         if (new_str) then
-            call add_stratum(this,this%top%top_elevation + h_tot,0._dp,dh_h2oice,dh_dust_h2o,h_pore,0._dp)
+            call add_stratum(this,this%top%top_elevation + h_tot,0._dp,dh_h2oice,dh_dust_h2o,h_pore,coef_poreice_0,current_date)
             new_str = .false.
         else
@@ -1625 +2008 @@
             end if
             ! Ice sublimation
-            if (h2subl_co2ice > 0._dp) call sublimate_ice(this,current,h2subl_co2ice,0._dp,new_lag,zlag)
+            if (h2subl_co2ice > 0._dp) call sublimate_ice(this,current,h2subl_co2ice,0._dp,new_lag,zlag,current_date)
 
             ! Is there dust to lift?
             h2lift = 0._dp
-            if (is_dust_lag(current) .and. h2lift_tot > 0._dp) then
+            if (is_pure_dust(current) .and. h2lift_tot > 0._dp) then
                 h2lift = min(h2lift_tot,current%h_dust)
                 h2lift_tot = h2lift_tot - h2lift
@@ -1635 +2018 @@
                 if (h2lift > 0._dp) call lift_dust(this,current,h2lift)
             else if (associated(current%up)) then
-                if (is_dust_lag(current%up) .and. h2lift_tot > 0.) then
+                if (is_pure_dust(current%up) .and. h2lift_tot > 0.) then
                     h2lift = min(h2lift_tot,current%up%h_dust)
                     h2lift_tot = h2lift_tot - h2lift
                     ! Dust lifting
-                    if (h2lift > 0.) call lift_dust(this,current%up,h2lift)
+                    if (h2lift > 0._dp) call lift_dust(this,current%up,h2lift)
                 end if
             end if
@@ -1662 +2045 @@
 
 if (unexpected) then
-    call print_msg('    dh_co2ice [m/y] = '//real2str(dh_co2ice),LVL_ERR)
-    call print_msg('    dh_h2oice [m/y] = '//real2str(dh_h2oice),LVL_ERR)
-    call print_msg('    dh_dust   [m/y] = '//real2str(dh_dust),LVL_ERR)
+    call print_msg('    dh_co2ice  [m/y] = '//real2str(dh_co2ice),LVL_ERR)
+    call print_msg('    dh_h2oice  [m/y] = '//real2str(dh_h2oice),LVL_ERR)
+    call print_msg('    dh_dust    [m/y] = '//real2str(dh_dust),LVL_ERR)
+    call print_msg('    flux_ssice [m/y] = '//real2str(flux_ssice),LVL_ERR)
     call stop_clean(__FILE__,__LINE__,'situation for the layering construction not managed!',1)
 end if

trunk/LMDZ.COMMON/libf/evolution/maths.F90

@@ -416 +416 @@
 
 !=======================================================================
+FUNCTION clamp(x,xmin,xmax) RESULT(y)
+!-----------------------------------------------------------------------
+! NAME
+!     clamp
+!
+! DESCRIPTION
+!     Clamp a scalar value into an admissible interval.
+!
+! AUTHORS & DATE
+!     JB Clement, 04/2026
+!
+! NOTES
+!
+!-----------------------------------------------------------------------
+
+! DECLARATION
+! -----------
+implicit none
+
+! ARGUMENTS
+! ---------
+real(dp), intent(in) :: x
+real(dp), intent(in) :: xmin, xmax
+
+! LOCAL VARIABLES
+! ---------------
+real(dp) :: y
+
+! CODE
+! ----
+if (xmin > xmax) then
+    y = x
+else
+    y = max(xmin,min(xmax,x))
+end if
+
+END FUNCTION clamp
+!=======================================================================
+
+!=======================================================================
+SUBROUTINE solve_linear_system(a,b,x)
+!-----------------------------------------------------------------------
+! NAME
+!     solve_linear_system
+!
+! DESCRIPTION
+!     Solve a dense linear system by Gaussian elimination.
+!
+! AUTHORS & DATE
+!     JB Clement, 04/2026
+!
+! NOTES
+!     Uses partial pivoting and returns ok = .false. if singular.
+!-----------------------------------------------------------------------
+
+! DEPENDENCIES
+! ------------
+use stoppage, only: stop_clean
+
+! DECLARATION
+! -----------
+implicit none
+
+! ARGUMENTS
+! ---------
+real(dp), dimension(:,:), intent(in)  :: a
+real(dp), dimension(:),   intent(in)  :: b
+real(dp), dimension(:),   intent(out) :: x
+
+! LOCAL VARIABLES
+! ---------------
+integer(di)                           :: n, i, j, k, ipiv
+real(dp)                              :: pivot, fac
+real(dp), dimension(:),   allocatable :: rhs, rowtmp
+real(dp), dimension(:,:), allocatable :: m
+
+! CODE
+! ----
+! Initialization
+n = size(b)
+x(:) = 0._dp
+if (size(a,1) /= n .or. size(a,2) /= n) call stop_clean(__FILE__,__LINE__,"size of a must be n x n!",1)
+if (size(x) /= n) call stop_clean(__FILE__,__LINE__,"size of x must match size of b!",1)
+
+allocate(m(n,n),rhs(n),rowtmp(n))
+m(:,:) = a(:,:)
+rhs(:) = b(:)
+
+! Gaussian elimination with partial pivoting
+do k = 1,n - 1
+    ipiv = k
+    pivot = abs(m(k,k))
+    do i = k + 1,n
+        if (abs(m(i,k)) > pivot) then
+            pivot = abs(m(i,k))
+            ipiv = i
+        end if
+    end do
+    if (pivot <= eps) call stop_clean(__FILE__,__LINE__,"matrix is singular!",1)
+
+    if (ipiv /= k) then
+        rowtmp(:) = m(k,:)
+        m(k,:) = m(ipiv,:)
+        m(ipiv,:) = rowtmp(:)
+        fac = rhs(k)
+        rhs(k) = rhs(ipiv)
+        rhs(ipiv) = fac
+    end if
+
+    do i = k + 1,n
+        fac = m(i,k)/m(k,k)
+        m(i,k:n) = m(i,k:n) - fac*m(k,k:n)
+        rhs(i) = rhs(i) - fac*rhs(k)
+    end do
+end do
+
+if (abs(m(n,n)) <= eps) call stop_clean(__FILE__,__LINE__,"matrix is singular!",1)
+
+do i = n,1,-1
+    x(i) = rhs(i)
+    do j = i + 1,n
+        x(i) = x(i) - m(i,j)*x(j)
+    end do
+    if (abs(m(i,i)) <= eps) call stop_clean(__FILE__,__LINE__,"matrix is singular!",1)
+    x(i) = x(i)/m(i,i)
+end do
+
+END SUBROUTINE solve_linear_system
+!=======================================================================
+
+!=======================================================================
+SUBROUTINE coef2prof(coef,x,fx)
+!-----------------------------------------------------------------------
+! NAME
+!     coef2prof
+!
+! DESCRIPTION
+!     Evaluate polynomial coefficients on a 1D profile.
+!
+! AUTHORS & DATE
+!     JB Clement, 04/2026
+!
+! NOTES
+!     Coefficients are ordered from degree 0 to degree n-1.
+!-----------------------------------------------------------------------
+
+! DEPENDENCIES
+! -----------
+use stoppage, only: stop_clean
+
+! DECLARATION
+! -----------
+implicit none
+
+! ARGUMENTS
+! ---------
+real(dp), dimension(:), intent(in)  :: coef
+real(dp), dimension(:), intent(in)  :: x
+real(dp), dimension(:), intent(out) :: fx
+
+! LOCAL VARIABLES
+! ---------------
+integer(di) :: i, j, ncoef
+real(dp)    :: yloc
+
+! CODE
+! ----
+! Initialization
+fx(:) = 0._dp
+if (size(fx) /= size(x)) call stop_clean(__FILE__,__LINE__,"size of fx must match size of x!",1)
+
+ncoef = int(size(coef),di)
+
+! Evaluate the polynomial at each point using Horner's method for numerical stability
+do i = 1,size(x)
+    yloc = coef(ncoef)
+    do j = ncoef - 1,1,-1
+        yloc = yloc*x(i) + coef(j)
+    end do
+    fx(i) = yloc
+end do
+
+END SUBROUTINE coef2prof
+!=======================================================================
+
+!=======================================================================
+SUBROUTINE prof2coef(x,fx,coef)
+!-----------------------------------------------------------------------
+! NAME
+!     prof2coef
+!
+! DESCRIPTION
+!     Fit polynomial coefficients from a 1D profile by least squares.
+!
+! AUTHORS & DATE
+!     JB Clement, 04/2026
+!
+! NOTES
+!     Coefficients are ordered from degree 0 to degree n-1.
+!-----------------------------------------------------------------------
+
+! DEPENDENCIES
+! ------------
+use stoppage, only: stop_clean
+
+! DECLARATION
+! -----------
+implicit none
+
+! ARGUMENTS
+! ---------
+real(dp), dimension(:), intent(in)  :: x
+real(dp), dimension(:), intent(in)  :: fx
+real(dp), dimension(:), intent(out) :: coef
+
+! LOCAL VARIABLES
+! ---------------
+integer(di)                           :: i, j, k, ncoef
+real(dp), dimension(:),   allocatable :: rhs, basis
+real(dp), dimension(:,:), allocatable :: normal
+
+! CODE
+! ----
+! Initialization
+coef(:) = 0._dp
+ncoef = int(size(coef),di)
+if (size(fx) /= size(x)) call stop_clean(__FILE__,__LINE__,"size of fx must match size of x!",1)
+if (size(x) < ncoef) call stop_clean(__FILE__,__LINE__,"size of x must be at least as large as the number of coefficients!",1)
+
+allocate(rhs(ncoef),basis(ncoef),normal(ncoef,ncoef))
+normal(:,:) = 0._dp
+rhs(:) = 0._dp
+
+! Set up the normal equations for least squares fitting
+do i = 1,size(x)
+    basis(1) = 1._dp
+    do j = 2,ncoef
+        basis(j) = basis(j - 1)*x(i)
+    end do
+
+    do j = 1,ncoef
+        rhs(j) = rhs(j) + basis(j)*fx(i)
+        do k = 1,ncoef
+            normal(j,k) = normal(j,k) + basis(j)*basis(k)
+        end do
+    end do
+end do
+
+! Solve the linear system
+call solve_linear_system(normal,rhs,coef)
+
+END SUBROUTINE prof2coef
+!=======================================================================
+
+!=======================================================================
+FUNCTION interp_profile(x_nodes,y_nodes,x_eval) RESULT(y_eval)
+!-----------------------------------------------------------------------
+! NAME
+!     interp_profile
+!
+! DESCRIPTION
+!     Interpolate a profile given at nodes.
+!
+! AUTHORS & DATE
+!     JB Clement, 04/2026
+!
+! NOTES
+!     Performs piecewise linear interpolation.
+!-----------------------------------------------------------------------
+
+! DECLARATION
+! -----------
+implicit none
+
+! ARGUMENTS
+! ---------
+real(dp), dimension(:), intent(in) :: x_nodes, y_nodes
+real(dp),               intent(in) :: x_eval
+
+! LOCAL VARIABLES
+! ---------------
+real(dp)    :: y_eval, weight
+integer(di) :: i, n
+
+! CODE
+! ----
+! Initialization
+n = int(size(x_nodes),di)
+y_eval = 0._dp
+if (n <= 0_di .or. size(y_nodes) /= n) return
+
+! Handle out-of-bounds cases by clamping to the nearest node value
+if (x_eval <= x_nodes(1)) then
+    y_eval = y_nodes(1)
+    return
+end if
+if (x_eval >= x_nodes(n)) then
+    y_eval = y_nodes(n)
+    return
+end if
+
+! Find the interval containing x_eval and perform linear interpolation
+do i = 1,n - 1
+    if (x_eval <= x_nodes(i + 1)) then
+        if (abs(x_nodes(i + 1) - x_nodes(i)) <= eps) then ! Avoid division by zero if nodes are too close
+            y_eval = y_nodes(i)
+        else
+            weight = (x_eval - x_nodes(i))/(x_nodes(i + 1) - x_nodes(i))
+            y_eval = weight*(y_nodes(i + 1) - y_nodes(i)) + y_nodes(i)
+        end if
+        return
+    end if
+end do
+
+y_eval = y_nodes(n)
+
+END FUNCTION interp_profile
+!=======================================================================
+
+!=======================================================================
 FUNCTION limit_slope(dqL,dqR,limiter) RESULT(limited_slope)
 !-----------------------------------------------------------------------
@@ -426 +746 @@
 !
 ! AUTHORS & DATE
-!     JB Clement, 02/s2026
+!     JB Clement, 02/2026
 !
 ! NOTES
@@ -480 +800 @@
 !
 ! AUTHORS & DATE
-!     JB Clement, 02/s2026
+!     JB Clement, 02/2026
 !
 ! NOTES

trunk/LMDZ.COMMON/libf/evolution/pem.F90

@@ -36 +36 @@
 use glaciers,           only: h2oice_flow, co2ice_flow, flow_co2glaciers, flow_h2oglaciers
 use ice_table,          only: icetable_equilibrium, icetable_dynamic, evolve_ice_table
-use layered_deposits,   only: do_layering, del_layering, evolve_layering, ptrarray, layering2surfice, surfice2layering, print_layering
+use layered_deposits,   only: do_layering, del_layering, evolve_layering, ptrarray, layering2surfice, surfice2layering, &
+                              print_layering, compute_poreice_dustlag
 use maths,              only: pi
 use numerics,           only: dp, qp, di, li, k4, eps, eps_qp
@@ -238 +239 @@
             do i = 1,ngrid
                 call evolve_layering(layerings_map(i,islope),d_co2ice(i,islope),d_h2oice(i,islope),h2oice_depth(i,islope),flux_ssice_avg(i,islope), &
-                                     new_str(i,islope),zshift_surf(i,islope),new_lag(i,islope),zlag(i,islope),current(i,islope)%p,i,islope)
+                                     new_str(i,islope),zshift_surf(i,islope),new_lag(i,islope),zlag(i,islope),current(i,islope)%p,                  &
+                                     pem_ini_date + n_yr_sim + n_yr_run + dt,i,islope)
+                !if (do_soil .and. icetable_dynamic) then
+                !    call compute_poreice_dustlag(layerings_map(i,islope))
+                !end if
                 if (is_lvl_enabled(LVL_DBG)) call print_layering(layerings_map(i,islope))
             end do

trunk/LMDZ.COMMON/libf/evolution/soil.F90

@@ -287 +287 @@
 if (size(soildepth_PCM_in) /= nsoil_PCM) call stop_clean(__FILE__,__LINE__,'soil size mismatch between "startfi.nc" and geometry dimensions!',1)
 
-! Set PCM soil depths read from "startfi.nc" to the mid-layer depth array
+! Set PCM soil layer depths read from "startfi.nc" to the mid-layer depth array
 mlayer(0:nsoil_PCM - 1) = soildepth_PCM_in(:)
 
@@ -335 +335 @@
 ! ---------------
 integer(di) :: ig, iloop, islope
-real(dp)    :: alpha, lay1       ! Coefficients for building layers
-real(dp)    :: index_powerlaw    ! Coefficient to match the PEM power law grid to the PCM grid
+real(dp)    :: alpha, lay1    ! Coefficients to build soil layer depths
+real(dp)    :: index_powerlaw ! Coefficient to match the PEM power law grid to the PCM grid
 
 ! CODE

trunk/LMDZ.COMMON/libf/evolution/subsurface_ice.F90

@@ -14 +14 @@
 !
 ! NOTES
-!    Based on Norbert Schorgofer's code. Updated for PEM.
+!    Wrapper for Norbert Schorgofer's code that has been slightly updated
+!    to work with the PEM.
 !-----------------------------------------------------------------------
 
@@ -39 +40 @@
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
-
-!=======================================================================
-!============================ dyn_ss_ice_m =============================
-!=======================================================================
-
-
+!=======================================================================
+!-----------------------------------------------------------------------
+! NAME
+!     dyn_ss_ice_m
+!
+! DESCRIPTION
+!     Advance subsurface ice state for one column (NP=1) using the fast
+!     dynamic method: update pore filling, ice table depth, and soil
+!     temperatures over a short integration.
+!
+! INPUTS
+!     ssi_depth_in      Initial ice-table depth [m]
+!     T1                Mean/representative surface temperature [K]
+!     Tb(nz)            Initial soil temperature profile [K]
+!     nz                Number of vertical levels
+!     thIn(1)           Thermal inertia [SI]
+!     p0(1)             Total atmospheric pressure [Pa]
+!     pfrost(1)         H2O vapor partial pressure proxy [Pa]
+!     porefill_in(nz,1) Initial pore-ice filling fraction [0-1]
+!
+! OUTPUTS
+!     ssi_depth         Updated ice-table depth [m]
+!     porefill(nz,1)    Updated pore-ice filling fraction [0-1]
+!     Tb(nz)            Updated soil temperature profile [K]
+!
+! NOTES
+!     - This wrapper is hard-coded for NP=1.
+!     - avrho1 is prescribed as pfrost/T1 in the underlying solver.
+!-----------------------------------------------------------------------
 SUBROUTINE dyn_ss_ice_m(ssi_depth_in,T1,Tb,nz,thIn,p0,pfrost,porefill_in,porefill,ssi_depth)
-
 !***********************************************************************
 ! Retreat and growth of subsurface ice on Mars
@@ -53 +76 @@
   use orbit, only: sol_len_s
   use maths, only: pi
-  ! DECLARATION
-! -----------
-implicit none
+  implicit none
   integer, parameter :: NP=1   ! # of sites
   integer nz, i, k, iloop
@@ -208 +229 @@
 
 end subroutine dyn_ss_ice_m
-
-!=======================================================================
-
+!=======================================================================
+
+!=======================================================================
 subroutine soilthprop(porosity,fill,rhocobs,tiobs,layertype, &
      &     newrhoc,newti,icefrac)
@@ -226 +247 @@
 !     output are newti and newrhoc
 !***********************************************************************
-
-! DEPENDENCIES
-! ------------
-use stoppage, only: stop_clean
-
-! DECLARATION
-! -----------
-implicit none
+  use stoppage, only: stop_clean
+  implicit none
   integer, intent(IN) :: layertype
   real(8), intent(IN) :: porosity, fill, rhocobs, tiobs
@@ -295 +310 @@
 
 end subroutine soilthprop
-
-
-!=======================================================================
-
+!=======================================================================
+
+!=======================================================================
       real*8 function frostpoint(p)
 !     inverse of psv
 !     input is partial pressure [Pascal]
 !     output is temperature [Kelvin]
-      ! DECLARATION
-! -----------
-implicit none
+      implicit none
       real*8 p
 
@@ -326 +338 @@
 
       end function frostpoint
-
-
-!=======================================================================
-!=========================== fast_subs_mars ============================
-!=======================================================================
-
-
+!=======================================================================
+
+!=======================================================================
 !*****************************************************
 ! Subroutines for fast method
 ! written by Norbert Schorghofer 2007-2011
 !*****************************************************
-
-
 subroutine icelayer_mars(bigstep,nz,NP,thIn,rhoc,z,porosity,pfrost, &
      & Tb,zdepthF,zdepthE,porefill,Tmean1,Tmean3,zdepthG, &
@@ -347 +353 @@
 ! latitude  [degree]
 !*************************************************************************
-  !use ice_table,      only: rho_ice
+  !use ice_table, only: rho_ice
   !use omp_lib
-  ! DECLARATION
-! -----------
-implicit none
+  implicit none
   integer, intent(IN) :: nz, NP
   real(8), intent(IN) :: bigstep
@@ -453 +457 @@
   !$omp end parallel
 end subroutine icelayer_mars
-
-
-
+!=======================================================================
+
+!=======================================================================
 subroutine ajsub_mars(typeT, albedo0, pfrost, nz, z, ti, rhocv, &
      &     fracIR, fracDust, patm, avdrho, avdrhoP, avrho1, &
@@ -467 +471 @@
 !***********************************************************************
   use orbit, only: yr_len_sol, sol_len_s
-  ! DECLARATION
-! -----------
-implicit none
+  implicit none
   integer, intent(IN) :: nz, typeT
 !  real(8), intent(IN) :: latitude  ! in radians
@@ -616 +618 @@
 
 end subroutine ajsub_mars
-
-
-
+!=======================================================================
+
+!=======================================================================
 real*8 function tfrostco2(p)
 !     the inverse of function psvco2
 !     input is pressure in Pascal, output is temperature in Kelvin
-! DECLARATION
-! -----------
 implicit none
 real*8 p
 tfrostco2 = 3182.48/(23.3494+log(100./p))
 end function
-
-!=======================================================================
-
+!=======================================================================
+
+!=======================================================================
       real*8 function psv(T)
 !     saturation vapor pressure of H2O ice [Pascal]
 !     input is temperature [Kelvin]
-      ! DECLARATION
-! -----------
-implicit none
+      implicit none
       real*8 T
 
@@ -659 +657 @@
 
       end function psv
-
-
-!=======================================================================
-!=========================== fast_subs_univ ============================
-!=======================================================================
-
-
+!=======================================================================
+
+!=======================================================================
 !*****************************************************
 ! Commonly used subroutines for fast method
 ! written by Norbert Schorghofer 2007-2011
 !*****************************************************
-
 pure function zint(y1,y2,z1,z2)
   ! interpolate between two values, y1*y2<0
-  ! DECLARATION
-! -----------
-implicit none
+  implicit none
   real(8), intent(IN) :: y1,y2,z1,z2
   real(8) zint
   zint = (y1*z2 - y2*z1)/(y1-y2)
 end function zint
-
-
-
+!=======================================================================
+
+!=======================================================================
 pure function equildepth(nz, z, rhosatav, rhosatav0, avrho1)
 !***********************************************************************
@@ -688 +679 @@
 !  this is not the true (final) equilibrium depth
 !***********************************************************************
-  ! DECLARATION
-! -----------
-implicit none
+  implicit none
   integer, intent(IN) :: nz
   real(8), intent(IN) :: z(nz), rhosatav(nz)
@@ -711 +700 @@
   if (equildepth>z(nz)) equildepth=-9999.   ! happens very rarely
 end function equildepth
-
-
-
+!=======================================================================
+
+!=======================================================================
 subroutine depths_avmeth(typeT, nz, z, rhosatav, rhosatav0, rlow, avrho1,  &
      & porefill, typeF, zdepthF, B, ypp, typeG, zdepthG)
@@ -724 +713 @@
 !  B = Diff*bigstep/(porosity*icedensity)  [SI units]
 !***********************************************************************
-  use maths,     only: deriv2_simple, deriv1_onesided, deriv1, colint
-  ! DECLARATION
-! -----------
-implicit none
+  use maths, only: deriv2_simple, deriv1_onesided, deriv1, colint
+  implicit none
   integer, intent(IN) :: nz, typeT
   real(8), intent(IN), dimension(nz) :: z, rhosatav, porefill
@@ -846 +833 @@
   ! if typeG>0, then all ice at and below typeG should be erased
 end subroutine depths_avmeth
-
-
-
+!=======================================================================
+
+!=======================================================================
 pure subroutine icechanges(nz,z,typeF,avdrho,avdrhoP,ypp, &
      & Diff,porosity,icefrac,bigstep,zdepthT,porefill,typeG)
@@ -856 +843 @@
 ! a crucial subroutine
 !***********************************************************
-  use maths,     only: colint
+  use maths, only: colint
   !use ice_table, only: rho_ice
-  ! DECLARATION
-! -----------
-implicit none
+  implicit none
   integer, intent(IN) :: nz, typeF, typeG
   real(8), intent(IN) :: z(nz), ypp(nz), avdrho, avdrhoP
@@ -940 +925 @@
   end if
 end subroutine icechanges
+!=======================================================================
 
 !=======================================================================

trunk/LMDZ.COMMON/libf/misc/job.F90

@@ -9 +9 @@
 !
 ! AUTHORS & DATE
-!     JB Clement, 07/07/2025
+!     JB Clement, 10/06/2024
 !
 ! NOTES
@@ -38 +38 @@
 !
 ! AUTHORS & DATE
-!     JB Clement, 07/07/2025
+!     JB Clement, 10/06/2024
 !
 ! NOTES
@@ -90 +90 @@
 !
 ! AUTHORS & DATE
-!     JB Clement, 07/07/2025
+!     JB Clement, 10/06/2024
 !
 ! NOTES

trunk/LMDZ.COMMON/makelmdz_fcm

@@ -659 +659 @@
 !
 ! AUTHORS & DATE
-!     JB Clement, 07/07/2025
+!     JB Clement, 13/01/2025
 !
 ! NOTES
@@ -703 +703 @@
 !
 ! AUTHORS & DATE
-!     JB Clement, 07/07/2025
+!     JB Clement, 13/01/2025
 !
 ! NOTES