Changeset 3809 for trunk


Ignore:
Timestamp:
Jun 16, 2025, 5:11:50 PM (2 weeks ago)
Author:
jbclement
Message:

PEM:

  • Correction to detect whether the stratum below the dust lag is made of ice or not.
  • Correction for the initialization of surface ice and 'h2o_ice_depth' according to the layerings map.
  • Correction to handle flow of glaciers with the layering algorithm.
  • Simplification of "recomp_tend_h2o", making it more robust.
  • Making the removing of a stratum more robust.
  • Plotting the orbital variations over the simulated time in "visu_evol_layering.py".
  • Lowering the value of dust tendency.
  • Addition of "_co2" in the variables name for the CO2 flux correction.

JBC

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

Legend:

Unmodified
Added
Removed

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

    r3792 r3809  
    704704== 04/06/2025 == JBC
    705705Updates in the Python scripts to visualize the layering structure.
     706
     707== 16/06/2025 == JBC
     708- Correction to detect whether the stratum below the dust lag is made of ice or not.
     709- Correction for the initialization of surface ice and 'h2o_ice_depth' according to the layerings map.
     710- Correction to handle flow of glaciers with the layering algorithm.
     711- Simplification of "recomp_tend_h2o", making it more robust.
     712- Making the removing of a stratum more robust.
     713- Plotting the orbital variations over the simulated time in "visu_evol_layering.py".
     714- Lowering the value of dust tendency.
     715- Addition of "_co2" in the variables name for the CO2 flux correction.

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

    r3792 r3809  
    11#!/usr/bin/env python3
    22#######################################################################################################
    3 ### Python script to output the stratification data over time from the "restartpem#.nc" files files ###
     3### Python script to output stratification data over time from "restartpem#.nc" files               ###
     4### and to plot orbital parameters from "obl_ecc_lsp.asc"                                           ###
    45#######################################################################################################
    5 
    66
    77import os
    … …  
    4747        ).strip() or "info_PEM.txt"
    4848
    49     return folder_path, base_name, infofile
     49    orbfile = input(
     50        "Enter the name of the orbital parameters ASCII file "
     51        "(press Enter for default [obl_ecc_lsp.asc]): "
     52    ).strip() or "obl_ecc_lsp.asc"
     53    while not os.path.isfile(orbfile):
     54        print(f"  » \"{orbfile}\" does not exist or is not a file.")
     55        orbfile = input(
     56            "Enter a valid orbital parameters ASCII filename (press Enter for default [obl_ecc_lsp.asc]): "
     57        ).strip() or "info_PEM.txt"
     58
     59    return folder_path, base_name, infofile, orbfile
    5060
    5161
    … …  
    195205    """
    196206    Given raw_prop_arrays for 'co2_ice', 'h2o_ice', 'dust', 'pore' (in meters),
    197     normalize each set of strata so that the sum of those four = 1 per strata.
     207    normalize each set of strata so that the sum of those four = 1 per cell.
    198208    Returns:
    199209      - frac_arrays: dict mapping same keys -> 4D arrays of fractions (0..1).
    … …  
    228238    """
    229239    Reads "info_PEM.txt". Expects:
    230       - First line: parameters (ignored).
    231       - Each subsequent line: floats where first value is timestamp.
    232     Returns: 1D numpy array of timestamps.
     240      - First line: parameters where the 3rd value is martian_to_earth conversion factor.
     241      - Each subsequent line: floats where first value is simulation timestamp (in Mars years).
     242    Returns:
     243      - date_time: 1D numpy array of timestamps (Mars years)
     244      - martian_to_earth: float conversion factor
    233245    """
    234246    date_time = []
    235247    with open(file_name, 'r') as fp:
    236         fp.readline()
     248        first = fp.readline().split()
     249        martian_to_earth = float(first[2])
    237250        for line in fp:
    238251            parts = line.strip().split()
    … …  
    243256            except ValueError:
    244257                continue
    245     return np.array(date_time, dtype=np.float64)
     258    return np.array(date_time, dtype=np.float64), martian_to_earth
    246259
    247260
    … …  
    293306):
    294307    """
    295     Build a reference grid and interpolate strata fractions (0..1) onto it.
    296 
    297     Also returns a 'top_index' array of shape (ngrid, ntime, nslope) that
    298     indicates, for each (ig, t_idx, isl), the number of ref_grid levels
    299     covered by the topmost valid stratum.
    300 
    301     Args:
    302       - heights_data: list of lists where heights_data[t][isl] is a 2D array
    303           (ngrid, n_strata_current) of top_elevation values.
    304       - prop_arrays: dict mapping each property_name to a 4D array of shape
    305           (ngrid, ntime, nslope, max_nb_str) holding fractions [0..1].
    306       - min_base_for_interp: float; if exclude_sub=True, this is 0.0.
    307       - max_top_elev: float
    308       - dz: float
    309       - exclude_sub: bool. If True, ignore strata with elevation < 0.
     308    Build a reference elevation grid and interpolate strata fractions onto it.
    310309
    311310    Returns:
    312311      - ref_grid: 1D array of elevations (nz,)
    313       - gridded_data: dict mapping each property_name to a 4D array of shape
    314           (ngrid, ntime, nslope, nz) with interpolated fractions.
    315       - top_index: 3D array (ngrid, ntime, nslope) of ints: number of levels
    316           of ref_grid covered by the topmost stratum.
    317     """
    318     # Build ref_grid, ensuring at least two points if surface-only and dz > max_top_elev
     312      - gridded_data: dict mapping each property_name to 4D array
     313        (ngrid, ntime, nslope, nz) with interpolated fractions.
     314      - top_index: 3D array (ngrid, ntime, nslope) of ints:
     315        number of levels covered by the topmost stratum.
     316    """
    319317    if exclude_sub and (dz > max_top_elev):
    320318        ref_grid = np.array([0.0, max_top_elev], dtype=np.float32)
    321319    else:
    322         ref_grid = np.arange(min_base_for_interp, max_top_elev + dz / 2, dz)
     320        ref_grid = np.arange(min_base_for_interp, max_top_elev + dz/2, dz)
    323321    nz = len(ref_grid)
    324322    print(f"> Number of reference grid points = {nz}")
    325323
    326     # Dimensions
    327324    sample_prop = next(iter(prop_arrays.values()))
    328     ngrid, ntime, nslope, max_nb_str = sample_prop.shape[0:4]
    329 
    330     # Prepare outputs
     325    ngrid, ntime, nslope, max_nb_str = sample_prop.shape
     326
    331327    gridded_data = {
    332328        prop: np.full((ngrid, ntime, nslope, nz), -1.0, dtype=np.float32)
    … …  
    342338                    continue
    343339
    344                 raw_h = h_mat[ig, :]  # (n_strata_current,)
    345                 # Create h_all of length max_nb_str, fill with NaN
     340                raw_h = h_mat[ig, :]
    346341                h_all = np.full((max_nb_str,), np.nan, dtype=np.float32)
    347342                n_strata_current = raw_h.shape[0]
    … …  
    359354                h_valid = h_all[valid_mask]
    360355                top_h = np.max(h_valid)
    361 
    362                 # Find i_zmax = number of ref_grid levels z <= top_h
    363356                i_zmax = np.searchsorted(ref_grid, top_h, side='right')
    364357                top_index[ig, t_idx, isl] = i_zmax
    365 
    366358                if i_zmax == 0:
    367                     # top_h < ref_grid[0], skip interpolation
    368359                    continue
    369360
    370361                for prop, arr in prop_arrays.items():
    371                     prop_profile_all = arr[ig, t_idx, isl, :]    # (max_nb_str,)
    372                     prop_profile = prop_profile_all[valid_mask]  # (n_valid_strata,)
     362                    prop_profile_all = arr[ig, t_idx, isl, :]
     363                    prop_profile = prop_profile_all[valid_mask]
    373364                    if prop_profile.size == 0:
    374365                        continue
    375366
    376                     # Step‐wise interpolation (kind='next')
    377367                    f_interp = interp1d(
    378368                        h_valid,
    … …  
    382372                        fill_value=-1.0
    383373                    )
    384 
    385                     # Evaluate for ref_grid[0:i_zmax]
    386374                    gridded_data[prop][ig, t_idx, isl, :i_zmax] = f_interp(ref_grid[:i_zmax])
    387375
    … …  
    399387):
    400388    """
    401     For each grid point (ig) and slope (isl), generate a 2×2 figure:
     389    For each grid point and slope, generate a 2×2 figure of:
    402390      - CO2 ice fraction
    403391      - H2O ice fraction
    404392      - Dust fraction
    405393      - Pore fraction
    406 
    407     Fractions are in [0..1]. Values < 0 (fill) are masked.
    408     Using top_index, any elevation above the last stratum is forced to NaN (white).
    409 
    410     Additionally, draw horizontal violet line segments at each stratum top elevation
    411     only over the interval [date_time[t_idx], date_time[t_idx+1]] where that stratum
    412     exists at time t_idx. This way, boundaries appear only where the strata exist.
    413     """
    414     import numpy as np
    415     import matplotlib.pyplot as plt
    416 
     394    """
    417395    prop_names = ['co2_ice', 'h2o_ice', 'dust', 'pore']
    418396    titles = ["CO2 ice", "H2O ice", "Dust", "Pore"]
    … …  
    426404    if exclude_sub:
    427405        positive_indices = np.where(ref_grid >= 0.0)[0]
    428         if positive_indices.size == 0:
    429             print("Warning: no positive elevations in ref_grid → nothing to display.")
    430             return
    431406        sub_ref_grid = ref_grid[positive_indices]
    432407    else:
    … …  
    438413            fig, axes = plt.subplots(2, 2, figsize=(10, 8))
    439414            fig.suptitle(
    440                 f"Content variation over time for (Grid Point {ig+1}, Slope {isl+1})",
     415                f"Content variation over time for (Grid point {ig+1}, Slope {isl+1})",
    441416                fontsize=14
    442417            )
    443418
    444             # For each time step t_idx, gather this stratum's valid tops
    445             # and draw a line segment from date_time[t_idx] to date_time[t_idx+1].
    446             # We'll skip t_idx = ntime - 1 since no next point.
    447             # Precompute, for each t_idx, the array of valid top elevations:
     419            # Precompute valid stratum tops per time
    448420            valid_tops_per_time = []
    449421            for t_idx in range(ntime):
    450                 raw_h = heights_data[t_idx][isl][ig, :]  # (n_strata_current,)
    451                 # Exclude NaNs or zeros
     422                raw_h = heights_data[t_idx][isl][ig, :]
    452423                h_all = raw_h[~np.isnan(raw_h)]
    453424                if exclude_sub:
    … …  
    457428            for idx, prop in enumerate(prop_names):
    458429                ax = axes.flat[idx]
    459                 data_3d = gridded_data[prop][ig, :, isl, :]    # shape (ntime, nz)
    460                 mat_full = data_3d.T                           # shape (nz, ntime)
    461                 mat = mat_full[positive_indices, :].copy()     # (nz_pos, ntime)
    462 
    463                 # Mask fill values (< 0) as NaN
     430                data_3d = gridded_data[prop][ig, :, isl, :]
     431                mat_full = data_3d.T
     432                mat = mat_full[positive_indices, :].copy()
    464433                mat[mat < 0.0] = np.nan
    465434
    466                 # Mask everything above the top stratum using top_index
     435                # Mask above top stratum
    467436                for t_idx in range(ntime):
    468437                    i_zmax = top_index[ig, t_idx, isl]
    … …  
    473442                        mat[count_z:, t_idx] = np.nan
    474443
    475                 # Draw pcolormesh
    476444                im = ax.pcolormesh(
    477445                    date_time,
    … …  
    484452                )
    485453                ax.set_title(titles[idx], fontsize=12)
    486                 ax.set_xlabel("Time (y)")
     454                ax.set_xlabel("Time (Mars years)")
    487455                ax.set_ylabel("Elevation (m)")
    488456
    489                 # Draw horizontal violet segments only where strata exist
    490                 for t_idx in range(ntime - 1):
    491                     h_vals = valid_tops_per_time[t_idx]
    492                     if h_vals.size == 0:
    493                         continue
    494                     t_left = date_time[t_idx]
    495                     t_right = date_time[t_idx + 1]
    496                     for h in h_vals:
    497                         # Only draw if h falls within sub_ref_grid
    498                         if h < sub_ref_grid[0] or h > sub_ref_grid[-1]:
    499                             continue
    500                         ax.hlines(
    501                             y=h,
    502                             xmin=t_left,
    503                             xmax=t_right,
    504                             color='violet',
    505                             linewidth=0.7,
    506                             linestyle='-'
    507                         )
    508 
    509             # Reserve extra space on the right for the colorbar
    510457            fig.subplots_adjust(right=0.88)
    511 
    512             # Place a single shared colorbar in its own axes
    513458            cbar_ax = fig.add_axes([0.90, 0.15, 0.02, 0.7])
    514             fig.colorbar(
    515                 im,
    516                 cax=cbar_ax,
    517                 orientation='vertical',
    518                 label="Content"
    519             )
    520 
    521             # Tight layout excluding the region we reserved (0.88)
     459            fig.colorbar(im, cax=cbar_ax, orientation='vertical', label="Content")
    522460            fig.tight_layout(rect=[0, 0, 0.88, 1.0])
    523461
    … …  
    526464            )
    527465            fig.savefig(fname, dpi=150)
    528             plt.show()
    529             plt.close(fig)
     466
     467
     468def plot_strata_count_and_total_height(heights_data, date_time, output_folder="."):
     469    """
     470    For each grid point and slope, plot:
     471      - Number of strata vs time
     472      - Total deposit height vs time
     473    """
     474    ntime = len(heights_data)
     475    nslope = len(heights_data[0])
     476    ngrid = heights_data[0][0].shape[0]
     477
     478    for ig in range(ngrid):
     479        for isl in range(nslope):
     480            n_strata_t = np.zeros(ntime, dtype=int)
     481            total_height_t = np.zeros(ntime, dtype=float)
     482
     483            for t_idx in range(ntime):
     484                h_mat = heights_data[t_idx][isl]
     485                raw_h = h_mat[ig, :]
     486                valid_mask = (~np.isnan(raw_h)) & (raw_h != 0.0)
     487                if np.any(valid_mask):
     488                    h_valid = raw_h[valid_mask]
     489                    n_strata_t[t_idx] = h_valid.size
     490                    total_height_t[t_idx] = np.max(h_valid)
     491
     492            fig, axes = plt.subplots(2, 1, figsize=(8, 6), sharex=True)
     493            fig.suptitle(
     494                f"Strata count & total height over time for (Grid point {ig+1}, Slope {isl+1})",
     495                fontsize=14
     496            )
     497
     498            axes[0].plot(date_time, n_strata_t, marker='+', linestyle='-')
     499            axes[0].set_ylabel("Number of strata")
     500            axes[0].grid(True)
     501
     502            axes[1].plot(date_time, total_height_t, marker='+', linestyle='-')
     503            axes[1].set_xlabel("Time (Mars years)")
     504            axes[1].set_ylabel("Total height (m)")
     505            axes[1].grid(True)
     506
     507            fig.tight_layout(rect=[0, 0, 1, 0.95])
     508            fname = os.path.join(
     509                output_folder, f"strata_count_height_ig{ig+1}_is{isl+1}.png"
     510            )
     511            fig.savefig(fname, dpi=150)
     512
     513
     514def read_orbital_data(orb_file, martian_to_earth):
     515    """
     516    Read the .asc file containing obliquity, eccentricity and Ls p.
     517    Columns:
     518      0 = time in thousand Martian years
     519      1 = obliquity (deg)
     520      2 = eccentricity
     521      3 = Ls p (deg)
     522    Converts times to Earth years.
     523    """
     524    data = np.loadtxt(orb_file)
     525    dates_mka = data[:, 0]
     526    dates_yr = dates_mka * 1e3 / martian_to_earth
     527    obliquity = data[:, 1]
     528    eccentricity = data[:, 2]
     529    lsp = data[:, 3]
     530    return dates_yr, obliquity, eccentricity, lsp
     531
     532
     533def plot_orbital_parameters(infofile, orb_file, date_time, output_folder="."):
     534    """
     535    Plot the evolution of obliquity, eccentricity and Ls p
     536    versus simulated time.
     537    """
     538    # Read conversion factor from infofile
     539    _, martian_to_earth = read_infofile(infofile)
     540
     541    # Read orbital data
     542    dates_yr, obl, ecc, lsp = read_orbital_data(orb_file, martian_to_earth)
     543
     544    # Interpolate orbital parameters at simulation dates (date_time)
     545    obl_interp = interp1d(dates_yr, obl, kind='linear', bounds_error=False, fill_value="extrapolate")(date_time)
     546    ecc_interp = interp1d(dates_yr, ecc, kind='linear', bounds_error=False, fill_value="extrapolate")(date_time)
     547    lsp_interp = interp1d(dates_yr, lsp, kind='linear', bounds_error=False, fill_value="extrapolate")(date_time)
     548
     549    # Plot
     550    fig, axes = plt.subplots(3, 1, figsize=(8, 10), sharex=True)
     551    fig.suptitle("Orbital Parameters vs Simulated Time", fontsize=14)
     552
     553    axes[0].plot(date_time, obl_interp, 'r+', linestyle='-')
     554    axes[0].set_ylabel("Obliquity (°)")
     555    axes[0].grid(True)
     556
     557    axes[1].plot(date_time, ecc_interp, 'b+', linestyle='-')
     558    axes[1].set_ylabel("Eccentricity")
     559    axes[1].grid(True)
     560
     561    axes[2].plot(date_time, lsp_interp, 'g+', linestyle='-')
     562    axes[2].set_ylabel("Ls p (°)")
     563    axes[2].set_xlabel("Time (Mars years)")
     564    axes[2].grid(True)
     565
     566    plt.tight_layout(rect=[0, 0, 1, 0.96])
     567    fname = os.path.join(output_folder, "orbital_parameters.png")
     568    fig.savefig(fname, dpi=150)
    530569
    531570
    532571def main():
    533572    # 1) Get user inputs
    534     folder_path, base_name, infofile = get_user_inputs()
     573    folder_path, base_name, infofile, orbfile = get_user_inputs()
    535574
    536575    # 2) List and verify NetCDF files
    537576    files = list_netcdf_files(folder_path, base_name)
    538577    if not files:
    539         print(f"No NetCDF files named \"{base_name}#.nc\" found in \"{folder_path}\". Exiting.")
     578        print(f"No NetCDF files named \"{base_name}#.nc\" found in \"{folder_path}\".")
    540579        sys.exit(1)
    541     nfile = len(files)
    542     print(f"> Found {nfile} NetCDF file(s) matching \"{base_name}#.nc\".")
    543 
    544     # 3) Open one sample to get ngrid, nslope, lon/lat
     580    print(f"> Found {len(files)} NetCDF file(s).")
     581
     582    # 3) Open one sample to get grid dimensions & coordinates
    545583    sample_file = files[0]
    546584    ngrid, nslope, longitude, latitude = open_sample_dataset(sample_file)
    547     print(f"> ngrid  = {ngrid}")
    548     print(f"> nslope = {nslope}")
    549 
    550     # 4) Scan all files to collect variable info + global min/max elevations
    551     var_info, max_nb_str, min_base_elev, max_top_elev = collect_stratification_variables(
    552         files, base_name
    553     )
    554     print(f"> max(nb_str_max)      = {max_nb_str}")
    555     print(f"> min(base_elevation)  = {min_base_elev:.3f}")
    556     print(f"> max(top_elevation)   = {max_top_elev:.3f}")
    557 
    558     # 5) Open all datasets for extraction
     585    print(f"> ngrid  = {ngrid}, nslope = {nslope}")
     586
     587    # 4) Collect variable info + global min/max elevations
     588    var_info, max_nb_str, min_base_elev, max_top_elev = collect_stratification_variables(files, base_name)
     589    print(f"> max strata per slope = {max_nb_str}")
     590    print(f"> min base elev = {min_base_elev:.3f} m, max top elev = {max_top_elev:.3f} m")
     591
     592    # 5) Load full datasets
    559593    datasets = load_full_datasets(files)
    560594
    561     # 6) Extract raw stratification data
    562     heights_data, raw_prop_arrays, ntime = extract_stratification_data(
    563         datasets, var_info, ngrid, nslope, max_nb_str
    564     )
    565 
    566     # 7) Close all datasets
     595    # 6) Extract stratification data
     596    heights_data, raw_prop_arrays, ntime = extract_stratification_data(datasets, var_info, ngrid, nslope, max_nb_str)
     597
     598    # 7) Close datasets
    567599    for ds in datasets:
    568600        ds.close()
    569601
    570     # 8) Normalize raw prop arrays to volume fractions
     602    # 8) Normalize to fractions
    571603    frac_arrays = normalize_to_fractions(raw_prop_arrays)
    572604
    573     # 9) Ask whether to show subsurface
     605    # 9) Ask whether to include subsurface
    574606    show_subsurface = get_yes_no_input("Show subsurface layers?")
    575607    exclude_sub = not show_subsurface
    576 
    577608    if exclude_sub:
    578609        min_base_for_interp = 0.0
    579         print("> Will interpolate only elevations >= 0 m (surface strata).")
     610        print("> Interpolating only elevations >= 0 m (surface strata).")
    580611    else:
    581612        min_base_for_interp = min_base_elev
    582         print(f"> Will interpolate full depth (min base = {min_base_elev:.3f} m).")
    583 
    584     # 10) Prompt for discretization step
     613        print(f"> Interpolating full depth down to {min_base_elev:.3f} m.")
     614
     615    # 10) Prompt discretization step
    585616    dz = prompt_discretization_step(max_top_elev)
    586617
    587     # 11) Build reference grid and interpolate (returns top_index as well)
     618    # 11) Build reference grid and interpolate
    588619    ref_grid, gridded_data, top_index = interpolate_data_on_refgrid(
    589         heights_data,
    590         frac_arrays,
    591         min_base_for_interp,
    592         max_top_elev,
    593         dz,
    594         exclude_sub=exclude_sub
     620        heights_data, frac_arrays, min_base_for_interp, max_top_elev, dz, exclude_sub=exclude_sub
    595621    )
    596622
    597     # 12) Read time stamps from "info_PEM.txt"
    598     date_time = read_infofile(infofile)
     623    # 12) Read timestamps and conversion factor from infofile
     624    date_time, martian_to_earth = read_infofile(infofile)
    599625    if date_time.size != ntime:
    600         print(
    601             "Warning: number of timestamps does not match number of NetCDF files "
    602             f"({date_time.size} vs {ntime})."
    603         )
    604 
    605     # 13) Plot and save figures (passing top_index and heights_data)
     626        print(f"Warning: {date_time.size} timestamps vs {ntime} NetCDF files.")
     627
     628    # 13) Plot stratification over time
    606629    plot_stratification_over_time(
    607         gridded_data,
    608         ref_grid,
    609         top_index,
    610         heights_data,
    611         date_time,
    612         exclude_sub=exclude_sub,
    613         output_folder="."
     630        gridded_data, ref_grid, top_index, heights_data, date_time,
     631        exclude_sub=exclude_sub, output_folder="."
    614632    )
     633
     634    # 14) Plot strata count & total height vs time
     635    plot_strata_count_and_total_height(heights_data, date_time, output_folder=".")
     636
     637    # 15) Plot orbital parameters
     638    plot_orbital_parameters(infofile, orbfile, date_time, output_folder=".")
     639
     640    # 16) Show all figures
     641    plt.show()
    615642
    616643

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

    r3789 r3809  
    1919
    2020! Physical parameters
    21 real, parameter :: d_dust            = 5.78e-2 ! Tendency of dust [kg.m-2.y-1]
    22 real, parameter :: dry_lag_porosity  = 0.4     ! Porosity of dust lag
    23 real, parameter :: wet_lag_porosity  = 0.1     ! Porosity of water ice lag
    24 real, parameter :: regolith_porosity = 0.4     ! Porosity of regolith
    25 real, parameter :: breccia_porosity  = 0.1     ! Porosity of breccia
    26 real, parameter :: co2ice_porosity   = 0.      ! Porosity of CO2 ice
    27 real, parameter :: h2oice_porosity   = 0.      ! Porosity of H2O ice
    28 real, parameter :: rho_dust          = 2500.   ! Density of dust [kg.m-3]
    29 real, parameter :: rho_rock          = 3200.   ! Density of rock [kg.m-3]
    30 real, parameter :: h_patchy_dust     = 5.e-4   ! Height under which a dust lag is considered patchy [m]
    31 real, parameter :: h_patchy_ice      = 5.e-4   ! Height under which a H2O ice lag is considered patchy [m]
     21real, parameter :: d_dust            = 1.e-3 ! Tendency of dust [kg.m-2.y-1]
     22real, parameter :: dry_lag_porosity  = 0.4   ! Porosity of dust lag
     23real, parameter :: wet_lag_porosity  = 0.1   ! Porosity of water ice lag
     24real, parameter :: regolith_porosity = 0.4   ! Porosity of regolith
     25real, parameter :: breccia_porosity  = 0.1   ! Porosity of breccia
     26real, parameter :: co2ice_porosity   = 0.    ! Porosity of CO2 ice
     27real, parameter :: h2oice_porosity   = 0.    ! Porosity of H2O ice
     28real, parameter :: rho_dust          = 2500. ! Density of dust [kg.m-3]
     29real, parameter :: rho_rock          = 3200. ! Density of rock [kg.m-3]
     30real, parameter :: h_patchy_dust     = 5.e-4 ! Height under which a dust lag is considered patchy [m]
     31real, parameter :: h_patchy_ice      = 5.e-4 ! Height under which a H2O ice lag is considered patchy [m]
    3232
    3333! Parameters for CO2 flux correction (inspired by Levrard et al. 2007)
    34 real, parameter :: hmin_lag  = 0.5 ! Minimal height of the lag deposit to reduce the sublimation rate [m]
    35 real, parameter :: fred_subl = 0.1 ! Reduction factor of sublimation rate
     34real, parameter :: hmin_lag_co2 = 0.5 ! Minimal height of the lag deposit to reduce the sublimation rate [m]
     35real, parameter :: fred_sublco2 = 0.1 ! Reduction factor of sublimation rate
    3636
    3737! Stratum = node of the linked list
    … …  
    216216
    217217!---- Local variables
    218 type(stratum), pointer :: new_str
     218type(stratum), pointer :: new_str ! To make the argument point towards something
    219219
    220220!---- Code
    … …  
    225225this%nb_str = this%nb_str - 1
    226226
     227! Making the new links
     228if (associated(str%down)) then ! If there is a lower stratum
     229    str%down%up => str%up
     230else ! If it was the first stratum (bottom of layering)
     231    this%bottom => str%up
     232endif
     233if (associated(str%up)) then ! If there is an upper stratum
     234    str%up%down => str%down
     235else ! If it was the last stratum (top of layering)
     236    this%top => str%down
     237endif
     238
    227239! Remove the stratum from the layering
    228 str%down%up => str%up
    229 if (associated(str%up)) then ! If it is not the last one at the top of the layering
    230     str%up%down => str%down
    231 else
    232     this%top => str%down
    233 endif
    234240new_str => str%down
    235241nullify(str%up,str%down)
    … …  
    550556!=======================================================================
    551557! To get data about possible subsurface ice in a layering
    552 SUBROUTINE subsurface_ice_layering(this,h_toplag,h2o_ice)
     558SUBROUTINE subsurface_ice_layering(this,h_toplag,h2o_ice,co2_ice)
    553559
    554560implicit none
    … …  
    556562!---- Arguments
    557563type(layering), intent(in) :: this
    558 real,    intent(out) :: h2o_ice, h_toplag
     564real, intent(out) :: h2o_ice, co2_ice, h_toplag
    559565
    560566!---- Local variables
    … …  
    564570h_toplag = 0.
    565571h2o_ice = 0.
     572co2_ice = 0.
    566573current => this%top
    567574! Is the considered stratum a dust lag?
    568 if (is_dust_lag(current)) return
     575if (.not. is_dust_lag(current)) return
    569576do while (is_dust_lag(current) .and. associated(current))
    570577    h_toplag = h_toplag + thickness(current)
    … …  
    572579enddo
    573580! Is the stratum below the dust lag made of ice?
    574 if (.not. is_h2oice_str(current) .or. h_toplag < h_patchy_dust) then
    575     h_toplag = 0.
    576     h2o_ice = current%h_h2oice
     581if (current%h_h2oice > 0. .and. current%h_h2oice > current%h_co2ice) then ! There is more H2O ice than CO2 ice
     582    if (h_toplag < h_patchy_dust) then ! But the dust lag is too thin to considered ice as subsurface ice
     583        h_toplag = 0.
     584        h2o_ice = current%h_h2oice
     585    endif
     586else if (current%h_co2ice > 0. .and. current%h_co2ice > current%h_h2oice) then ! There is more CO2 ice than H2O ice
     587    h_toplag = 0. ! Because it matters only for H2O ice
     588    if (h_toplag < h_patchy_ice) then ! But the dust lag is too thin to considered ice as subsurface ice
     589        co2_ice = current%h_co2ice
     590    endif
    577591endif
    578592
    … …  
    759773                else if (0. < currentloc%h_h2oice .and. currentloc%h_h2oice < h_patchy_ice) then
    760774                    h_toplag = h_toplag + thickness(currentloc)
    761                     R_subl = fred_subl**(h_toplag/hmin_lag)
     775                    R_subl = fred_sublco2**(h_toplag/hmin_lag_co2)
    762776                else
    763                     R_subl = fred_subl**(h_toplag/hmin_lag)
     777                    R_subl = fred_sublco2**(h_toplag/hmin_lag_co2)
    764778                endif
    765779            endif

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

    r3789 r3809  
    670670    do ig = 1,ngrid
    671671        do islope = 1,nslope
    672             co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice
    673             h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice
     672            if (is_co2ice_str(layerings_map(ig,islope)%top)) then
     673                co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice
     674            else if (is_h2oice_str(layerings_map(ig,islope)%top)) then
     675                h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice
     676            else
     677                call subsurface_ice_layering(layerings_map(ig,islope),h2o_ice_depth(ig,islope),h2o_ice(ig,islope),co2_ice(ig,islope))
     678            endif
    674679        enddo
    675680    enddo
    … …  
    846851                    h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice
    847852                else
    848                     call subsurface_ice_layering(layerings_map(ig,islope),h2o_ice_depth(ig,islope),h2o_ice(ig,islope))
     853                    call subsurface_ice_layering(layerings_map(ig,islope),h2o_ice_depth(ig,islope),h2o_ice(ig,islope),co2_ice(ig,islope))
    849854                endif
    850855            enddo
    851856        enddo
    852         call print_layering(layerings_map(1,1))
    853857    else
    854858        zlag = 0.
    … …  
    862866!------------------------
    863867    allocate(flag_co2flow(ngrid,nslope),flag_h2oflow(ngrid,nslope))
    864     if (co2ice_flow .and. nslope > 1) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys, &
    865                                                             ps_timeseries,ps_avg_global_old,ps_avg_global_new,co2_ice,flag_co2flow)
    866     if (h2oice_flow .and. nslope > 1) call flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow)
    867     if (layering_algo) then
    868         do islope = 1,nslope
    869             do ig = 1,ngrid
    870                 layerings_map(ig,islope)%top%h_co2ice = co2_ice(ig,islope)
    871                 layerings_map(ig,islope)%top%h_h2oice = h2o_ice(ig,islope)
     868    if (co2ice_flow .and. nslope > 1) then
     869        call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys, &
     870                              ps_timeseries,ps_avg_global_old,ps_avg_global_new,co2_ice,flag_co2flow)
     871        if (layering_algo) then
     872            do islope = 1,nslope
     873                do ig = 1,ngrid
     874                    layerings_map(ig,islope)%top%h_co2ice = co2_ice(ig,islope)
     875                enddo
    872876            enddo
    873         enddo
     877        endif
     878    endif
     879    if (h2oice_flow .and. nslope > 1) then
     880        call flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow)
     881        if (layering_algo) then
     882            do islope = 1,nslope
     883                do ig = 1,ngrid
     884!~                     layerings_map(ig,islope)%top%h_h2oice = h2o_ice(ig,islope)
     885                enddo
     886            enddo
     887        endif
    874888    endif
    875889
    … …  
    10781092    if (soil_pem) deallocate(icetable_depth_old,tsoil_PEM_timeseries_old)
    10791093    deallocate(vmr_co2_PEM_phys)
    1080     write(*,*) "> Updating the H2O sub-surface ice depth"
    10811094
    10821095!------------------------

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

    r3789 r3809  
    9292Rz_old = h2oice_depth_old*zcdv/coef_diff ! Old resistance from PCM
    9393Rz_new = h2oice_depth_new*zcdv/coef_diff ! New resistance based on new depth
    94 R_dec = 1.
    95 if (Rz_new >= Rz_old) R_dec = Rz_new/Rz_old ! Decrease because of resistance
     94R_dec = Rz_old/Rz_new ! Decrease because of resistance
    9695
    9796! The maxmimum of the daily averages over one year for the saturation vapor pressure at the ice table location
    … …  
    104103
    105104! Lower humidity due to growing lag layer (higher depth)
    106 hum_dec = psv_max_old/psv_max_new ! Decrease because of lower water vapor pressure at the new depth
     105if (abs(psv_max_old) < 1.e2*epsilon(1.)) then
     106    hum_dec = 1.
     107else
     108    hum_dec = psv_max_new/psv_max_old ! Decrease because of lower water vapor pressure at the new depth
     109endif
    107110
    108111! Flux correction (decrease)
    109 d_h2oice = d_h2oice/R_dec/hum_dec
     112d_h2oice = d_h2oice*R_dec*hum_dec
    110113
    111114END SUBROUTINE recomp_tend_h2o
Note: See TracChangeset for help on using the changeset viewer.