Index: /trunk/LMDZ.COMMON/libf/evolution/changelog.txt =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3925) +++ /trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3926) @@ -776,2 +776,5 @@ == 03/09/2025 == JBC Managing automatically the slope variables outputted by XIOS according to 'nslope' to avoid subsequent errors. + +== 07/10/2025 == JBC +Updates of files in the deftank + parameter values for the layering algorithm. Index: /trunk/LMDZ.COMMON/libf/evolution/deftank/clean.sh =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/deftank/clean.sh (revision 3925) +++ /trunk/LMDZ.COMMON/libf/evolution/deftank/clean.sh (revision 3926) @@ -27,4 +27,5 @@ rm -f restart1D*.txt rm -f start1D*.txt +rm -f profile_temp.out # Reset the initial starting files to prepare a new simulation if [ -f "starts/startfi.nc" ]; then Index: /trunk/LMDZ.COMMON/libf/evolution/deftank/launchPEM.sh =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/deftank/launchPEM.sh (revision 3925) +++ /trunk/LMDZ.COMMON/libf/evolution/deftank/launchPEM.sh (revision 3926) @@ -79,5 +79,5 @@ errlaunch fi - echo "The relaunch is initialized with a specific previous successful run." + echo "The relaunch is initialized with a previous successful run to be specified." while true; do echo "Do you want to relaunch from a 'PCM' or 'PEM' run?" Index: /trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def (revision 3925) +++ /trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def (revision 3926) @@ -105,12 +105,12 @@ # layering=.false. -# Value of the dust tendency. Default = 1.e-3 kg.m-2.y-1 -# d_dust=1.e-3 +# Value of the dust tendency. Default = 5.78e-2 kg.m-2.y-1 +# d_dust=5.78e-2 # Do you want to impose a dust-to-ice ratio instead of a dust tendency? Default = .false. # impose_dust_ratio=.false. -# If impose_dust_ratio=.true., value of the dust-to-ice ratio. Default = 0.01 -# dust2ice_ratio=0.01 +# If impose_dust_ratio=.true., value of the dust-to-ice ratio. Default = 0.1 +# dust2ice_ratio=0.1 # Some definitions for the physics, in file 'callphys.def' Index: /trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py (revision 3925) +++ /trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py (revision 3926) @@ -13,4 +13,5 @@ from mpl_toolkits.axes_grid1.inset_locator import inset_axes from matplotlib.colors import LinearSegmentedColormap, LogNorm +from matplotlib.ticker import FuncFormatter from scipy.interpolate import interp1d @@ -576,5 +577,4 @@ # Compute RGB stratification over time rgb = np.ones((nz, ntime, 3), dtype=float) - frac_all = np.full((nz, ntime, 3), np.nan, dtype=float) # store fH2O, fCO2, fDust for t in range(ntime): @@ -590,7 +590,10 @@ fCO2 = cCO2 / total fDust = cDust / total - frac_all[:depth, t, :] = np.stack([fH2O, fCO2, fDust], axis=1) - mix = np.outer(fH2O, blue) + np.outer(fCO2, violet) + np.outer(fDust, orange) - rgb[:depth, t, :] = np.clip(mix, 0, 1) + frac_all[:depth, t, 0] = fH2O + frac_all[:depth, t, 1] = fCO2 + frac_all[:depth, t, 2] = fDust + rgb[:depth, t, 0] = fH2O * blue[0] + fCO2 * violet[0] + fDust * orange[0] + rgb[:depth, t, 1] = fH2O * blue[1] + fCO2 * violet[1] + fDust * orange[1] + rgb[:depth, t, 2] = fH2O * blue[2] + fCO2 * violet[2] + fDust * orange[2] # Mask elevation @@ -701,14 +704,7 @@ """ h2o = gridded_data['h2o_ice'] + co2 = gridded_data['co2_ice'] dust = gridded_data['dust'] ngrid, ntime, nslope, nz = h2o.shape - - # Elevation mask - if exclude_sub: - elevation_mask = (ref_grid >= 0.0) - elev = ref_grid[elevation_mask] - else: - elevation_mask = np.ones_like(ref_grid, dtype=bool) - elev = ref_grid # Define custom blue-to-orange colormap @@ -729,6 +725,8 @@ log_ratio_array = np.full((nz, ntime), np.nan, dtype=np.float32) for t in range(ntime): - if ti[t] <= 0: + if t > 0 and ti[t] <= 0: ti[t] = ti[t-1] + elif ti[t] <= 0: + continue zmax = ti[t] if zmax <= 0: @@ -750,9 +748,9 @@ log_ratio_array[:zmax, t] = log_ratio - # Mask and prepare for display ratio_array = 10**log_ratio_array - ratio_display = ratio_array[elevation_mask, :] - x_edges = np.concatenate([date_time, [date_time[-1] + (date_time[-1] - date_time[-2])]]) - y_edges = np.concatenate([elev, [elev[-1] + (elev[-1] - elev[-2])]]) + + # Compute edges for pcolormesh + x_edges = np.concatenate([date_time, [date_time[-1] + (date_time[-1] - date_time[-2])]]) * martian_to_earth + y_edges = np.concatenate([ref_grid, [ref_grid[-1] + (ref_grid[-1] - ref_grid[-2])]]) # Plot @@ -761,10 +759,10 @@ date_time, elev, - ratio_display, + ratio_array, shading='auto', cmap='managua_r', norm=LogNorm(vmin=10**vmin, vmax=10**vmax), ) - attach_format_coord(ax, ratio_display, x_edges, y_edges, is_pcolormesh=True) + attach_format_coord(ax, ratio_array, x_edges, y_edges, is_pcolormesh=True) ax.set_title(f"Dust-to-Ice ratio over time (Grid point {ig+1}, Slope {isl+1})", fontweight='bold') ax.set_xlabel('Time (Mars years)') @@ -966,5 +964,5 @@ ls_perihelion = mars_ls(date_peri_day,0.,eccentricity,year_day) - return dates_yr, obliquity, eccentricity, ls_perihelion + return dates_yr, obliquity, eccentricity, ls_perihelion, martian_to_earth @@ -972,8 +970,9 @@ """ Plot the evolution of obliquity, eccentricity and Ls p coming from simulation data + versus simulated time, plus an additional figure of sin(eccentricity)*Lsp. versus simulated time. """ # Read orbital data - times_yr, obl, ecc, lsp = read_orbital_data_nc(starts_folder, infofile) + times_yr, obl, ecc, lsp, martian_to_earth = read_orbital_data_nc(starts_folder, infofile) fargs = dict(kind='linear', bounds_error=False, fill_value='extrapolate') @@ -982,4 +981,6 @@ lsp_i = interp1d(times_yr, lsp, **fargs)(date_time) + date_time = date_time * martian_to_earth / 1e6 + fig, axes = plt.subplots(3,1, figsize=(8,10), sharex=True) fig.suptitle("Orbital parameters vs simulated time", fontsize=14, fontweight='bold') @@ -996,5 +997,5 @@ axes[2].plot(date_time, lsp_i, 'g-', marker='+') axes[2].set_ylabel("Ls of perihelion (°)") - axes[2].set_xlabel("Time (Mars years)") + axes[2].set_xlabel("Time (Myr)") axes[2].grid(True) @@ -1002,4 +1003,18 @@ outname = os.path.join(output_folder, "orbital_parameters_simu.png") fig.savefig(outname, dpi=150) + + eps_sin_lsp = ecc_i * np.sin(np.radians(lsp_i)) + + fig2, ax2 = plt.subplots(figsize=(8,5)) + fig2.suptitle(r"$\epsilon \times \sin(L_{sp})$", fontweight='bold') + + ax2.plot(date_time, eps_sin_lsp, 'm-', marker='+') + ax2.set_ylabel(r"$\epsilon \cdot \sin(L_{sp})$") + ax2.set_xlabel("Time (Myr)") + ax2.grid(True) + + plt.tight_layout(rect=[0,0,1,0.95]) + outname2 = os.path.join(output_folder, "sin_ecc_times_Lsp.png") + fig2.savefig(outname2, dpi=150) @@ -1025,13 +1040,23 @@ # Read orbital data - times_yr, obl, _, _ = read_orbital_data_nc(starts_folder, infofile) + times_yr, obl, _, _, martian_to_earth = read_orbital_data_nc(starts_folder, infofile) fargs = dict(kind='linear', bounds_error=False, fill_value='extrapolate') obliquity = interp1d(times_yr, obl, **fargs)(date_time) - # Computed total height + # Define custom blue-to-orange colormap + blue = np.array([0, 0, 255], dtype=float) / 255 + orange = np.array([255, 165, 0], dtype=float) / 255 + custom_cmap = LinearSegmentedColormap.from_list('BlueOrange', [blue, orange], N=256) + color_map = { 1: 'green', -1: 'red', 0: 'orange' } + + # Log‑ratio bounds and small epsilon to avoid log(0) + vmin, vmax = -2, 1 + epsilon = 1e-6 + + # Loop over grids and slopes for ig in range(ngrid): for isl in range(nslope): + # Compute total height time series total_height_t = np.zeros(ntime, dtype=float) - for t_idx in range(ntime): h_mat = heights_data[t_idx][isl] @@ -1042,29 +1067,12 @@ total_height_t[t_idx] = np.max(h_valid) - # Compute the per-interval sign of height change - dh = np.diff(total_height_t) - signs = np.sign(dh) - color_map = { 1: 'green', -1: 'red', 0: 'orange' } - - # Elevation mask - if exclude_sub: - elevation_mask = (ref_grid >= 0.0) - elev = ref_grid[elevation_mask] - else: - elevation_mask = np.ones_like(ref_grid, dtype=bool) - elev = ref_grid - - # Define custom blue-to-orange colormap - blue = np.array([0, 0, 255], dtype=float) / 255 - orange = np.array([255, 165, 0], dtype=float) / 255 - custom_cmap = LinearSegmentedColormap.from_list('BlueOrange', [blue, orange], N=256) - - # Log‑ratio bounds and small epsilon to avoid log(0) - vmin, vmax = -2, 1 - epsilon = 1e-6 - - # Loop over grids and slopes - for ig in range(ngrid): - for isl in range(nslope): + # Compute the per-interval sign of height change + if ntime > 1: + dh = np.diff(total_height_t) + signs = np.sign(dh).astype(int) + else: + dh = np.array([], dtype=float) + signs = np.array([], dtype=int) + # Prepare fraction and ratio arrays ti = top_index[ig, :, isl].copy().astype(int) @@ -1072,6 +1080,8 @@ frac_all = np.full((nz, ntime, 3), np.nan, dtype=float) # store fH2O, fCO2, fDust for t in range(ntime): - if ti[t] <= 0: + if t > 0 and ti[t] <= 0: ti[t] = ti[t-1] + elif ti[t] <= 0: + continue zmax = ti[t] if zmax <= 0: @@ -1086,11 +1096,10 @@ fCO2 = cCO2 / total fDust = cDust / total - frac_all[:zmax, t, :] = np.stack([fH2O, fCO2, fDust], axis=1) + frac_all[:zmax, t, 0] = fH2O + frac_all[:zmax, t, 1] = fCO2 + frac_all[:zmax, t, 2] = fDust with np.errstate(divide='ignore', invalid='ignore'): - ratio = np.where( - cH2O > 0, - cDust / cH2O, - 10**(vmax + 1) + ratio = np.where(cH2O > 0, cDust / cH2O, 10**(vmax + 1) ) log_ratio = np.log10(ratio + epsilon) @@ -1099,15 +1108,10 @@ log_ratio_array[:zmax, t] = log_ratio - # Mask elevation ratio_array = 10**log_ratio_array - ratio_display = ratio_array[elevation_mask, :] - display_frac = frac_all[elevation_mask, :, :] # Compute edges for pcolormesh dt = date_time[1] - date_time[0] if len(date_time) > 1 else 1 - x_edges = np.concatenate([date_time, [date_time[-1] + dt]]) - d_e = np.diff(elev) - last_e = elev[-1] + (d_e[-1] if len(d_e)>0 else 1) - y_edges = np.concatenate([elev, [last_e]]) + x_edges = np.concatenate([date_time, [date_time[-1] + (date_time[-1] - date_time[-2])]]) * martian_to_earth + y_edges = np.concatenate([ref_grid, [ref_grid[-1] + (ref_grid[-1] - ref_grid[-2])]]) # Plot @@ -1116,9 +1120,13 @@ x_edges, y_edges, - ratio_display, + ratio_array, shading='auto', cmap='managua_r', norm=LogNorm(vmin=10**vmin, vmax=10**vmax), ) + + # Custom formatter for millions of Earth years + def millions_formatter(x, pos): + return f"{x/1e6:.1f}" def format_coord_custom(x_input, y_input): @@ -1135,13 +1143,14 @@ i = np.searchsorted(x_edges, x) - 1 j = np.searchsorted(y_edges, y) - 1 - i = np.clip(i, 0, ratio_display.shape[1] - 1) - j = np.clip(j, 0, ratio_display.shape[0] - 1) + i = np.clip(i, 0, ratio_array.shape[1] - 1) + j = np.clip(j, 0, ratio_array.shape[0] - 1) # get fractions and obliquity - fH2O, fCO2, fDust = display_frac[j, i] - obl = np.interp(x, date_time, obliquity) + fH2O, fCO2, fDust = frac_all[j, i] + obl = np.interp(x / martian_to_earth, date_time, obliquity) return f"Time={x:.2f}, Elev={y:.2f}, H2O={fH2O:.4f}, Dust={fDust:.4f}, Obl={obl:.2f}°" ax.set_title(f"Dust-to-Ice ratio over time (Grid point {ig+1}, Slope {isl+1})", fontweight='bold') - ax.set_xlabel('Time (Mars years)') + ax.xaxis.set_major_formatter(FuncFormatter(millions_formatter)) + ax.set_xlabel('Time (Myr)') ax.set_ylabel('Elevation (m)') @@ -1156,5 +1165,5 @@ for i in range(len(dh)): ax2.plot( - [date_time[i], date_time[i+1]], + [date_time[i] * martian_to_earth, date_time[i+1] * martian_to_earth], [obliquity[i], obliquity[i+1]], color=color_map[signs[i]], Index: /trunk/LMDZ.COMMON/libf/evolution/deftank/visu_layering.py =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/deftank/visu_layering.py (revision 3925) +++ /trunk/LMDZ.COMMON/libf/evolution/deftank/visu_layering.py (revision 3926) @@ -232,4 +232,33 @@ +def save_dust_profile_ascii(height: np.ndarray, contents: np.ndarray, filename: str = "dust_layering_profile.txt") -> None: + """ + Save a simple ASCII file with height corresponding to the middle of each layer (m) + and the dust content (0–1) in the second column. + """ + dust_content = contents[2, :] + #max_dust = np.max(dust_content) + #if max_dust <= 0: + # print("Warning: Dust content is zero everywhere, nothing to normalize.") + # normalized_dust = dust_content + #else: + # normalized_dust = dust_content / max_dust + + height_mid = 0.5 * (height[:-1] + height[1:]) + + if len(dust_content) != len(height_mid): + print("Warning: layer count mismatch; adjusting to shortest length.") + n = min(len(dust_content), len(height_mid)) + dust_content = dust_content[:n] + height_mid = height_mid[:n] + + #data = np.column_stack((height_mid, normalized_dust)) + data = np.column_stack((height_mid, dust_content)) + header = "MidLayer_Height_m Dust_Content" + + np.savetxt(filename, data, fmt="%.6e", header=header) + print(f"Dust profile saved to '{filename}'") + + def main(): if len(sys.argv) > 1: @@ -307,4 +336,6 @@ ) + #save_dust_profile_ascii(height_arr, contents_arr) + component_labels = ("CO2 Ice", "H2O Ice", "Dust", "Air", "Pore ice") Index: /trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90 (revision 3925) +++ /trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90 (revision 3926) @@ -20,12 +20,12 @@ ! Physical parameters logical :: impose_dust_ratio = .false. ! Impose dust-to-ice ratio instead of dust tendency (see 'dust2ice_ratio') -real :: dust2ice_ratio = 0.01 ! Dust-to-ice ratio when ice condenses (1% by default) -real :: d_dust = 1.e-3 ! Tendency of dust [kg.m-2.y-1] -real, parameter :: dry_lag_porosity = 0.4 ! Porosity of dust lag -real, parameter :: wet_lag_porosity = 0.1 ! Porosity of water ice lag +real :: dust2ice_ratio = 0.1 ! Dust-to-ice ratio when ice condenses (10% by default) +real :: d_dust = 5.78e-2 ! Tendency of dust [kg.m-2.y-1] (1.e-9 kg.m-2.s-1 by default) +real, parameter :: dry_lag_porosity = 0.2 ! Porosity of dust lag +real, parameter :: wet_lag_porosity = 0.15 ! Porosity of water ice lag real, parameter :: regolith_porosity = 0.4 ! Porosity of regolith real, parameter :: breccia_porosity = 0.1 ! Porosity of breccia -real, parameter :: co2ice_porosity = 0. ! Porosity of CO2 ice -real, parameter :: h2oice_porosity = 0. ! Porosity of H2O ice +real, parameter :: co2ice_porosity = 0.05 ! Porosity of CO2 ice +real, parameter :: h2oice_porosity = 0.1 ! Porosity of H2O ice real, parameter :: rho_dust = 2500. ! Density of dust [kg.m-3] real, parameter :: rho_rock = 3200. ! Density of rock [kg.m-3] @@ -603,5 +603,4 @@ END SUBROUTINE subsurface_ice_layering -!======================================================================= !======================================================================= ! To lift dust in a stratum