Changeset 3850


Ignore:
Timestamp:
Jul 15, 2025, 2:43:41 PM (21 hours ago)
Author:
jbclement
Message:

PEM:

  • 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.

JBC

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

Legend:

Unmodified
Added
Removed

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

    r3842 r3850  
    733733- Correction of the incorporation of the sub-surface sublimation tendency into the overall tendency to be more robust.
    734734- Revision of the layering initialization in the case where there is no "startpem.nc" file.
     735
     736== 15/07/2025 == JBC
     737- Handling the case where the new date is positive (present-day climate) in "recomp_orb_param_mod.F90", especially because of rounding errors.
     738- Correction in "visu_evol_layering.py" to get the proper extent of the stratication + replacing 'imshow' by 'pcolormesh' to better process the data.
     739- 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.
     740- Correction of numbering in "lib_launchPEM.sh" when requesting a relaunch.
     741- Update of "PEMrun.job" due to new parsing of arguments introduced in r3836/r3837.

  • trunk/LMDZ.COMMON/libf/evolution/deftank/PEMrun.job

    r3667 r3850  
    2626
    2727# Argument for the PEM execution (for SLURM: "$SLURM_JOB_ID" | for PBD/TORQUE: "$PBS_JOBID" | "" when the script is not run as a job):
    28 arg_pem="$SLURM_JOB_ID"
     28arg_pem="--jobid $SLURM_JOB_ID"
    2929########################################################################
    3030

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

    r3840 r3850  
    172172    fi
    173173    if [ $(echo "$i_myear < $n_myear" | bc -l) -eq 1 ]; then
    174         echo "Run \"PCM $iPCM\" ($ii/$3)..."
     174        echo "Run \"PCM $iPCM\" ($ii/$3)"
    175175        if [ $1 -eq 0 ]; then # Mode: processing scripts
    176176            sed -i "s/^k=-\?[0-9]\+$/k=$(echo "$ii - $3 + 2" | bc)/" PCMrun.job
    … …  
    199199    for ((i = $ii; i <= $3; i++)); do
    200200        if [ $(echo "$i_myear < $n_myear" | bc -l) -eq 1 ]; then
    201             echo "Run \"PCM $iPCM\" ($i/$3)..."
     201            echo "Run \"PCM $iPCM\" ($i/$3)"
    202202            if [ $1 -eq 0 ]; then # Mode: processing scripts
    203203                sed -i "s/^k=-\?[0-9]\+$/k=$(echo "$i - $3 + 2" | bc)/" PCMrun.job
    … …  
    366366        if [ $il -eq $(($nPCM - 1)) ]; then # Second to last PCM run
    367367            cp diags/data2reshape${irelaunch}.nc data2reshape_Y1.nc
    368             cyclelaunch $1 $2 $nPCM $il
     368            cyclelaunch $1 $2 $nPCM $(($il + 1))
    369369        elif [ $il -eq $nPCM ]; then # Last PCM run so the next job is a PEM run
    370370            cp diags/data2reshape$(($irelaunch - 1)).nc data2reshape_Y1.nc
    … …  
    372372            submitPEM $1
    373373        else
    374             cyclelaunch $1 $2 $nPCM $il
     374            cyclelaunch $1 $2 $nPCM $(($il + 1))
    375375        fi
    376376    fi

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

    r3840 r3850  
    525525    Includes a triangular legend showing the mix proportions.
    526526    """
    527 
    528527    # Define constant colors
    529     violet = np.array([255, 0, 255], dtype=float) / 255
    530     blue   = np.array([0, 0, 255], dtype=float) / 255
    531     orange = np.array([255, 165, 0], dtype=float) / 255
    532 
    533     # Prepare elevation mask
    534     mask_elev = (ref_grid >= 0.0) if exclude_sub else np.ones_like(ref_grid, dtype=bool)
    535     elev = ref_grid[mask_elev]
    536 
    537     # Generate legend image once
     528    violet = np.array([255,   0, 255], dtype=float) / 255
     529    blue   = np.array([  0,   0, 255], dtype=float) / 255
     530    orange = np.array([255, 165,   0], dtype=float) / 255
     531
     532    # Elevation mask and array
     533    if exclude_sub:
     534        elevation_mask = (ref_grid >= 0.0)
     535        elev = ref_grid[elevation_mask]
     536    else:
     537        elevation_mask = np.ones_like(ref_grid, dtype=bool)
     538        elev = ref_grid
     539
     540    # Pre-compute legend triangle
    538541    res = 300
    539542    u = np.linspace(0, 1, res)
    … …  
    544547    W_bary = 1 - U_bary - V_bary
    545548    mask_triangle = (U_bary >= 0) & (V_bary >= 0) & (W_bary >= 0)
    546 
    547549    legend_rgb = (
    548550        U_bary[..., None] * violet
    … …  
    555557    legend_rgba[..., 3] = mask_triangle.astype(float)
    556558
    557     # Loop over grid and slope
     559    # Extract data arrays
    558560    h2o = gridded_data['h2o_ice']
    559561    co2 = gridded_data['co2_ice']
    … …  
    561563    ngrid, ntime, nslope, nz = h2o.shape
    562564
     565    # Fill missing depths
     566    ti = top_index.copy().astype(int)
     567    for ig in range(ngrid):
     568        for isl in range(nslope):
     569            for t in range(1, ntime):
     570                if ti[ig, t, isl] <= 0:
     571                    ti[ig, t, isl] = ti[ig, t-1, isl]
     572
     573    # Loop over grid and slope
    563574    for ig in range(ngrid):
    564575        for isl in range(nslope):
    565576            # Compute RGB stratification over time
    566577            rgb = np.ones((nz, ntime, 3), dtype=float)
     578
     579            frac_all = np.zeros((nz, ntime, 3), dtype=float)  # store fH2O, fCO2, fDust
    567580            for t in range(ntime):
    568                 mask_z = np.arange(nz) < top_index[ig, t, isl]
    569                 if not mask_z.any():
     581                depth = ti[ig, t, isl]
     582                if depth <= 0:
    570583                    continue
    571                 cH2O = np.clip(h2o[ig, t, isl, mask_z], 0, None)
    572                 cCO2 = np.clip(co2[ig, t, isl, mask_z], 0, None)
    573                 cDust = np.clip(dust[ig, t, isl, mask_z], 0, None)
     584                cH2O = np.clip(h2o[ig, t, isl, :depth], 0, None)
     585                cCO2 = np.clip(co2[ig, t, isl, :depth], 0, None)
     586                cDust = np.clip(dust[ig, t, isl, :depth], 0, None)
    574587                total = cH2O + cCO2 + cDust
    575588                total[total == 0] = 1.0
    … …  
    577590                fCO2 = cCO2 / total
    578591                fDust = cDust / total
    579                 mix = (
    580                     np.outer(fH2O, blue)
    581                     + np.outer(fCO2, violet)
    582                     + np.outer(fDust, orange)
    583                 )
    584                 mix = np.clip(mix, 0.0, 1.0)
    585                 rgb[mask_z, t, :] = mix
    586 
    587             display_rgb = rgb[mask_elev, :, :]
     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)
     595
     596            # Mask elevation
     597            display_rgb = rgb[elevation_mask, :, :]
     598            display_frac = frac_all[elevation_mask, :, :]
     599
     600            display_rgb = rgb[elevation_mask, :, :]
     601
     602            # Compute edges for pcolormesh
     603            dt = date_time[1] - date_time[0] if len(date_time) > 1 else 1
     604            x_edges = np.concatenate([date_time, [date_time[-1] + dt]])
     605            d_e = np.diff(elev)
     606            last_e = elev[-1] + (d_e[-1] if len(d_e)>0 else 1)
     607            y_edges = np.concatenate([elev, [last_e]])
    588608
    589609            # Create figure with legend
    590610            fig, (ax_main, ax_leg) = plt.subplots(
    591                 1, 2, figsize=(12, 5), dpi=200,
     611                1, 2, figsize=(8, 4), dpi=200,
    592612                gridspec_kw={'width_ratios': [5, 1]}
    593613            )
    594614
    595615            # Main stratification panel
    596             ax_main.imshow(
     616            mesh = ax_main.pcolormesh(
     617                x_edges,
     618                y_edges,
    597619                display_rgb,
    598                 aspect='auto',
    599                 extent=[date_time[0], date_time[-1], elev.min(), elev.max()],
    600                 interpolation='nearest',
    601                 origin='lower'
     620                shading='auto',
     621                edgecolors='none'
    602622            )
    603             x_centers = np.linspace(date_time[0], date_time[-1], display_rgb.shape[1])
    604             y_centers = np.linspace(elev.min(), elev.max(), display_rgb.shape[0])
    605             attach_format_coord(ax_main, display_rgb, x_centers, y_centers, is_pcolormesh=False)
     623
     624            # Custom coordinate formatter: show time, elevation, and mixture fractions
     625            def main_format(x, y):
     626                # check bounds
     627                if x < x_edges[0] or x > x_edges[-1] or y < y_edges[0] or y > y_edges[-1]:
     628                    return ''
     629                # locate cell
     630                i = np.searchsorted(x_edges, x) - 1
     631                j = np.searchsorted(y_edges, y) - 1
     632                i = np.clip(i, 0, display_rgb.shape[1]-1)
     633                j = np.clip(j, 0, display_rgb.shape[0]-1)
     634                # get fractions
     635                fH2O, fCO2, fDust = display_frac[j, i]
     636                return f"Time={x:.2f}, Elev={y:.2f}, H2O={fH2O:.2f}, CO2={fCO2:.2f}, Dust={fDust:.2f}"
     637            ax_main.format_coord = main_format
    606638            ax_main.set_facecolor('white')
    607639            ax_main.set_title(f"Ternary mix over time (Grid point {ig+1}, Slope {isl+1})", fontweight='bold')
    … …  
    609641            ax_main.set_ylabel("Elevation (m)")
    610642
    611             # Legend panel
    612             ax_leg.imshow(
     643            # Legend panel using proper edges
     644            u_edges = np.linspace(0, 1, res+1)
     645            v_edges = np.linspace(0, np.sqrt(3)/2, res+1)
     646            ax_leg.pcolormesh(
     647                u_edges,
     648                v_edges,
    613649                legend_rgba,
    614                 extent=[0, 1, 0, np.sqrt(3)/2],
    615                 origin='lower',
    616                 interpolation='nearest'
     650                shading='auto',
     651                edgecolors='none'
    617652            )
    618             attach_format_coord(ax_leg, legend_rgba, np.array([0, 1]), np.array([0, np.sqrt(3)/2]), is_pcolormesh=False)
    619 
    620             # Draw triangle border
     653            ax_leg.set_aspect('equal')
     654
     655            # Custom coordinate formatter for legend: show barycentric fractions
     656            def legend_format(x, y):
     657                # compute barycentric coords from cartesian (x,y)
     658                V = 2 * y / np.sqrt(3)
     659                U = x - 0.5 * V
     660                W = 1 - U - V
     661                if U >= 0 and V >= 0 and W >= 0:
     662                    return f"H2O: {W:.2f}, Dust: {V:.2f}, CO2: {U:.2f}"
     663                else:
     664                    return ''
     665            ax_leg.format_coord = legend_format
     666
     667            # Draw triangle border and gridlines
    621668            triangle = np.array([[0, 0], [1, 0], [0.5, np.sqrt(3)/2], [0, 0]])
    622             ax_leg.plot(triangle[:, 0], triangle[:, 1], 'k-', linewidth=1)
    623 
    624             # Dashed gridlines
     669            ax_leg.plot(triangle[:, 0], triangle[:, 1], 'k-', linewidth=1, clip_on=False, zorder=10)
    625670            ticks = np.linspace(0.25, 0.75, 3)
    626671            for f in ticks:
    627                 ax_leg.plot([1 - f, 0.5 * (1 - f)], [0, (1 - f)*np.sqrt(3)/2], '--', color='k', linewidth=0.5)
    628                 ax_leg.plot([f, f + 0.5 * (1 - f)], [0, (1 - f)*np.sqrt(3)/2], '--', color='k', linewidth=0.5)
     672                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)
     673                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)
    629674                y = (np.sqrt(3)/2) * f
    630                 ax_leg.plot([0.5 * f, 1 - 0.5 * f], [y, y], '--', color='k', linewidth=0.5)
     675                ax_leg.plot([0.5 * f, 1 - 0.5 * f], [y, y], '--', color='k', linewidth=0.5, clip_on=False, zorder=9)
    631676
    632677            # Legend labels
    … …  
    636681            ax_leg.axis('off')
    637682
     683            # Save figure
    638684            plt.tight_layout()
    639 
    640             # Save figure
    641685            fname = os.path.join(output_folder, f"layering_rgb_evolution_ig{ig+1}_is{isl+1}.png")
    642686            fig.savefig(fname, dpi=150, bbox_inches='tight')
    … …  
    678722    vmin, vmax = -2, 1
    679723    epsilon = 1e-6
    680 
     724   
    681725    # Loop over grids and slopes
    682726    for ig in range(ngrid):
    683727        for isl in range(nslope):
     728            ti = top_index[ig, :, isl].copy().astype(int)
     729            for t in range(1, ntime):
     730                if ti[t] <= 0:
     731                    ti[t] = ti[t-1]
     732
     733            # Compute log10(dust/ice) profile at each time step
    684734            log_ratio_array = np.full((nz, ntime), np.nan, dtype=np.float32)
    685 
    686             # Compute log10(dust/ice) profile at each time step
    687735            for t in range(ntime):
    688                 zmax = top_index[ig, t, isl]
     736                zmax = ti[t]
    689737                if zmax <= 0:
    690738                    continue
    … …  
    710758            # Plot
    711759            fig, ax = plt.subplots(figsize=(8, 6), dpi=150)
     760            im = ax.pcolormesh(
     761                date_time,
     762                elev,
     763                ratio_display,
     764                shading='auto',
     765                cmap='managua_r',
     766                norm=LogNorm(vmin=10**vmin, vmax=10**vmax),
     767            )
     768            x_edges = np.concatenate([date_time, [date_time[-1] + (date_time[-1]-date_time[-2])]])
     769            attach_format_coord(ax, ratio_display, x_edges, np.concatenate([elev, [elev[-1] + (elev[-1]-elev[-2])]]), is_pcolormesh=True)
    712770            ax.set_title(f"Dust-to-Ice ratio over time (Grid point {ig+1}, Slope {isl+1})", fontweight='bold')
    713             im = ax.imshow(
    714                 ratio_display,
    715                 aspect='auto',
    716                 extent=[date_time[0], date_time[-1], elev.min(), elev.max()],
    717                 origin='lower',
    718                 interpolation='nearest',
    719                 cmap='managua_r',
    720                 norm=LogNorm(vmin=10**vmin, vmax=10**vmax)
    721             )
    722             x_centers = np.linspace(date_time[0], date_time[-1], ratio_display.shape[1])
    723             y_centers = np.linspace(elev.min(), elev.max(), ratio_display.shape[0])
    724             attach_format_coord(ax, ratio_display, x_centers, y_centers, is_pcolormesh=False)
     771            ax.set_xlabel('Time (Mars years)')
     772            ax.set_ylabel('Elevation (m)')
    725773
    726774            # Add colorbar with simplified ratio labels
    … …  
    827875    # Plot
    828876    fig, axes = plt.subplots(3, 1, figsize=(8, 10), sharex=True)
    829     fig.suptitle("Orbital Parameters vs Simulated Time", fontsize=14, fontweight='bold')
    830 
    831     axes[0].plot(date_time, obl_interp, 'r+', linestyle='-')
     877    fig.suptitle("Orbital parameters vs simulated time", fontsize=14, fontweight='bold')
     878
     879    axes[0].plot(date_time, obl_interp, 'r-', marker='+')
    832880    axes[0].set_ylabel("Obliquity (°)")
    833881    axes[0].grid(True)
    834882
    835     axes[1].plot(date_time, ecc_interp, 'b+', linestyle='-')
     883    axes[1].plot(date_time, ecc_interp, 'b-', marker='+')
    836884    axes[1].set_ylabel("Eccentricity")
    837885    axes[1].grid(True)
    838886
    839     axes[2].plot(date_time, lsp_interp, 'g+', linestyle='-')
    840     axes[2].set_ylabel("Ls p (°)")
     887    axes[2].plot(date_time, lsp_interp, 'g-', marker='+')
     888    axes[2].set_ylabel("Ls of perihelion (°)")
    841889    axes[2].set_xlabel("Time (Mars years)")
    842890    axes[2].grid(True)
    843891
    844892    plt.tight_layout(rect=[0, 0, 1, 0.96])
    845     fname = os.path.join(output_folder, "orbital_parameters.png")
     893    fname = os.path.join(output_folder, "orbital_parameters_laskar.png")
    846894    fig.savefig(fname, dpi=150)
     895
     896
     897def mars_ls(pday, peri_day, e_elips, year_day, lsperi=0.0):
     898    """
     899    Compute solar longitude (Ls) in radians for a given Mars date array 'pday'.
     900    Returns Ls in degrees [0, 360).
     901    """
     902    zz = (pday - peri_day) / year_day
     903    zanom = 2 * np.pi * (zz - np.round(zz))
     904    xref = np.abs(zanom)
     905
     906    # Solve Kepler's equation via Newton–Raphson
     907    zx0 = xref + e_elips * np.sin(xref)
     908    for _ in range(10):
     909        f  = zx0 - e_elips * np.sin(zx0) - xref
     910        fp = 1 - e_elips * np.cos(zx0)
     911        dz = -f / fp
     912        zx0 += dz
     913        if np.all(np.abs(dz) <= 1e-7):
     914            break
     915
     916    zx0 = np.where(zanom < 0, -zx0, zx0)
     917    zteta = 2 * np.arctan(
     918        np.sqrt((1 + e_elips) / (1 - e_elips)) * np.tan(zx0 / 2)
     919    )
     920    psollong = np.mod(zteta + lsperi, 2 * np.pi)
     921
     922    return np.degrees(psollong)
     923
     924
     925def read_orbital_data_nc(starts_folder, infofile=None):
     926    """
     927    Read orbital parameters from restartfi_postPEM*.nc files in starts_folder.
     928    """
     929    if not os.path.isdir(starts_folder):
     930        raise ValueError(f"Invalid starts_folder '{starts_folder}': not a directory.")
     931
     932    # Read simulation time mapping if provided
     933    if infofile:
     934        dates_yr, martian_to_earth = read_infofile(infofile)
     935    else:
     936        dates_yr = None
     937
     938    pattern = os.path.join(starts_folder, "restartfi_postPEM*.nc")
     939    files = glob(pattern)
     940    if not files:
     941        raise FileNotFoundError(f"No NetCDF restart files found matching {pattern}")
     942
     943    def extract_number(path):
     944        name = os.path.basename(path)
     945        prefix = 'restartfi_postPEM'
     946        if name.startswith(prefix) and name.endswith('.nc'):
     947            num_str = name[len(prefix):-3]
     948            if num_str.isdigit():
     949                return int(num_str)
     950        return float('inf')
     951
     952    files = sorted(files, key=extract_number)
     953
     954    all_year_day, all_peri, all_aphe, all_date_peri, all_obl = [], [], [], [], []
     955    for nc_path in files:
     956        with Dataset(nc_path, 'r') as nc:
     957            ctrl = nc.variables['controle'][:]
     958            all_year_day.append(ctrl[13])
     959            all_peri.append(ctrl[14])
     960            all_aphe.append(ctrl[15])
     961            all_date_peri.append(ctrl[16])
     962            all_obl.append(ctrl[17])
     963
     964    year_day      = np.array(all_year_day)
     965    perihelion    = np.array(all_peri)
     966    aphelion      = np.array(all_aphe)
     967    date_peri_day = np.array(all_date_peri)
     968    obliquity     = np.array(all_obl)
     969
     970    eccentricity  = (aphelion - perihelion) / (aphelion + perihelion)
     971    ls_perihelion = mars_ls(
     972        date_peri_day,
     973        date_peri_day,
     974        eccentricity,
     975        year_day
     976    )
     977
     978    return dates_yr, obliquity, eccentricity, ls_perihelion
     979
     980
     981def plot_orbital_parameters_nc(starts_folder, infofile, date_time, output_folder="."):
     982    """
     983    Plot the evolution of obliquity, eccentricity and Ls p coming from simulation data
     984    versus simulated time.
     985    """
     986    # Read orbital data
     987    times_yr, obl, ecc, lsp = read_orbital_data_nc(starts_folder, infofile)
     988
     989    fargs = dict(kind='linear', bounds_error=False, fill_value='extrapolate')
     990    obl_i = interp1d(times_yr, obl, **fargs)(date_time)
     991    ecc_i = interp1d(times_yr, ecc, **fargs)(date_time)
     992    lsp_i = interp1d(times_yr, lsp, **fargs)(date_time)
     993
     994    fig, axes = plt.subplots(3,1, figsize=(8,10), sharex=True)
     995    fig.suptitle("Orbital parameters vs simulated time", fontsize=14, fontweight='bold')
     996
     997    # Plot
     998    axes[0].plot(date_time, obl_i, 'r-', marker='+')
     999    axes[0].set_ylabel("Obliquity (°)")
     1000    axes[0].grid(True)
     1001
     1002    axes[1].plot(date_time, ecc_i, 'b-', marker='+')
     1003    axes[1].set_ylabel("Eccentricity")
     1004    axes[1].grid(True)
     1005
     1006    axes[2].plot(date_time, lsp_i, 'g-', marker='+')
     1007    axes[2].set_ylabel("Ls of perihelion (°)")
     1008    axes[2].set_xlabel("Time (Mars years)")
     1009    axes[2].grid(True)
     1010
     1011    plt.tight_layout(rect=[0,0,1,0.96])
     1012    outname = os.path.join(output_folder, "orbital_parameters_simu.png")
     1013    fig.savefig(outname, dpi=150)
     1014
     1015
     1016def plot_dust_to_ice_ratio_with_obliquity(
     1017    starts_folder,
     1018    infofile,
     1019    gridded_data,
     1020    ref_grid,
     1021    top_index,
     1022    heights_data,
     1023    date_time,
     1024    exclude_sub=False,
     1025    output_folder="."
     1026):
     1027    """
     1028    Plot the dust-to-ice ratio over time as a heatmap, and overlay the evolution of
     1029    obliquity on a secondary y-axis.
     1030    """
     1031    h2o = gridded_data['h2o_ice']
     1032    dust = gridded_data['dust']
     1033    ngrid, ntime, nslope, nz = h2o.shape
     1034
     1035    # Read orbital data
     1036    times_yr, obl, _, _ = read_orbital_data_nc(starts_folder, infofile)
     1037    fargs = dict(kind='linear', bounds_error=False, fill_value='extrapolate')
     1038    obliquity = interp1d(times_yr, obl, **fargs)(date_time)
     1039
     1040    # Computed total height
     1041    for ig in range(ngrid):
     1042        for isl in range(nslope):
     1043            total_height_t = np.zeros(ntime, dtype=float)
     1044
     1045            for t_idx in range(ntime):
     1046                h_mat = heights_data[t_idx][isl]
     1047                raw_h = h_mat[ig, :]
     1048                valid_mask = (~np.isnan(raw_h)) & (raw_h != 0.0)
     1049                if np.any(valid_mask):
     1050                    h_valid = raw_h[valid_mask]
     1051                    total_height_t[t_idx] = np.max(h_valid)
     1052
     1053    # Compute the per-interval sign of height change
     1054    dh = np.diff(total_height_t)
     1055    signs = np.sign(dh)
     1056    color_map = { 1: 'green', -1: 'red', 0: 'orange' }
     1057
     1058    # Elevation mask
     1059    if exclude_sub:
     1060        elevation_mask = (ref_grid >= 0.0)
     1061        elev = ref_grid[elevation_mask]
     1062    else:
     1063        elevation_mask = np.ones_like(ref_grid, dtype=bool)
     1064        elev = ref_grid
     1065
     1066    # Custom colormap: blue (ice) to orange (dust)
     1067    blue = np.array([0, 0, 255]) / 255
     1068    orange = np.array([255, 165, 0]) / 255
     1069    custom_cmap = LinearSegmentedColormap.from_list('BlueOrange', [blue, orange], N=256)
     1070
     1071    # Log‑ratio bounds and small epsilon to avoid log(0)
     1072    vmin, vmax = -2, 1
     1073    epsilon = 1e-6
     1074
     1075    # Loop over grids and slopes
     1076    for ig in range(ngrid):
     1077        for isl in range(nslope):
     1078            ti = top_index[ig, :, isl].copy().astype(int)
     1079            for t in range(1, ntime):
     1080                if ti[t] <= 0:
     1081                    ti[t] = ti[t-1]
     1082
     1083            # Compute log10(dust/ice) profile at each time step
     1084            log_ratio_array = np.full((nz, ntime), np.nan, dtype=np.float32)
     1085            for t in range(ntime):
     1086                zmax = ti[t]
     1087                if zmax <= 0:
     1088                    continue
     1089
     1090                h2o_profile = np.clip(h2o[ig, t, isl, :zmax], 0, None)
     1091                dust_profile = np.clip(dust[ig, t, isl, :zmax], 0, None)
     1092
     1093                with np.errstate(divide='ignore', invalid='ignore'):
     1094                    ratio_profile = np.where(
     1095                        h2o_profile > 0,
     1096                        dust_profile / h2o_profile,
     1097                        10**(vmax + 1)
     1098                    )
     1099                    log_ratio = np.log10(ratio_profile + epsilon)
     1100                    log_ratio = np.clip(log_ratio, vmin, vmax)
     1101
     1102                log_ratio_array[:zmax, t] = log_ratio
     1103
     1104            # Convert back to linear ratio and mask
     1105            ratio_array = 10**log_ratio_array
     1106            ratio_display = ratio_array[elevation_mask, :]
     1107
     1108            # Plot
     1109            fig, ax = plt.subplots(figsize=(8, 6), dpi=150)
     1110            im = ax.pcolormesh(
     1111                date_time,
     1112                elev,
     1113                ratio_display,
     1114                shading='auto',
     1115                cmap='managua_r',
     1116                norm=LogNorm(vmin=10**vmin, vmax=10**vmax),
     1117            )
     1118            x_edges = np.concatenate([date_time, [date_time[-1] + (date_time[-1]-date_time[-2])]])
     1119            attach_format_coord(ax, ratio_display, x_edges, np.concatenate([elev, [elev[-1] + (elev[-1]-elev[-2])]]), is_pcolormesh=True)
     1120            ax.set_title(f"Dust-to-Ice ratio over time (Grid point {ig+1}, Slope {isl+1})", fontweight='bold')
     1121            ax.set_xlabel('Time (Mars years)')
     1122            ax.set_ylabel('Elevation (m)')
     1123
     1124            # Add colorbar
     1125            cbar = fig.colorbar(im, ax=ax, orientation='vertical', pad=0.15)
     1126            cbar.set_label('Dust / H₂O ice (ratio)')
     1127            cbar.set_ticks([1e-2, 1e-1, 1, 1e1])
     1128            cbar.set_ticklabels(['1:100', '1:10', '1:1', '10:1'])
     1129
     1130            # Overlay obliquity on secondary y-axis
     1131            ax2 = ax.twinx()
     1132            for i in range(len(dh)):
     1133                ax2.plot(
     1134                    [date_time[i], date_time[i+1]],
     1135                    [obliquity[i], obliquity[i+1]],
     1136                    color=color_map[signs[i]],
     1137                    marker='+',
     1138                    linewidth=1.5
     1139                )
     1140            ax2.set_ylabel('Obliquity (°)')
     1141            ax2.tick_params(axis='y')
     1142            ax2.grid(False)
     1143
     1144            # Save
     1145            os.makedirs(output_folder, exist_ok=True)
     1146            outname = os.path.join(
     1147                output_folder,
     1148                f'dust_ice_obliquity_grid{ig+1}_slope{isl+1}.png'
     1149            )
     1150            plt.tight_layout()
     1151            fig.savefig(outname, dpi=150)
    8471152
    8481153
    … …  
    9131218        exclude_sub=exclude_sub, output_folder="."
    9141219    )
    915     plot_dust_to_ice_ratio_over_time(
     1220    #plot_dust_to_ice_ratio_over_time(
     1221    #    gridded_data, ref_grid, top_index, heights_data, date_time,
     1222    #    exclude_sub=exclude_sub, output_folder="."
     1223    #)
     1224    plot_dust_to_ice_ratio_with_obliquity(
     1225        folder_path, infofile,
    9161226        gridded_data, ref_grid, top_index, heights_data, date_time,
    9171227        exclude_sub=exclude_sub, output_folder="."
    9181228    )
    919     plot_strata_count_and_total_height(heights_data, date_time, output_folder=".")
     1229    #plot_strata_count_and_total_height(heights_data, date_time, output_folder=".")
    9201230
    9211231    # 14) Plot orbital parameters
    922     plot_orbital_parameters(infofile, orbfile, date_time, output_folder=".")
     1232    #plot_orbital_parameters(infofile, orbfile, date_time, output_folder=".")
     1233    plot_orbital_parameters_nc(folder_path, infofile, date_time, output_folder=".")
    9231234
    9241235    # 15) Show all figures

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

    r3842 r3850  
    6969use writediagpem_mod,           only: writediagpem, writediagsoilpem
    7070use co2condens_mod,             only: CO2cond_ps
    71 use layering_mod,               only: layering, del_layering, make_layering, layering_algo, subsurface_ice_layering,            &
    72                                       ptrarray, stratum, get_nb_str_max, nb_str_max, is_dust_lag, is_co2ice_str, is_h2oice_str, print_layering
     71use layering_mod,               only: layering, del_layering, make_layering, layering_algo, subsurface_ice_layering, &
     72                                      ptrarray, stratum, get_nb_str_max, nb_str_max, is_dust_lag, is_co2ice_str, is_h2oice_str
    7373use dyn_ss_ice_m_mod,           only: dyn_ss_ice_m
    7474use parse_args_mod,             only: parse_args

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

    r3842 r3850  
    393393                do islope = 1,nslope
    394394                    call ini_layering(layerings_map(ig,islope))
    395                     print*, 'coucou', 1.05*ini_huge_h2oice/rho_h2oice, layerings_map(ig,islope)%top%top_elevation
    396395                    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.)
    397396                enddo

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

    r3742 r3850  
    3131! Input Variables
    3232!--------------------------------------------------------
    33 real, intent(in) :: i_myear   ! Number of simulated Martian years
     33real, intent(in) :: i_myear     ! Number of simulated Martian years
    3434real, intent(in) :: i_myear_leg ! Number of iterations of the PEM
    3535
    … …  
    6262Lsp = lsperi*180./pi         ! We convert in degree
    6363
     64print*, year_bp_ini,i_myear,i_myear_leg
     65
    6466write(*,*) 'recomp_orb_param, Old year =', Year - i_myear_leg
    6567write(*,*) 'recomp_orb_param, Old obl. =', obliquit
    … …  
    6971
    7072found_year = .false.
    71 do ilask = last_ilask,2,-1
    72     xa = yearlask(ilask)
    73     xb = yearlask(ilask - 1)
    74     if (xa <= Year .and. Year < xb) then
    75         ! Obliquity
    76         if (var_obl) then
    77             ya = obllask(ilask)
    78             yb = obllask(ilask - 1)
    79             obliquit = (Year - xa)*(yb - ya)/(xb - xa) + ya
     73if (Year < 0.) then
     74    do ilask = last_ilask,2,-1
     75        xa = yearlask(ilask)
     76        xb = yearlask(ilask - 1)
     77        if (xa <= Year .and. Year < xb) then
     78            ! Obliquity
     79            if (var_obl) then
     80                ya = obllask(ilask)
     81                yb = obllask(ilask - 1)
     82                obliquit = (Year - xa)*(yb - ya)/(xb - xa) + ya
     83            endif
     84
     85            ! Eccentricity
     86            if (var_ecc) then
     87                ya = ecclask(ilask)
     88                yb = ecclask(ilask - 1)
     89                e_elips = (Year - xa)*(yb - ya)/(xb - xa) + ya
     90            endif
     91
     92            ! Lsp
     93            if (var_lsp) then
     94                ya = lsplask(ilask)
     95                yb = lsplask(ilask - 1)
     96                if (abs(yb - ya) > 300.) then ! If modulo is "crossed" through the interval
     97                    if (ya < yb) then ! Lsp should be decreasing
     98                        ya = ya + 360.
     99                    else ! Lsp should be increasing
     100                        yb = yb + 360.
     101                    endif
     102                endif
     103                Lsp = modulo((Year - xa)*(yb - ya)/(xb - xa) + ya,360.)
     104            endif
     105            found_year = .true.
     106            exit ! The loop is left as soon as the right interval is found
    80107        endif
    81 
    82         ! Eccentricity
    83         if (var_ecc) then
    84             ya = ecclask(ilask)
    85             yb = ecclask(ilask - 1)
    86             e_elips = (Year - xa)*(yb - ya)/(xb - xa) + ya
    87         endif
    88 
    89         ! Lsp
    90         if (var_lsp) then
    91             ya = lsplask(ilask)
    92             yb = lsplask(ilask - 1)
    93             if (abs(yb - ya) > 300.) then ! If modulo is "crossed" through the interval
    94                 if (ya < yb) then ! Lsp should be decreasing
    95                     ya = ya + 360.
    96                 else ! Lsp should be increasing
    97                     yb = yb + 360.
    98                 endif
    99             endif
    100             Lsp = modulo((Year - xa)*(yb - ya)/(xb - xa) + ya,360.)
    101         endif
    102         found_year = .true.
    103         exit ! The loop is left as soon as the right interval is found
    104     endif
    105 enddo
     108    enddo
     109else if (Year >= 0 .and. Year < 100.) then ! For present orbital characteristics
     110    if (var_obl) obliquit = 25.19
     111    if (var_ecc) e_elips = 0.0934
     112    if (var_lsp) Lsp = 251.
     113    found_year = .true.
     114endif
    106115
    107116if (.not. found_year) error stop 'The new year could not be found in the file "obl_ecc_lsp.asc".'
     117
     118write(*,*) 'recomp_orb_param, New obl. =', obliquit
     119write(*,*) 'recomp_orb_param, New ecc. =', e_elips
     120write(*,*) 'recomp_orb_param, New Lsp  =', Lsp
    108121
    109122halfaxe = 227.94
    … …  
    118131#endif
    119132
    120 write(*,*) 'recomp_orb_param, New obl. =', obliquit
    121 write(*,*) 'recomp_orb_param, New ecc. =', e_elips
    122 write(*,*) 'recomp_orb_param, New Lsp  =', Lsp
    123 
    124133END SUBROUTINE recomp_orb_param
    125134
Note: See TracChangeset for help on using the changeset viewer.