Changeset 3926


Ignore:
Timestamp:
Oct 7, 2025, 4:31:36 PM (3 months ago)
Author:
jbclement
Message:

PEM:
Updates of files in the deftank + parameter values for the layering algorithm.
JBC

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

Legend:

Unmodified
Added
Removed

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

    r3909 r3926  
    776776== 03/09/2025 == JBC
    777777Managing automatically the slope variables outputted by XIOS according to 'nslope' to avoid subsequent errors.
     778
     779== 07/10/2025 == JBC
     780Updates of files in the deftank + parameter values for the layering algorithm.

  • trunk/LMDZ.COMMON/libf/evolution/deftank/clean.sh

    r3861 r3926  
    2727rm -f restart1D*.txt
    2828rm -f start1D*.txt
     29rm -f profile_temp.out
    2930# Reset the initial starting files to prepare a new simulation
    3031if [ -f "starts/startfi.nc" ]; then

  • trunk/LMDZ.COMMON/libf/evolution/deftank/launchPEM.sh

    r3861 r3926  
    7979            errlaunch
    8080        fi
    81         echo "The relaunch is initialized with a specific previous successful run."
     81        echo "The relaunch is initialized with a previous successful run to be specified."
    8282        while true; do
    8383            echo "Do you want to relaunch from a 'PCM' or 'PEM' run?"

  • trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def

    r3842 r3926  
    105105# layering=.false.
    106106 
    107 # Value of the dust tendency. Default = 1.e-3 kg.m-2.y-1
    108 # d_dust=1.e-3
     107# Value of the dust tendency. Default = 5.78e-2 kg.m-2.y-1
     108# d_dust=5.78e-2
    109109 
    110110# Do you want to impose a dust-to-ice ratio instead of a dust tendency? Default = .false.
    111111# impose_dust_ratio=.false.
    112112 
    113 # If impose_dust_ratio=.true., value of the dust-to-ice ratio. Default = 0.01
    114 # dust2ice_ratio=0.01
     113# If impose_dust_ratio=.true., value of the dust-to-ice ratio. Default = 0.1
     114# dust2ice_ratio=0.1
    115115
    116116# Some definitions for the physics, in file 'callphys.def'

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

    r3883 r3926  
    1313from mpl_toolkits.axes_grid1.inset_locator import inset_axes
    1414from matplotlib.colors import LinearSegmentedColormap, LogNorm
     15from matplotlib.ticker import FuncFormatter
    1516from scipy.interpolate import interp1d
    1617
    … …  
    576577            # Compute RGB stratification over time
    577578            rgb = np.ones((nz, ntime, 3), dtype=float)
    578 
    579579            frac_all = np.full((nz, ntime, 3), np.nan, dtype=float)  # store fH2O, fCO2, fDust
    580580            for t in range(ntime):
    … …  
    590590                fCO2 = cCO2 / total
    591591                fDust = cDust / total
    592                 frac_all[:depth, t, :] = np.stack([fH2O, fCO2, fDust], axis=1)
    593                 mix = np.outer(fH2O, blue) + np.outer(fCO2, violet) + np.outer(fDust, orange)
    594                 rgb[:depth, t, :] = np.clip(mix, 0, 1)
     592                frac_all[:depth, t, 0] = fH2O
     593                frac_all[:depth, t, 1] = fCO2
     594                frac_all[:depth, t, 2] = fDust
     595                rgb[:depth, t, 0] = fH2O * blue[0] + fCO2 * violet[0] + fDust * orange[0]
     596                rgb[:depth, t, 1] = fH2O * blue[1] + fCO2 * violet[1] + fDust * orange[1]
     597                rgb[:depth, t, 2] = fH2O * blue[2] + fCO2 * violet[2] + fDust * orange[2]
    595598
    596599            # Mask elevation
    … …  
    701704    """
    702705    h2o = gridded_data['h2o_ice']
     706    co2 = gridded_data['co2_ice']
    703707    dust = gridded_data['dust']
    704708    ngrid, ntime, nslope, nz = h2o.shape
    705 
    706     # Elevation mask
    707     if exclude_sub:
    708         elevation_mask = (ref_grid >= 0.0)
    709         elev = ref_grid[elevation_mask]
    710     else:
    711         elevation_mask = np.ones_like(ref_grid, dtype=bool)
    712         elev = ref_grid
    713709
    714710    # Define custom blue-to-orange colormap
    … …  
    729725            log_ratio_array = np.full((nz, ntime), np.nan, dtype=np.float32)
    730726            for t in range(ntime):
    731                 if ti[t] <= 0:
     727                if t > 0 and ti[t] <= 0:
    732728                    ti[t] = ti[t-1]
     729                elif ti[t] <= 0:
     730                    continue
    733731                zmax = ti[t]
    734732                if zmax <= 0:
    … …  
    750748                log_ratio_array[:zmax, t] = log_ratio
    751749
    752             # Mask and prepare for display
    753750            ratio_array = 10**log_ratio_array
    754             ratio_display = ratio_array[elevation_mask, :]
    755             x_edges = np.concatenate([date_time, [date_time[-1] + (date_time[-1] - date_time[-2])]])
    756             y_edges = np.concatenate([elev, [elev[-1] + (elev[-1] - elev[-2])]])
     751
     752            # Compute edges for pcolormesh
     753            x_edges = np.concatenate([date_time, [date_time[-1] + (date_time[-1] - date_time[-2])]]) * martian_to_earth
     754            y_edges = np.concatenate([ref_grid, [ref_grid[-1] + (ref_grid[-1] - ref_grid[-2])]])
    757755
    758756            # Plot
    … …  
    761759                date_time,
    762760                elev,
    763                 ratio_display,
     761                ratio_array,
    764762                shading='auto',
    765763                cmap='managua_r',
    766764                norm=LogNorm(vmin=10**vmin, vmax=10**vmax),
    767765            )
    768             attach_format_coord(ax, ratio_display, x_edges, y_edges, is_pcolormesh=True)
     766            attach_format_coord(ax, ratio_array, x_edges, y_edges, is_pcolormesh=True)
    769767            ax.set_title(f"Dust-to-Ice ratio over time (Grid point {ig+1}, Slope {isl+1})", fontweight='bold')
    770768            ax.set_xlabel('Time (Mars years)')
    … …  
    966964    ls_perihelion = mars_ls(date_peri_day,0.,eccentricity,year_day)
    967965
    968     return dates_yr, obliquity, eccentricity, ls_perihelion
     966    return dates_yr, obliquity, eccentricity, ls_perihelion, martian_to_earth
    969967
    970968
    … …  
    972970    """
    973971    Plot the evolution of obliquity, eccentricity and Ls p coming from simulation data
     972    versus simulated time, plus an additional figure of sin(eccentricity)*Lsp.
    974973    versus simulated time.
    975974    """
    976975    # Read orbital data
    977     times_yr, obl, ecc, lsp = read_orbital_data_nc(starts_folder, infofile)
     976    times_yr, obl, ecc, lsp, martian_to_earth = read_orbital_data_nc(starts_folder, infofile)
    978977
    979978    fargs = dict(kind='linear', bounds_error=False, fill_value='extrapolate')
    … …  
    982981    lsp_i = interp1d(times_yr, lsp, **fargs)(date_time)
    983982
     983    date_time = date_time * martian_to_earth / 1e6
     984
    984985    fig, axes = plt.subplots(3,1, figsize=(8,10), sharex=True)
    985986    fig.suptitle("Orbital parameters vs simulated time", fontsize=14, fontweight='bold')
    … …  
    996997    axes[2].plot(date_time, lsp_i, 'g-', marker='+')
    997998    axes[2].set_ylabel("Ls of perihelion (°)")
    998     axes[2].set_xlabel("Time (Mars years)")
     999    axes[2].set_xlabel("Time (Myr)")
    9991000    axes[2].grid(True)
    10001001
    … …  
    10021003    outname = os.path.join(output_folder, "orbital_parameters_simu.png")
    10031004    fig.savefig(outname, dpi=150)
     1005
     1006    eps_sin_lsp = ecc_i * np.sin(np.radians(lsp_i))
     1007
     1008    fig2, ax2 = plt.subplots(figsize=(8,5))
     1009    fig2.suptitle(r"$\epsilon \times \sin(L_{sp})$", fontweight='bold')
     1010
     1011    ax2.plot(date_time, eps_sin_lsp, 'm-', marker='+')
     1012    ax2.set_ylabel(r"$\epsilon \cdot \sin(L_{sp})$")
     1013    ax2.set_xlabel("Time (Myr)")
     1014    ax2.grid(True)
     1015
     1016    plt.tight_layout(rect=[0,0,1,0.95])
     1017    outname2 = os.path.join(output_folder, "sin_ecc_times_Lsp.png")
     1018    fig2.savefig(outname2, dpi=150)
    10041019
    10051020
    … …  
    10251040
    10261041    # Read orbital data
    1027     times_yr, obl, _, _ = read_orbital_data_nc(starts_folder, infofile)
     1042    times_yr, obl, _, _, martian_to_earth = read_orbital_data_nc(starts_folder, infofile)
    10281043    fargs = dict(kind='linear', bounds_error=False, fill_value='extrapolate')
    10291044    obliquity = interp1d(times_yr, obl, **fargs)(date_time)
    10301045
    1031     # Computed total height
     1046    # Define custom blue-to-orange colormap
     1047    blue = np.array([0, 0, 255], dtype=float) / 255
     1048    orange = np.array([255, 165, 0], dtype=float) / 255
     1049    custom_cmap = LinearSegmentedColormap.from_list('BlueOrange', [blue, orange], N=256)
     1050    color_map = { 1: 'green', -1: 'red', 0: 'orange' }
     1051
     1052    # Log‑ratio bounds and small epsilon to avoid log(0)
     1053    vmin, vmax = -2, 1
     1054    epsilon = 1e-6
     1055
     1056    # Loop over grids and slopes
    10321057    for ig in range(ngrid):
    10331058        for isl in range(nslope):
     1059            # Compute total height time series
    10341060            total_height_t = np.zeros(ntime, dtype=float)
    1035 
    10361061            for t_idx in range(ntime):
    10371062                h_mat = heights_data[t_idx][isl]
    … …  
    10421067                    total_height_t[t_idx] = np.max(h_valid)
    10431068
    1044     # Compute the per-interval sign of height change
    1045     dh = np.diff(total_height_t)
    1046     signs = np.sign(dh)
    1047     color_map = { 1: 'green', -1: 'red', 0: 'orange' }
    1048 
    1049     # Elevation mask
    1050     if exclude_sub:
    1051         elevation_mask = (ref_grid >= 0.0)
    1052         elev = ref_grid[elevation_mask]
    1053     else:
    1054         elevation_mask = np.ones_like(ref_grid, dtype=bool)
    1055         elev = ref_grid
    1056 
    1057     # Define custom blue-to-orange colormap
    1058     blue = np.array([0, 0, 255], dtype=float) / 255
    1059     orange = np.array([255, 165, 0], dtype=float) / 255
    1060     custom_cmap = LinearSegmentedColormap.from_list('BlueOrange', [blue, orange], N=256)
    1061 
    1062     # Log‑ratio bounds and small epsilon to avoid log(0)
    1063     vmin, vmax = -2, 1
    1064     epsilon = 1e-6
    1065 
    1066     # Loop over grids and slopes
    1067     for ig in range(ngrid):
    1068         for isl in range(nslope):
     1069            # Compute the per-interval sign of height change
     1070            if ntime > 1:
     1071                dh = np.diff(total_height_t)
     1072                signs = np.sign(dh).astype(int)
     1073            else:
     1074                dh = np.array([], dtype=float)
     1075                signs = np.array([], dtype=int)
     1076
    10691077            # Prepare fraction and ratio arrays
    10701078            ti = top_index[ig, :, isl].copy().astype(int)
    … …  
    10721080            frac_all = np.full((nz, ntime, 3), np.nan, dtype=float)  # store fH2O, fCO2, fDust
    10731081            for t in range(ntime):
    1074                 if ti[t] <= 0:
     1082                if t > 0 and ti[t] <= 0:
    10751083                    ti[t] = ti[t-1]
     1084                elif ti[t] <= 0:
     1085                    continue
    10761086                zmax = ti[t]
    10771087                if zmax <= 0:
    … …  
    10861096                fCO2 = cCO2 / total
    10871097                fDust = cDust / total
    1088                 frac_all[:zmax, t, :] = np.stack([fH2O, fCO2, fDust], axis=1)
     1098                frac_all[:zmax, t, 0] = fH2O
     1099                frac_all[:zmax, t, 1] = fCO2
     1100                frac_all[:zmax, t, 2] = fDust
    10891101
    10901102                with np.errstate(divide='ignore', invalid='ignore'):
    1091                     ratio = np.where(
    1092                         cH2O > 0,
    1093                         cDust / cH2O,
    1094                         10**(vmax + 1)
     1103                    ratio = np.where(cH2O > 0, cDust / cH2O, 10**(vmax + 1)
    10951104                    )
    10961105                    log_ratio = np.log10(ratio + epsilon)
    … …  
    10991108                log_ratio_array[:zmax, t] = log_ratio
    11001109
    1101             # Mask elevation
    11021110            ratio_array = 10**log_ratio_array
    1103             ratio_display = ratio_array[elevation_mask, :]
    1104             display_frac = frac_all[elevation_mask, :, :]
    11051111
    11061112            # Compute edges for pcolormesh
    11071113            dt = date_time[1] - date_time[0] if len(date_time) > 1 else 1
    1108             x_edges = np.concatenate([date_time, [date_time[-1] + dt]])
    1109             d_e = np.diff(elev)
    1110             last_e = elev[-1] + (d_e[-1] if len(d_e)>0 else 1)
    1111             y_edges = np.concatenate([elev, [last_e]])
     1114            x_edges = np.concatenate([date_time, [date_time[-1] + (date_time[-1] - date_time[-2])]]) * martian_to_earth
     1115            y_edges = np.concatenate([ref_grid, [ref_grid[-1] + (ref_grid[-1] - ref_grid[-2])]])
    11121116
    11131117            # Plot
    … …  
    11161120                x_edges,
    11171121                y_edges,
    1118                 ratio_display,
     1122                ratio_array,
    11191123                shading='auto',
    11201124                cmap='managua_r',
    11211125                norm=LogNorm(vmin=10**vmin, vmax=10**vmax),
    11221126            )
     1127
     1128            # Custom formatter for millions of Earth years
     1129            def millions_formatter(x, pos):
     1130                return f"{x/1e6:.1f}"
    11231131
    11241132            def format_coord_custom(x_input, y_input):
    … …  
    11351143                i = np.searchsorted(x_edges, x) - 1
    11361144                j = np.searchsorted(y_edges, y) - 1
    1137                 i = np.clip(i, 0, ratio_display.shape[1] - 1)
    1138                 j = np.clip(j, 0, ratio_display.shape[0] - 1)
     1145                i = np.clip(i, 0, ratio_array.shape[1] - 1)
     1146                j = np.clip(j, 0, ratio_array.shape[0] - 1)
    11391147                # get fractions and obliquity
    1140                 fH2O, fCO2, fDust = display_frac[j, i]
    1141                 obl   = np.interp(x, date_time, obliquity)
     1148                fH2O, fCO2, fDust = frac_all[j, i]
     1149                obl   = np.interp(x / martian_to_earth, date_time, obliquity)
    11421150                return f"Time={x:.2f}, Elev={y:.2f}, H2O={fH2O:.4f}, Dust={fDust:.4f}, Obl={obl:.2f}°"
    11431151
    11441152            ax.set_title(f"Dust-to-Ice ratio over time (Grid point {ig+1}, Slope {isl+1})", fontweight='bold')
    1145             ax.set_xlabel('Time (Mars years)')
     1153            ax.xaxis.set_major_formatter(FuncFormatter(millions_formatter))
     1154            ax.set_xlabel('Time (Myr)')
    11461155            ax.set_ylabel('Elevation (m)')
    11471156
    … …  
    11561165            for i in range(len(dh)):
    11571166                ax2.plot(
    1158                     [date_time[i], date_time[i+1]],
     1167                    [date_time[i] * martian_to_earth, date_time[i+1] * martian_to_earth],
    11591168                    [obliquity[i], obliquity[i+1]],
    11601169                    color=color_map[signs[i]],

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

    r3819 r3926  
    232232
    233233
     234def save_dust_profile_ascii(height: np.ndarray, contents: np.ndarray, filename: str = "dust_layering_profile.txt") -> None:
     235    """
     236    Save a simple ASCII file with height corresponding to the middle of each layer (m)
     237    and the dust content (0–1) in the second column.
     238    """
     239    dust_content = contents[2, :]
     240    #max_dust = np.max(dust_content)
     241    #if max_dust <= 0:
     242    #    print("Warning: Dust content is zero everywhere, nothing to normalize.")
     243    #    normalized_dust = dust_content
     244    #else:
     245    #    normalized_dust = dust_content / max_dust
     246
     247    height_mid = 0.5 * (height[:-1] + height[1:])
     248   
     249    if len(dust_content) != len(height_mid):
     250        print("Warning: layer count mismatch; adjusting to shortest length.")
     251        n = min(len(dust_content), len(height_mid))
     252        dust_content = dust_content[:n]
     253        height_mid = height_mid[:n]
     254   
     255    #data = np.column_stack((height_mid, normalized_dust))
     256    data = np.column_stack((height_mid, dust_content))
     257    header = "MidLayer_Height_m  Dust_Content"
     258
     259    np.savetxt(filename, data, fmt="%.6e", header=header)
     260    print(f"Dust profile saved to '{filename}'")
     261
     262
    234263def main():
    235264    if len(sys.argv) > 1:
    … …  
    307336        )
    308337
     338        #save_dust_profile_ascii(height_arr, contents_arr)
     339
    309340        component_labels = ("CO2 Ice", "H2O Ice", "Dust", "Air", "Pore ice")
    310341

  • trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90

    r3851 r3926  
    2020! Physical parameters
    2121logical         :: impose_dust_ratio = .false. ! Impose dust-to-ice ratio instead of dust tendency (see 'dust2ice_ratio')
    22 real            :: dust2ice_ratio    = 0.01    ! Dust-to-ice ratio when ice condenses (1% by default)
    23 real            :: d_dust            = 1.e-3   ! Tendency of dust [kg.m-2.y-1]
    24 real, parameter :: dry_lag_porosity  = 0.4     ! Porosity of dust lag
    25 real, parameter :: wet_lag_porosity  = 0.1     ! Porosity of water ice lag
     22real            :: dust2ice_ratio    = 0.1     ! Dust-to-ice ratio when ice condenses (10% by default)
     23real            :: d_dust            = 5.78e-2 ! Tendency of dust [kg.m-2.y-1] (1.e-9 kg.m-2.s-1 by default)
     24real, parameter :: dry_lag_porosity  = 0.2     ! Porosity of dust lag
     25real, parameter :: wet_lag_porosity  = 0.15    ! Porosity of water ice lag
    2626real, parameter :: regolith_porosity = 0.4     ! Porosity of regolith
    2727real, parameter :: breccia_porosity  = 0.1     ! Porosity of breccia
    28 real, parameter :: co2ice_porosity   = 0.      ! Porosity of CO2 ice
    29 real, parameter :: h2oice_porosity   = 0.      ! Porosity of H2O ice
     28real, parameter :: co2ice_porosity   = 0.05    ! Porosity of CO2 ice
     29real, parameter :: h2oice_porosity   = 0.1     ! Porosity of H2O ice
    3030real, parameter :: rho_dust          = 2500.   ! Density of dust [kg.m-3]
    3131real, parameter :: rho_rock          = 3200.   ! Density of rock [kg.m-3]
    … …  
    603603END SUBROUTINE subsurface_ice_layering
    604604
    605 !=======================================================================
    606605!=======================================================================
    607606! To lift dust in a stratum
Note: See TracChangeset for help on using the changeset viewer.