Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3842) +++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3850) @@ -733,2 +733,9 @@ - Correction of the incorporation of the sub-surface sublimation tendency into the overall tendency to be more robust. - Revision of the layering initialization in the case where there is no "startpem.nc" file. + +== 15/07/2025 == JBC +- Handling the case where the new date is positive (present-day climate) in "recomp_orb_param_mod.F90", especially because of rounding errors. +- Correction in "visu_evol_layering.py" to get the proper extent of the stratication + replacing 'imshow' by 'pcolormesh' to better process the data. +- Addition of a graphic to plot the changes in obliquity directly onto the evolution of the stratification over time for better analysis. Obliquity curve is colored according to stratification gain/loss. +- Correction of numbering in "lib_launchPEM.sh" when requesting a relaunch. +- Update of "PEMrun.job" due to new parsing of arguments introduced in r3836/r3837. Index: trunk/LMDZ.COMMON/libf/evolution/deftank/PEMrun.job =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/deftank/PEMrun.job (revision 3842) +++ trunk/LMDZ.COMMON/libf/evolution/deftank/PEMrun.job (revision 3850) @@ -26,5 +26,5 @@ # Argument for the PEM execution (for SLURM: "$SLURM_JOB_ID" | for PBD/TORQUE: "$PBS_JOBID" | "" when the script is not run as a job): -arg_pem="$SLURM_JOB_ID" +arg_pem="--jobid $SLURM_JOB_ID" ######################################################################## Index: trunk/LMDZ.COMMON/libf/evolution/deftank/lib_launchPEM.sh =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/deftank/lib_launchPEM.sh (revision 3842) +++ trunk/LMDZ.COMMON/libf/evolution/deftank/lib_launchPEM.sh (revision 3850) @@ -172,5 +172,5 @@ fi if [ $(echo "$i_myear < $n_myear" | bc -l) -eq 1 ]; then - echo "Run \"PCM $iPCM\" ($ii/$3)..." + echo "Run \"PCM $iPCM\" ($ii/$3)" if [ $1 -eq 0 ]; then # Mode: processing scripts sed -i "s/^k=-\?[0-9]\+$/k=$(echo "$ii - $3 + 2" | bc)/" PCMrun.job @@ -199,5 +199,5 @@ for ((i = $ii; i <= $3; i++)); do if [ $(echo "$i_myear < $n_myear" | bc -l) -eq 1 ]; then - echo "Run \"PCM $iPCM\" ($i/$3)..." + echo "Run \"PCM $iPCM\" ($i/$3)" if [ $1 -eq 0 ]; then # Mode: processing scripts sed -i "s/^k=-\?[0-9]\+$/k=$(echo "$i - $3 + 2" | bc)/" PCMrun.job @@ -366,5 +366,5 @@ if [ $il -eq $(($nPCM - 1)) ]; then # Second to last PCM run cp diags/data2reshape${irelaunch}.nc data2reshape_Y1.nc - cyclelaunch $1 $2 $nPCM $il + cyclelaunch $1 $2 $nPCM $(($il + 1)) elif [ $il -eq $nPCM ]; then # Last PCM run so the next job is a PEM run cp diags/data2reshape$(($irelaunch - 1)).nc data2reshape_Y1.nc @@ -372,5 +372,5 @@ submitPEM $1 else - cyclelaunch $1 $2 $nPCM $il + cyclelaunch $1 $2 $nPCM $(($il + 1)) fi fi Index: trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py (revision 3842) +++ trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py (revision 3850) @@ -525,15 +525,18 @@ Includes a triangular legend showing the mix proportions. """ - # Define constant colors - violet = np.array([255, 0, 255], dtype=float) / 255 - blue = np.array([0, 0, 255], dtype=float) / 255 - orange = np.array([255, 165, 0], dtype=float) / 255 - - # Prepare elevation mask - mask_elev = (ref_grid >= 0.0) if exclude_sub else np.ones_like(ref_grid, dtype=bool) - elev = ref_grid[mask_elev] - - # Generate legend image once + violet = np.array([255, 0, 255], dtype=float) / 255 + blue = np.array([ 0, 0, 255], dtype=float) / 255 + orange = np.array([255, 165, 0], dtype=float) / 255 + + # Elevation mask and array + 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 + + # Pre-compute legend triangle res = 300 u = np.linspace(0, 1, res) @@ -544,5 +547,4 @@ W_bary = 1 - U_bary - V_bary mask_triangle = (U_bary >= 0) & (V_bary >= 0) & (W_bary >= 0) - legend_rgb = ( U_bary[..., None] * violet @@ -555,5 +557,5 @@ legend_rgba[..., 3] = mask_triangle.astype(float) - # Loop over grid and slope + # Extract data arrays h2o = gridded_data['h2o_ice'] co2 = gridded_data['co2_ice'] @@ -561,15 +563,26 @@ ngrid, ntime, nslope, nz = h2o.shape + # Fill missing depths + ti = top_index.copy().astype(int) + for ig in range(ngrid): + for isl in range(nslope): + for t in range(1, ntime): + if ti[ig, t, isl] <= 0: + ti[ig, t, isl] = ti[ig, t-1, isl] + + # Loop over grid and slope for ig in range(ngrid): for isl in range(nslope): # Compute RGB stratification over time rgb = np.ones((nz, ntime, 3), dtype=float) + + frac_all = np.zeros((nz, ntime, 3), dtype=float) # store fH2O, fCO2, fDust for t in range(ntime): - mask_z = np.arange(nz) < top_index[ig, t, isl] - if not mask_z.any(): + depth = ti[ig, t, isl] + if depth <= 0: continue - cH2O = np.clip(h2o[ig, t, isl, mask_z], 0, None) - cCO2 = np.clip(co2[ig, t, isl, mask_z], 0, None) - cDust = np.clip(dust[ig, t, isl, mask_z], 0, None) + cH2O = np.clip(h2o[ig, t, isl, :depth], 0, None) + cCO2 = np.clip(co2[ig, t, isl, :depth], 0, None) + cDust = np.clip(dust[ig, t, isl, :depth], 0, None) total = cH2O + cCO2 + cDust total[total == 0] = 1.0 @@ -577,31 +590,50 @@ fCO2 = cCO2 / total fDust = cDust / total - mix = ( - np.outer(fH2O, blue) - + np.outer(fCO2, violet) - + np.outer(fDust, orange) - ) - mix = np.clip(mix, 0.0, 1.0) - rgb[mask_z, t, :] = mix - - display_rgb = rgb[mask_elev, :, :] + 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) + + # Mask elevation + display_rgb = rgb[elevation_mask, :, :] + display_frac = frac_all[elevation_mask, :, :] + + display_rgb = rgb[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]]) # Create figure with legend fig, (ax_main, ax_leg) = plt.subplots( - 1, 2, figsize=(12, 5), dpi=200, + 1, 2, figsize=(8, 4), dpi=200, gridspec_kw={'width_ratios': [5, 1]} ) # Main stratification panel - ax_main.imshow( + mesh = ax_main.pcolormesh( + x_edges, + y_edges, display_rgb, - aspect='auto', - extent=[date_time[0], date_time[-1], elev.min(), elev.max()], - interpolation='nearest', - origin='lower' + shading='auto', + edgecolors='none' ) - x_centers = np.linspace(date_time[0], date_time[-1], display_rgb.shape[1]) - y_centers = np.linspace(elev.min(), elev.max(), display_rgb.shape[0]) - attach_format_coord(ax_main, display_rgb, x_centers, y_centers, is_pcolormesh=False) + + # Custom coordinate formatter: show time, elevation, and mixture fractions + def main_format(x, y): + # check bounds + if x < x_edges[0] or x > x_edges[-1] or y < y_edges[0] or y > y_edges[-1]: + return '' + # locate cell + i = np.searchsorted(x_edges, x) - 1 + j = np.searchsorted(y_edges, y) - 1 + i = np.clip(i, 0, display_rgb.shape[1]-1) + j = np.clip(j, 0, display_rgb.shape[0]-1) + # get fractions + fH2O, fCO2, fDust = display_frac[j, i] + return f"Time={x:.2f}, Elev={y:.2f}, H2O={fH2O:.2f}, CO2={fCO2:.2f}, Dust={fDust:.2f}" + ax_main.format_coord = main_format ax_main.set_facecolor('white') ax_main.set_title(f"Ternary mix over time (Grid point {ig+1}, Slope {isl+1})", fontweight='bold') @@ -609,24 +641,37 @@ ax_main.set_ylabel("Elevation (m)") - # Legend panel - ax_leg.imshow( + # Legend panel using proper edges + u_edges = np.linspace(0, 1, res+1) + v_edges = np.linspace(0, np.sqrt(3)/2, res+1) + ax_leg.pcolormesh( + u_edges, + v_edges, legend_rgba, - extent=[0, 1, 0, np.sqrt(3)/2], - origin='lower', - interpolation='nearest' + shading='auto', + edgecolors='none' ) - attach_format_coord(ax_leg, legend_rgba, np.array([0, 1]), np.array([0, np.sqrt(3)/2]), is_pcolormesh=False) - - # Draw triangle border + ax_leg.set_aspect('equal') + + # Custom coordinate formatter for legend: show barycentric fractions + def legend_format(x, y): + # compute barycentric coords from cartesian (x,y) + V = 2 * y / np.sqrt(3) + U = x - 0.5 * V + W = 1 - U - V + if U >= 0 and V >= 0 and W >= 0: + return f"H2O: {W:.2f}, Dust: {V:.2f}, CO2: {U:.2f}" + else: + return '' + ax_leg.format_coord = legend_format + + # Draw triangle border and gridlines triangle = np.array([[0, 0], [1, 0], [0.5, np.sqrt(3)/2], [0, 0]]) - ax_leg.plot(triangle[:, 0], triangle[:, 1], 'k-', linewidth=1) - - # Dashed gridlines + ax_leg.plot(triangle[:, 0], triangle[:, 1], 'k-', linewidth=1, clip_on=False, zorder=10) ticks = np.linspace(0.25, 0.75, 3) for f in ticks: - ax_leg.plot([1 - f, 0.5 * (1 - f)], [0, (1 - f)*np.sqrt(3)/2], '--', color='k', linewidth=0.5) - ax_leg.plot([f, f + 0.5 * (1 - f)], [0, (1 - f)*np.sqrt(3)/2], '--', color='k', linewidth=0.5) + ax_leg.plot([1 - f, 0.5 * (1 - f)], [0, (1 - f)*np.sqrt(3)/2], '--', color='k', linewidth=0.5, clip_on=False, zorder=9) + ax_leg.plot([f, f + 0.5 * (1 - f)], [0, (1 - f)*np.sqrt(3)/2], '--', color='k', linewidth=0.5, clip_on=False, zorder=9) y = (np.sqrt(3)/2) * f - ax_leg.plot([0.5 * f, 1 - 0.5 * f], [y, y], '--', color='k', linewidth=0.5) + ax_leg.plot([0.5 * f, 1 - 0.5 * f], [y, y], '--', color='k', linewidth=0.5, clip_on=False, zorder=9) # Legend labels @@ -636,7 +681,6 @@ ax_leg.axis('off') + # Save figure plt.tight_layout() - - # Save figure fname = os.path.join(output_folder, f"layering_rgb_evolution_ig{ig+1}_is{isl+1}.png") fig.savefig(fname, dpi=150, bbox_inches='tight') @@ -678,13 +722,17 @@ vmin, vmax = -2, 1 epsilon = 1e-6 - + # Loop over grids and slopes for ig in range(ngrid): for isl in range(nslope): + ti = top_index[ig, :, isl].copy().astype(int) + for t in range(1, ntime): + if ti[t] <= 0: + ti[t] = ti[t-1] + + # Compute log10(dust/ice) profile at each time step log_ratio_array = np.full((nz, ntime), np.nan, dtype=np.float32) - - # Compute log10(dust/ice) profile at each time step for t in range(ntime): - zmax = top_index[ig, t, isl] + zmax = ti[t] if zmax <= 0: continue @@ -710,17 +758,17 @@ # Plot fig, ax = plt.subplots(figsize=(8, 6), dpi=150) + im = ax.pcolormesh( + date_time, + elev, + ratio_display, + shading='auto', + cmap='managua_r', + norm=LogNorm(vmin=10**vmin, vmax=10**vmax), + ) + x_edges = np.concatenate([date_time, [date_time[-1] + (date_time[-1]-date_time[-2])]]) + attach_format_coord(ax, ratio_display, x_edges, np.concatenate([elev, [elev[-1] + (elev[-1]-elev[-2])]]), is_pcolormesh=True) ax.set_title(f"Dust-to-Ice ratio over time (Grid point {ig+1}, Slope {isl+1})", fontweight='bold') - im = ax.imshow( - ratio_display, - aspect='auto', - extent=[date_time[0], date_time[-1], elev.min(), elev.max()], - origin='lower', - interpolation='nearest', - cmap='managua_r', - norm=LogNorm(vmin=10**vmin, vmax=10**vmax) - ) - x_centers = np.linspace(date_time[0], date_time[-1], ratio_display.shape[1]) - y_centers = np.linspace(elev.min(), elev.max(), ratio_display.shape[0]) - attach_format_coord(ax, ratio_display, x_centers, y_centers, is_pcolormesh=False) + ax.set_xlabel('Time (Mars years)') + ax.set_ylabel('Elevation (m)') # Add colorbar with simplified ratio labels @@ -827,22 +875,279 @@ # Plot fig, axes = plt.subplots(3, 1, figsize=(8, 10), sharex=True) - fig.suptitle("Orbital Parameters vs Simulated Time", fontsize=14, fontweight='bold') - - axes[0].plot(date_time, obl_interp, 'r+', linestyle='-') + fig.suptitle("Orbital parameters vs simulated time", fontsize=14, fontweight='bold') + + axes[0].plot(date_time, obl_interp, 'r-', marker='+') axes[0].set_ylabel("Obliquity (°)") axes[0].grid(True) - axes[1].plot(date_time, ecc_interp, 'b+', linestyle='-') + axes[1].plot(date_time, ecc_interp, 'b-', marker='+') axes[1].set_ylabel("Eccentricity") axes[1].grid(True) - axes[2].plot(date_time, lsp_interp, 'g+', linestyle='-') - axes[2].set_ylabel("Ls p (°)") + axes[2].plot(date_time, lsp_interp, 'g-', marker='+') + axes[2].set_ylabel("Ls of perihelion (°)") axes[2].set_xlabel("Time (Mars years)") axes[2].grid(True) plt.tight_layout(rect=[0, 0, 1, 0.96]) - fname = os.path.join(output_folder, "orbital_parameters.png") + fname = os.path.join(output_folder, "orbital_parameters_laskar.png") fig.savefig(fname, dpi=150) + + +def mars_ls(pday, peri_day, e_elips, year_day, lsperi=0.0): + """ + Compute solar longitude (Ls) in radians for a given Mars date array 'pday'. + Returns Ls in degrees [0, 360). + """ + zz = (pday - peri_day) / year_day + zanom = 2 * np.pi * (zz - np.round(zz)) + xref = np.abs(zanom) + + # Solve Kepler's equation via Newton–Raphson + zx0 = xref + e_elips * np.sin(xref) + for _ in range(10): + f = zx0 - e_elips * np.sin(zx0) - xref + fp = 1 - e_elips * np.cos(zx0) + dz = -f / fp + zx0 += dz + if np.all(np.abs(dz) <= 1e-7): + break + + zx0 = np.where(zanom < 0, -zx0, zx0) + zteta = 2 * np.arctan( + np.sqrt((1 + e_elips) / (1 - e_elips)) * np.tan(zx0 / 2) + ) + psollong = np.mod(zteta + lsperi, 2 * np.pi) + + return np.degrees(psollong) + + +def read_orbital_data_nc(starts_folder, infofile=None): + """ + Read orbital parameters from restartfi_postPEM*.nc files in starts_folder. + """ + if not os.path.isdir(starts_folder): + raise ValueError(f"Invalid starts_folder '{starts_folder}': not a directory.") + + # Read simulation time mapping if provided + if infofile: + dates_yr, martian_to_earth = read_infofile(infofile) + else: + dates_yr = None + + pattern = os.path.join(starts_folder, "restartfi_postPEM*.nc") + files = glob(pattern) + if not files: + raise FileNotFoundError(f"No NetCDF restart files found matching {pattern}") + + def extract_number(path): + name = os.path.basename(path) + prefix = 'restartfi_postPEM' + if name.startswith(prefix) and name.endswith('.nc'): + num_str = name[len(prefix):-3] + if num_str.isdigit(): + return int(num_str) + return float('inf') + + files = sorted(files, key=extract_number) + + all_year_day, all_peri, all_aphe, all_date_peri, all_obl = [], [], [], [], [] + for nc_path in files: + with Dataset(nc_path, 'r') as nc: + ctrl = nc.variables['controle'][:] + all_year_day.append(ctrl[13]) + all_peri.append(ctrl[14]) + all_aphe.append(ctrl[15]) + all_date_peri.append(ctrl[16]) + all_obl.append(ctrl[17]) + + year_day = np.array(all_year_day) + perihelion = np.array(all_peri) + aphelion = np.array(all_aphe) + date_peri_day = np.array(all_date_peri) + obliquity = np.array(all_obl) + + eccentricity = (aphelion - perihelion) / (aphelion + perihelion) + ls_perihelion = mars_ls( + date_peri_day, + date_peri_day, + eccentricity, + year_day + ) + + return dates_yr, obliquity, eccentricity, ls_perihelion + + +def plot_orbital_parameters_nc(starts_folder, infofile, date_time, output_folder="."): + """ + Plot the evolution of obliquity, eccentricity and Ls p coming from simulation data + versus simulated time. + """ + # Read orbital data + times_yr, obl, ecc, lsp = read_orbital_data_nc(starts_folder, infofile) + + fargs = dict(kind='linear', bounds_error=False, fill_value='extrapolate') + obl_i = interp1d(times_yr, obl, **fargs)(date_time) + ecc_i = interp1d(times_yr, ecc, **fargs)(date_time) + lsp_i = interp1d(times_yr, lsp, **fargs)(date_time) + + fig, axes = plt.subplots(3,1, figsize=(8,10), sharex=True) + fig.suptitle("Orbital parameters vs simulated time", fontsize=14, fontweight='bold') + + # Plot + axes[0].plot(date_time, obl_i, 'r-', marker='+') + axes[0].set_ylabel("Obliquity (°)") + axes[0].grid(True) + + axes[1].plot(date_time, ecc_i, 'b-', marker='+') + axes[1].set_ylabel("Eccentricity") + axes[1].grid(True) + + 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].grid(True) + + plt.tight_layout(rect=[0,0,1,0.96]) + outname = os.path.join(output_folder, "orbital_parameters_simu.png") + fig.savefig(outname, dpi=150) + + +def plot_dust_to_ice_ratio_with_obliquity( + starts_folder, + infofile, + gridded_data, + ref_grid, + top_index, + heights_data, + date_time, + exclude_sub=False, + output_folder="." +): + """ + Plot the dust-to-ice ratio over time as a heatmap, and overlay the evolution of + obliquity on a secondary y-axis. + """ + h2o = gridded_data['h2o_ice'] + dust = gridded_data['dust'] + ngrid, ntime, nslope, nz = h2o.shape + + # Read orbital data + times_yr, obl, _, _ = 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 + for ig in range(ngrid): + for isl in range(nslope): + total_height_t = np.zeros(ntime, dtype=float) + + for t_idx in range(ntime): + h_mat = heights_data[t_idx][isl] + raw_h = h_mat[ig, :] + valid_mask = (~np.isnan(raw_h)) & (raw_h != 0.0) + if np.any(valid_mask): + h_valid = raw_h[valid_mask] + 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 + + # Custom colormap: blue (ice) to orange (dust) + blue = np.array([0, 0, 255]) / 255 + orange = np.array([255, 165, 0]) / 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): + ti = top_index[ig, :, isl].copy().astype(int) + for t in range(1, ntime): + if ti[t] <= 0: + ti[t] = ti[t-1] + + # Compute log10(dust/ice) profile at each time step + log_ratio_array = np.full((nz, ntime), np.nan, dtype=np.float32) + for t in range(ntime): + zmax = ti[t] + if zmax <= 0: + continue + + h2o_profile = np.clip(h2o[ig, t, isl, :zmax], 0, None) + dust_profile = np.clip(dust[ig, t, isl, :zmax], 0, None) + + with np.errstate(divide='ignore', invalid='ignore'): + ratio_profile = np.where( + h2o_profile > 0, + dust_profile / h2o_profile, + 10**(vmax + 1) + ) + log_ratio = np.log10(ratio_profile + epsilon) + log_ratio = np.clip(log_ratio, vmin, vmax) + + log_ratio_array[:zmax, t] = log_ratio + + # Convert back to linear ratio and mask + ratio_array = 10**log_ratio_array + ratio_display = ratio_array[elevation_mask, :] + + # Plot + fig, ax = plt.subplots(figsize=(8, 6), dpi=150) + im = ax.pcolormesh( + date_time, + elev, + ratio_display, + shading='auto', + cmap='managua_r', + norm=LogNorm(vmin=10**vmin, vmax=10**vmax), + ) + x_edges = np.concatenate([date_time, [date_time[-1] + (date_time[-1]-date_time[-2])]]) + attach_format_coord(ax, ratio_display, x_edges, np.concatenate([elev, [elev[-1] + (elev[-1]-elev[-2])]]), 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)') + ax.set_ylabel('Elevation (m)') + + # Add colorbar + cbar = fig.colorbar(im, ax=ax, orientation='vertical', pad=0.15) + cbar.set_label('Dust / H₂O ice (ratio)') + cbar.set_ticks([1e-2, 1e-1, 1, 1e1]) + cbar.set_ticklabels(['1:100', '1:10', '1:1', '10:1']) + + # Overlay obliquity on secondary y-axis + ax2 = ax.twinx() + for i in range(len(dh)): + ax2.plot( + [date_time[i], date_time[i+1]], + [obliquity[i], obliquity[i+1]], + color=color_map[signs[i]], + marker='+', + linewidth=1.5 + ) + ax2.set_ylabel('Obliquity (°)') + ax2.tick_params(axis='y') + ax2.grid(False) + + # Save + os.makedirs(output_folder, exist_ok=True) + outname = os.path.join( + output_folder, + f'dust_ice_obliquity_grid{ig+1}_slope{isl+1}.png' + ) + plt.tight_layout() + fig.savefig(outname, dpi=150) @@ -913,12 +1218,18 @@ exclude_sub=exclude_sub, output_folder="." ) - plot_dust_to_ice_ratio_over_time( + #plot_dust_to_ice_ratio_over_time( + # gridded_data, ref_grid, top_index, heights_data, date_time, + # exclude_sub=exclude_sub, output_folder="." + #) + plot_dust_to_ice_ratio_with_obliquity( + folder_path, infofile, gridded_data, ref_grid, top_index, heights_data, date_time, exclude_sub=exclude_sub, output_folder="." ) - plot_strata_count_and_total_height(heights_data, date_time, output_folder=".") + #plot_strata_count_and_total_height(heights_data, date_time, output_folder=".") # 14) Plot orbital parameters - plot_orbital_parameters(infofile, orbfile, date_time, output_folder=".") + #plot_orbital_parameters(infofile, orbfile, date_time, output_folder=".") + plot_orbital_parameters_nc(folder_path, infofile, date_time, output_folder=".") # 15) Show all figures Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3842) +++ trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3850) @@ -69,6 +69,6 @@ use writediagpem_mod, only: writediagpem, writediagsoilpem use co2condens_mod, only: CO2cond_ps -use layering_mod, only: layering, del_layering, make_layering, layering_algo, subsurface_ice_layering, & - ptrarray, stratum, get_nb_str_max, nb_str_max, is_dust_lag, is_co2ice_str, is_h2oice_str, print_layering +use layering_mod, only: layering, del_layering, make_layering, layering_algo, subsurface_ice_layering, & + ptrarray, stratum, get_nb_str_max, nb_str_max, is_dust_lag, is_co2ice_str, is_h2oice_str use dyn_ss_ice_m_mod, only: dyn_ss_ice_m use parse_args_mod, only: parse_args Index: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 (revision 3842) +++ trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 (revision 3850) @@ -393,5 +393,4 @@ do islope = 1,nslope call ini_layering(layerings_map(ig,islope)) - print*, 'coucou', 1.05*ini_huge_h2oice/rho_h2oice, layerings_map(ig,islope)%top%top_elevation call add_stratum(layerings_map(ig,islope),1.05*ini_huge_h2oice/rho_h2oice,0.,ini_huge_h2oice/rho_h2oice,0.05*ini_huge_h2oice/rho_h2oice,0.,0.) enddo Index: trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90 (revision 3842) +++ trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90 (revision 3850) @@ -31,5 +31,5 @@ ! Input Variables !-------------------------------------------------------- -real, intent(in) :: i_myear ! Number of simulated Martian years +real, intent(in) :: i_myear ! Number of simulated Martian years real, intent(in) :: i_myear_leg ! Number of iterations of the PEM @@ -62,4 +62,6 @@ Lsp = lsperi*180./pi ! We convert in degree +print*, year_bp_ini,i_myear,i_myear_leg + write(*,*) 'recomp_orb_param, Old year =', Year - i_myear_leg write(*,*) 'recomp_orb_param, Old obl. =', obliquit @@ -69,41 +71,52 @@ found_year = .false. -do ilask = last_ilask,2,-1 - xa = yearlask(ilask) - xb = yearlask(ilask - 1) - if (xa <= Year .and. Year < xb) then - ! Obliquity - if (var_obl) then - ya = obllask(ilask) - yb = obllask(ilask - 1) - obliquit = (Year - xa)*(yb - ya)/(xb - xa) + ya +if (Year < 0.) then + do ilask = last_ilask,2,-1 + xa = yearlask(ilask) + xb = yearlask(ilask - 1) + if (xa <= Year .and. Year < xb) then + ! Obliquity + if (var_obl) then + ya = obllask(ilask) + yb = obllask(ilask - 1) + obliquit = (Year - xa)*(yb - ya)/(xb - xa) + ya + endif + + ! Eccentricity + if (var_ecc) then + ya = ecclask(ilask) + yb = ecclask(ilask - 1) + e_elips = (Year - xa)*(yb - ya)/(xb - xa) + ya + endif + + ! Lsp + if (var_lsp) then + ya = lsplask(ilask) + yb = lsplask(ilask - 1) + if (abs(yb - ya) > 300.) then ! If modulo is "crossed" through the interval + if (ya < yb) then ! Lsp should be decreasing + ya = ya + 360. + else ! Lsp should be increasing + yb = yb + 360. + endif + endif + Lsp = modulo((Year - xa)*(yb - ya)/(xb - xa) + ya,360.) + endif + found_year = .true. + exit ! The loop is left as soon as the right interval is found endif - - ! Eccentricity - if (var_ecc) then - ya = ecclask(ilask) - yb = ecclask(ilask - 1) - e_elips = (Year - xa)*(yb - ya)/(xb - xa) + ya - endif - - ! Lsp - if (var_lsp) then - ya = lsplask(ilask) - yb = lsplask(ilask - 1) - if (abs(yb - ya) > 300.) then ! If modulo is "crossed" through the interval - if (ya < yb) then ! Lsp should be decreasing - ya = ya + 360. - else ! Lsp should be increasing - yb = yb + 360. - endif - endif - Lsp = modulo((Year - xa)*(yb - ya)/(xb - xa) + ya,360.) - endif - found_year = .true. - exit ! The loop is left as soon as the right interval is found - endif -enddo + enddo +else if (Year >= 0 .and. Year < 100.) then ! For present orbital characteristics + if (var_obl) obliquit = 25.19 + if (var_ecc) e_elips = 0.0934 + if (var_lsp) Lsp = 251. + found_year = .true. +endif if (.not. found_year) error stop 'The new year could not be found in the file "obl_ecc_lsp.asc".' + +write(*,*) 'recomp_orb_param, New obl. =', obliquit +write(*,*) 'recomp_orb_param, New ecc. =', e_elips +write(*,*) 'recomp_orb_param, New Lsp =', Lsp halfaxe = 227.94 @@ -118,8 +131,4 @@ #endif -write(*,*) 'recomp_orb_param, New obl. =', obliquit -write(*,*) 'recomp_orb_param, New ecc. =', e_elips -write(*,*) 'recomp_orb_param, New Lsp =', Lsp - END SUBROUTINE recomp_orb_param