Changeset 3851


Ignore:
Timestamp:
Jul 16, 2025, 3:25:48 PM (4 days ago)
Author:
jbclement
Message:

PEM:

  • Making the computation of CO2 mass balance more robust, especially regarding 'CO2cond_ps'.
  • Small correction about the dust tendency management for the layering algorithm.
  • Small improvement of the visualization in "visu_evol_layering.py".
  • File "log_launchPEM.txt" renamed into "launchPEM.log".

JBC

Location:
trunk
Files:
10 edited

Legend:

Unmodified
Added
Removed

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

    r3850 r3851  
    740740- Correction of numbering in "lib_launchPEM.sh" when requesting a relaunch.
    741741- Update of "PEMrun.job" due to new parsing of arguments introduced in r3836/r3837.
     742
     743== 16/07/2025 == JBC
     744- Making the computation of CO2 mass balance more robust, especially regarding 'CO2cond_ps'.
     745- Small correction about the dust tendency management for the layering algorithm.
     746- Small improvement of the visualization in "visu_evol_layering.py".
     747- File "log_launchPEM.txt" renamed into "launchPEM.log".

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

    r3667 r3851  
    44### GENOA nodes accommodate 96 cores
    55#SBATCH --constraint=GENOA
    6 ### Number of Nodes to use
     6### Number of nodes/cores to use
    77#SBATCH --nodes=1
    88#SBATCH --ntasks-per-node=24

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

    r3850 r3851  
    44### GENOA nodes accommodate 96 cores
    55#SBATCH --constraint=GENOA
    6 ### Number of Nodes to use
     6### Number of nodes/cores to use
    77#SBATCH --nodes=1
    88#SBATCH --ntasks=1

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

    r3579 r3851  
    1717rm -f *.out
    1818rm -f kill_launchPEM.sh
    19 rm -f log_launchPEM.txt
     19rm -f launchPEM.log
    2020rm -f out_run*
    2121rm -f run.def

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

    r3639 r3851  
    4848    # Starting from scratch
    4949    echo "The launching script is starting!"
    50     echo "The output file is \"log_launchPEM.txt\"."
    51     exec > log_launchPEM.txt 2>&1
     50    echo "The output file is \"launchPEM.log\"."
     51    exec > launchPEM.log 2>&1
    5252    echo "Beginning of the launching script for the PEM simulation."
    5353    date
    … …  
    5959    # Starting a new cycle
    6060    if [ $1 = "new" ]; then
    61         exec >> log_launchPEM.txt 2>&1
     61        exec >> launchPEM.log 2>&1
    6262        echo
    6363        echo "This is a new cycle for the PEM simulation."
    … …  
    111111            fi
    112112        done
    113         exec >> log_launchPEM.txt 2>&1
     113        exec >> launchPEM.log 2>&1
    114114        echo
    115115        echo "This is a relaunch for the PEM simulation from the run \"$relaunch $irelaunch\"."
    … …  
    138138    # Continuing the PEM run
    139139    elif [ $1 = "cont" ]; then
    140         exec >> log_launchPEM.txt 2>&1
     140        exec >> launchPEM.log 2>&1
    141141        echo
    142142        echo "This is a continuation of the previous PEM run."

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

    r3850 r3851  
    10301030    """
    10311031    h2o = gridded_data['h2o_ice']
     1032    co2 = gridded_data['co2_ice']
    10321033    dust = gridded_data['dust']
    10331034    ngrid, ntime, nslope, nz = h2o.shape
    … …  
    10771078        for isl in range(nslope):
    10781079            ti = top_index[ig, :, isl].copy().astype(int)
     1080            frac_all = np.zeros((nz, ntime, 3), dtype=float)  # store fH2O, fCO2, fDust
    10791081            for t in range(1, ntime):
    10801082                if ti[t] <= 0:
    … …  
    10871089                if zmax <= 0:
    10881090                    continue
     1091                cH2O = np.clip(h2o[ig, t, isl, :zmax], 0, None)
     1092                cCO2 = np.clip(co2[ig, t, isl, :zmax], 0, None)
     1093                cDust = np.clip(dust[ig, t, isl, :zmax], 0, None)
     1094                total = cH2O + cCO2 + cDust
     1095                total[total == 0] = 1.0
     1096                fH2O = cH2O / total
     1097                fCO2 = cCO2 / total
     1098                fDust = cDust / total
     1099                frac_all[:zmax, t, :] = np.stack([fH2O, fCO2, fDust], axis=1)
    10891100
    10901101                h2o_profile = np.clip(h2o[ig, t, isl, :zmax], 0, None)
    … …  
    11051116            ratio_array = 10**log_ratio_array
    11061117            ratio_display = ratio_array[elevation_mask, :]
     1118            display_frac = frac_all[elevation_mask, :, :]
    11071119
    11081120            # Plot
    … …  
    11171129            )
    11181130            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)
     1131            y_edges = np.concatenate([elev, [elev[-1] + (elev[-1]-elev[-2])]])
     1132            def format_coord_all(x, y):
     1133                # check bounds
     1134                if x < x_edges[0] or x > x_edges[-1] or y < y_edges[0] or y > y_edges[-1]:
     1135                    return ''
     1136                # locate cell
     1137                i = np.searchsorted(x_edges, x) - 1
     1138                j = np.searchsorted(y_edges, y) - 1
     1139                i = np.clip(i, 0, display_frac.shape[1]-1)
     1140                j = np.clip(j, 0, display_frac.shape[0]-1)
     1141                # get fractions
     1142                fH2O  = display_frac[j, i, 0]
     1143                fDust = display_frac[j, i, 2]
     1144                obl   = np.interp(x, date_time, obliquity)
     1145                return f"Time={x:.2f}, Elev={y:.2f}, H2O={fH2O:.2f}, Dust={fDust:.2f}, Obl={obl:.2f}°"
     1146
     1147            ax.format_coord = format_coord_all
    11201148            ax.set_title(f"Dust-to-Ice ratio over time (Grid point {ig+1}, Slope {isl+1})", fontweight='bold')
    11211149            ax.set_xlabel('Time (Mars years)')
    … …  
    11381166                    linewidth=1.5
    11391167                )
     1168            ax2.format_coord = format_coord_all
    11401169            ax2.set_ylabel('Obliquity (°)')
    11411170            ax2.tick_params(axis='y')

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

    r3842 r3851  
    722722dh_h2oice = d_h2oice/rho_h2oice
    723723dh_dust = 0.
    724 if (dh_h2oice >= 0.) then ! To put a dust sedimentation tendency only when ice is condensing
     724if (dh_h2oice > 0.) then ! To put a dust sedimentation tendency only when ice is condensing
    725725    if (impose_dust_ratio) then
    726         if (dh_co2ice >= 0.) then
     726        if (dh_co2ice > 0.) then
    727727            dh_dust = dust2ice_ratio*(dh_h2oice + dh_co2ice)
    728728        else

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

    r3850 r3851  
    291291! Loop variables
    292292integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap
     293real    :: totmass_ini
    293294logical :: num_str
    294295
    … …  
    995996
    996997    ! Checking mass balance for CO2
    997     totmass_co2ice = 0.
    998     totmass_atmco2 = 0.
    999     do ig = 1,ngrid
    1000         totmass_atmco2 = totmass_atmco2 + cell_area(ig)*ps_avg(ig)/g
    1001         do islope = 1,nslope
    1002             totmass_co2ice = totmass_co2ice + co2_ice(ig,islope)*cell_area(ig)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
    1003         enddo
    1004     enddo
    1005     write(*,'(a,f8.3,a)') " > Relative total CO2 mass balance = ", 100.*(totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini), ' %'
    1006     if ((totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini) > 0.01) then
    1007         write(*,*) '  /!\ Warning: mass balance is not conseved!'
    1008         write(*,'(a,f8.3,a)') '       Atmospheric CO2 mass balance = ', 100.*(totmass_atmco2 - totmass_atmco2_ini)/totmass_atmco2_ini, ' %'
    1009         write(*,'(a,f8.3,a)') '       CO2 ice mass balance         = ', 100.*(totmass_co2ice - totmass_co2ice_ini)/totmass_co2ice_ini, ' %'
    1010         write(*,'(a,f8.3,a)') '       Adsorbed CO2 mass balance    = ', 100.*(totmass_adsco2 - totmass_adsco2_ini)/totmass_adsco2_ini, ' %'
     998    if (abs(CO2cond_ps - 1.) < 1.e-10) then
     999        totmass_co2ice = 0.
     1000        totmass_atmco2 = 0.
     1001        do ig = 1,ngrid
     1002            totmass_atmco2 = totmass_atmco2 + cell_area(ig)*ps_avg(ig)/g
     1003            do islope = 1,nslope
     1004                totmass_co2ice = totmass_co2ice + co2_ice(ig,islope)*cell_area(ig)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
     1005            enddo
     1006        enddo
     1007        totmass_ini = max(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini,1.e-10)
     1008        write(*,'(a,f8.3,a)') " > Relative total CO2 mass balance = ", 100.*(totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/totmass_ini, ' %'
     1009        if ((totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/totmass_ini > 0.01) then
     1010            write(*,*) '  /!\ Warning: mass balance is not conseved!'
     1011            totmass_ini = max(totmass_atmco2_ini,1.e-10)
     1012            write(*,'(a,f8.3,a)') '       Atmospheric CO2 mass balance = ', 100.*(totmass_atmco2 - totmass_atmco2_ini)/totmass_ini, ' %'
     1013            totmass_ini = max(totmass_co2ice_ini,1.e-10)
     1014            write(*,'(a,f8.3,a)') '       CO2 ice mass balance         = ', 100.*(totmass_co2ice - totmass_co2ice_ini)/totmass_ini, ' %'
     1015            totmass_ini = max(totmass_adsco2_ini,1.e-10)
     1016            write(*,'(a,f8.3,a)') '       Adsorbed CO2 mass balance    = ', 100.*(totmass_adsco2 - totmass_adsco2_ini)/totmass_ini, ' %'
     1017        endif
    10111018    endif
    10121019

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

    r3850 r3851  
    6161Year = year_bp_ini + i_myear ! We compute the current year
    6262Lsp = lsperi*180./pi         ! We convert in degree
    63 
    64 print*, year_bp_ini,i_myear,i_myear_leg
    6563
    6664write(*,*) 'recomp_orb_param, Old year =', Year - i_myear_leg

  • trunk/LMDZ.MARS/util/display_netcdf.py

    r3849 r3851  
    8484
    8585# Attempt to load MOLA topography
    86 try:
    87     MOLA = np.load('MOLA_1px_per_deg.npy') # shape (nlat, nlon) at 1° per pixel: lat from -90 to 90, lon from 0 to 360
     86if os.path.isfile(MOLA_NPY): # shape (nlat, nlon) at 1° per pixel: lat from -90 to 90, lon from 0 to 360
     87    MOLA = np.load(MOLA_NPY)
    8888    nlat, nlon = MOLA.shape
    8989    topo_lats = np.linspace(90 - 0.5, -90 + 0.5, nlat)
    … …  
    9191    topo_lon2d, topo_lat2d = np.meshgrid(topo_lons, topo_lats)
    9292    topo_loaded = True
    93 except Exception as e:
    94     print(f"Warning: '{MOLA_NPY}' not found: {e}")
     93else:
     94    print(f"Warning: '{MOLA_NPY}' not found! Topography contours disabled.")
    9595    topo_loaded = False
    9696
    … …  
    101101    csv_loaded = True
    102102else:
    103     print(f"Warning: '{MOLA_CSV}' not found. 3D view colors disabled.")
     103    print(f"Warning: '{MOLA_CSV}' not found! 3D view colors disabled.")
    104104    csv_loaded = False
    105105
    … …  
    585585        ax.set_ylabel(f"Latitude ({getattr(dataset.variables[latv], 'units', 'deg')})")
    586586        # Prompt for polar-stereo views if interactive
    587         if input("Display polar-stereo views? [y/n]: ").strip().lower() == "y":
     587        if input("Display polar-stereo views? [y = yes, anything else = no]: ").strip().lower() == "y":
    588588            units = getattr(var, 'units', None)
    589589            plot_polar_views(lon2d, lat2d, grid, colormap, varname, units)
    590590        # Prompt for 3D globe view if interactive
    591         if input("Display 3D globe view? [y/n]: ").strip().lower() == "y":
     591        if input("Display 3D globe view? [y = yes, anything else = no]: ").strip().lower() == "y":
    592592            units = getattr(var, 'units', None)
    593593            plot_3D_globe(lon2d, lat2d, grid, colormap, varname, units)
    … …  
    659659        ax.grid(True)
    660660        # Prompt for polar-stereo views if interactive
    661         if sys.stdin.isatty() and input("Display polar-stereo views? [y/n]: ").strip().lower() == "y":
     661        if sys.stdin.isatty() and input("Display polar-stereo views? [y = yes, anything else = no]: ").strip().lower() == "y":
    662662            units = getattr(dataset.variables[varname], "units", None)
    663663            plot_polar_views(x_coords, y_coords, plot_data, colormap, varname, units)
    664664        # Prompt for 3D globe view if interactive
    665         if sys.stdin.isatty() and input("Display 3D globe view? [y/n]: ").strip().lower() == "y":
     665        if sys.stdin.isatty() and input("Display 3D globe view? [y = yes, anything else = no]: ").strip().lower() == "y":
    666666            units = getattr(dataset.variables[varname], "units", None)
    667667            plot_3D_globe(x_coords, y_coords, plot_data, colormap, varname, units)
Note: See TracChangeset for help on using the changeset viewer.