Ignore:
Timestamp:
Jul 7, 2025, 6:09:01 PM (27 hours ago)
Author:
jbclement
Message:

PEM:

  • Correction of a bug in the launching script.
  • Update of "visu_evol_layering.py", in particular to show value at cursor for 2D heatmaps.
  • Few cleanings.

JBC

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

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.COMMON/libf/evolution/changelog.txt

    r3821 r3840  
    722722
    723723== 01/07/2025 == JBC
    724 Addition of the flag 'impose_dust_ratio' (default: .false.) to impose a dust-to-ice ratio instead of a dust tendency for the layering algorithm. The ratio value can be set in the "run_PEM.def" by 'dust_to_ice_ratio' (default: 0.01).
     724Addition of the flag 'impose_dust_ratio' (default: .false.) to impose a dust-to-ice ratio instead of a dust tendency for the layering algorithm. The ratio value can be set in the "run_PEM.def" by 'dust2ice_ratio' (default: 0.01).
     725
     726== 07/07/2025 == JBC
     727- Correction of a bug in the launching script.
     728- Update of "visu_evol_layering.py", in particular to show value at cursor for 2D heatmaps.
     729- Few cleanings.
  • trunk/LMDZ.COMMON/libf/evolution/deftank/lib_launchPEM.sh

    r3820 r3840  
    174174        echo "Run \"PCM $iPCM\" ($ii/$3)..."
    175175        if [ $1 -eq 0 ]; then # Mode: processing scripts
    176             sed -i "s/^k=[0-9]\+$/k=$(echo "$ii - $3 + 2" | bc)/" PCMrun.job
     176            sed -i "s/^k=-\?[0-9]\+$/k=$(echo "$ii - $3 + 2" | bc)/" PCMrun.job
    177177            ./PCMrun.job
    178178            if [ $? -ne 0 ]; then
     
    182182            cp PCMrun.job PCMrun${iPCM}.job
    183183            sed -i -E "/^$name_job/s/(.*[^0-9])([0-9]+)(_[^0-9]*)?$/\1${iPCM}\3/" PCMrun${iPCM}.job
    184             sed -i "s/^k=[0-9]\+$/k=$(echo "$ii - $3 + 2" | bc)/" PCMrun${iPCM}.job
     184            sed -i "s/^k=-\?[0-9]\+$/k=$(echo "$ii - $3 + 2" | bc)/" PCMrun${iPCM}.job
    185185            jobID=$(eval "$submit_job PCMrun${iPCM}.job")
    186186            # Create a file to cancel the dependent jobs of the cycle
     
    201201            echo "Run \"PCM $iPCM\" ($i/$3)..."
    202202            if [ $1 -eq 0 ]; then # Mode: processing scripts
    203                 sed -i "s/^k=[0-9]\+$/k=$(echo "$i - $3 + 2" | bc)/" PCMrun.job
     203                sed -i "s/^k=-\?[0-9]\+$/k=$(echo "$i - $3 + 2" | bc)/" PCMrun.job
    204204                ./PCMrun.job
    205205                if [ $? -ne 0 ]; then
     
    209209                cp PCMrun.job PCMrun${iPCM}.job
    210210                sed -i -E "/^$name_job/s/(.*[^0-9])([0-9]+)(_[^0-9]*)?$/\1${iPCM}\3/" PCMrun${iPCM}.job
    211                 sed -i "s/^k=[0-9]\+$/k=$(echo "$i - $3 + 2" | bc)/" PCMrun${iPCM}.job
     211                sed -i "s/^k=-\?[0-9]\+$/k=$(echo "$i - $3 + 2" | bc)/" PCMrun${iPCM}.job
    212212                jobID=$(eval "$submit_dependjob=afterok:${jobID} PCMrun${iPCM}.job")
    213213                echo $kill_job $jobID >> kill_launchPEM.sh
  • trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def

    r3821 r3840  
    1818
    1919# Do you want to vary the obliquity when following "obl_ecc_lsp.asc"? Default = .true.
    20 # var_obl=.false.
     20# var_obl=.true.
    2121
    2222# Do you want to vary the eccentricity when following "obl_ecc_lsp.asc"? Default = .true.
    23 # var_ecc=.false.
     23# var_ecc=.true.
    2424
    2525# Do you want to vary the ls perihelie when following "obl_ecc_lsp.asc"? Default = .true.
    26 # var_lsp=.false.
     26# var_lsp=.true.
    2727
    2828# Time step length of the PEM in Martian years? Default = 1.
     
    3131#---------- Stopping criteria parameters ----------#
    3232# If evol_orbit_pem=.false., maximal number of iterations if no stopping criterion is reached. Default=100000000
    33 # Max_iter_pem=10
     33# Max_iter_pem=100000000
    3434
    3535# Acceptance rate of sublimating H2O ice surface change. Default = 0.2
     
    5656
    5757# Do you want to run with adsoprtion in the PEM? Default = .false.
    58 # adsorption_pem=.true.
     58# adsorption_pem=.false.
    5959
    6060# Do you want to modify the soil thermal properties with the pressure? Default = .false.
     
    7171
    7272# Is the ice table computed at equilibrium? Default = .true.
    73 # icetable_equilibrium=.false.
     73# icetable_equilibrium=.true.
    7474
    7575# Is the ice table computed with the dynamic method? Default = .false.
     
    9999
    100100# Do you want the CO2 ice to flow along subslope inside a cell? Default = .true.
    101 # co2ice_flow=.false.
     101# co2ice_flow=.true.
    102102
    103103#---------- Layering parameters ----------#
    104104# Do you want to run with the layering algorithm? Default = .false.
    105 # layering=.true.
     105# layering=.false.
    106106 
    107107# Value of the dust tendency. Default = 1.e-3 kg.m-2.y-1
     
    109109 
    110110# Do you want to impose a dust-to-ice ratio instead of a dust tendency? Default = .false.
    111 # impose_dust_ratio=.true.
     111# impose_dust_ratio=.false.
    112112 
    113113# If impose_dust_ratio=.true., value of the dust-to-ice ratio. Default = 0.01
    114 # dust2ice_ratio=0.05
     114# dust2ice_ratio=0.01
    115115
    116116# Some definitions for the physics, in file 'callphys.def'
  • trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py

    r3819 r3840  
    379379
    380380
     381def attach_format_coord(ax, mat, x, y, is_pcolormesh=True):
     382    """
     383    Attach a format_coord function to the axes to display x, y, and value at cursor.
     384    Works for both pcolormesh and imshow style grids.
     385    """
     386    # Determine dimensions
     387    if mat.ndim == 2:
     388        ny, nx = mat.shape
     389    elif mat.ndim == 3 and mat.shape[2] in (3, 4):
     390        ny, nx, nc = mat.shape
     391    else:
     392        raise ValueError(f"Unsupported mat shape {mat.shape}")
     393    # Edges or extents
     394    if is_pcolormesh:
     395        xedges, yedges = x, y
     396    else:
     397        x0, x1 = x.min(), x.max()
     398        y0, y1 = y.min(), y.max()
     399
     400    def format_coord(xp, yp):
     401        # Map to indices
     402        if is_pcolormesh:
     403            col = np.searchsorted(xedges, xp) - 1
     404            row = np.searchsorted(yedges, yp) - 1
     405        else:
     406            col = int((xp - x0) / (x1 - x0) * nx)
     407            row = int((yp - y0) / (y1 - y0) * ny)
     408        # Within bounds?
     409        if 0 <= row < ny and 0 <= col < nx:
     410            if mat.ndim == 2:
     411                v = mat[row, col]
     412                return f"x={xp:.3g}, y={yp:.3g}, val={v:.3g}"
     413            else:
     414                vals = mat[row, col]
     415                txt = ", ".join(f"{vv:.3g}" for vv in vals[:3])
     416                return f"x={xp:.3g}, y={yp:.3g}, val=({txt})"
     417        return f"x={xp:.3g}, y={yp:.3g}"
     418
     419    ax.format_coord = format_coord
     420
     421
    381422def plot_stratification_over_time(
    382423    gridded_data,
     
    454495                    vmax=vmax
    455496                )
     497                x_edges = np.concatenate([date_time, [date_time[-1] + (date_time[-1]-date_time[-2])]])
     498                attach_format_coord(ax, mat, x_edges, np.concatenate([sub_ref_grid, [sub_ref_grid[-1] + (sub_ref_grid[-1]-sub_ref_grid[-2])]]), is_pcolormesh=True)
    456499                ax.set_title(titles[idx], fontsize=12)
    457500                ax.set_xlabel("Time (Mars years)")
     
    558601                origin='lower'
    559602            )
     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)
    560606            ax_main.set_facecolor('white')
    561607            ax_main.set_title(f"Ternary mix over time (Grid point {ig+1}, Slope {isl+1})", fontweight='bold')
     
    570616                interpolation='nearest'
    571617            )
     618            attach_format_coord(ax_leg, legend_rgba, np.array([0, 1]), np.array([0, np.sqrt(3)/2]), is_pcolormesh=False)
    572619
    573620            # Draw triangle border
     
    663710            # Plot
    664711            fig, ax = plt.subplots(figsize=(8, 6), dpi=150)
     712            ax.set_title(f"Dust-to-Ice ratio over time (Grid point {ig+1}, Slope {isl+1})", fontweight='bold')
    665713            im = ax.imshow(
    666714                ratio_display,
     
    672720                norm=LogNorm(vmin=10**vmin, vmax=10**vmax)
    673721            )
     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)
    674725
    675726            # Add colorbar with simplified ratio labels
  • trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90

    r3821 r3840  
    723723dh_co2ice = d_co2ice/rho_co2ice
    724724dh_h2oice = d_h2oice/rho_h2oice
    725 ! To disable dust sedimenation when there is ice sublimation (because dust tendency is fixed)
    726 if (dh_co2ice >= 0. .and. dh_h2oice >= 0.) then
     725dh_dust = 0.
     726if (dh_h2oice >= 0. .and. dh_co2ice >= 0.) then ! To put a dust sedimentation tendency only when ice is condensing
    727727    if (impose_dust_ratio) then
    728728        dh_dust = dust2ice_ratio*(dh_h2oice + dh_co2ice)
     
    730730        dh_dust = d_dust/rho_dust
    731731    endif
    732 else
    733     dh_dust = 0.
    734732endif
    735733zshift_surf = this%top%top_elevation
  • trunk/LMDZ.COMMON/libf/evolution/pem.F90

    r3836 r3840  
    7070use co2condens_mod,             only: CO2cond_ps
    7171use 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, &
    73                                       print_layering
     72                                      ptrarray, stratum, get_nb_str_max, nb_str_max, is_dust_lag, is_co2ice_str, is_h2oice_str
    7473use dyn_ss_ice_m_mod,           only: dyn_ss_ice_m
    7574use parse_args_mod,             only: parse_args
Note: See TracChangeset for help on using the changeset viewer.