Ignore:
Timestamp:
Dec 19, 2023, 3:14:48 PM (12 months ago)
Author:
jbclement
Message:

PEM:

  • Addition of flags defined in the "run_PEM.def" to decide to do or not CO2 & H2O ice metamorphism: 'metam_co2ice' and 'metam_h2oice' (default is false).
  • The variations of infinite reservoirs ('watercap') during the PCM years are now taken into account to update H2O ice at the PEM initialization.
  • 'ini_h2o_bigreservoir' is renamed into 'ini_huge_h2oice'.
  • Some cleanings, in particular for the main program "pem.F90".

JBC

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

Legend:

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

    r3159 r3161  
    179179- Addition of (CO2 and H2O) ice metamorphism: if the frost is above a given value, then the excess frost is transferred to the perennial ice.
    180180- Thresholds value for ice management can now be set in the "run_PEM.def".
     181
     182== 19/12/2023 == JBC
     183PEM:
     184- Addition of flags defined in the "run_PEM.def" to decide to do or not CO2 & H2O ice metamorphism: 'metam_co2ice' and 'metam_h2oice' (default is false).
     185- The variations of infinite reservoirs ('watercap') during the PCM years are now taken into account to update H2O ice at the PEM initialization.
     186- 'ini_h2o_bigreservoir' is renamed into 'ini_huge_h2oice'.
     187- Some cleanings, in particular for the main program "pem.F90".
  • trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90

    r3149 r3161  
    33implicit none
    44
    5 integer, parameter                        :: nsoilmx_PEM = 68     ! number of layers in the PEM
    6 real, allocatable, dimension(:),     save :: layer_PEM            ! soil layer depths [m]
    7 real, allocatable, dimension(:),     save :: mlayer_PEM           ! soil mid-layer depths [m]
    8 real, allocatable, dimension(:,:,:), save :: TI_PEM               ! soil thermal inertia [SI]
    9 real, allocatable, dimension(:,:),   save :: inertiedat_PEM       ! soil thermal inertia saved as reference for current climate [SI]
     5integer, parameter                        :: nsoilmx_PEM = 68   ! number of layers in the PEM
     6real, allocatable, dimension(:),     save :: layer_PEM          ! soil layer depths [m]
     7real, allocatable, dimension(:),     save :: mlayer_PEM         ! soil mid-layer depths [m]
     8real, allocatable, dimension(:,:,:), save :: TI_PEM             ! soil thermal inertia [SI]
     9real, allocatable, dimension(:,:),   save :: inertiedat_PEM     ! soil thermal inertia saved as reference for current climate [SI]
    1010! Variables (FC: built in firstcall in soil.F)
    11 real, allocatable, dimension(:,:,:), save :: tsoil_PEM            ! sub-surface temperatures [K]
    12 real, allocatable, dimension(:,:),   save :: mthermdiff_PEM       ! (FC) mid-layer thermal diffusivity [SI]
    13 real, allocatable, dimension(:,:),   save :: thermdiff_PEM        ! (FC) inter-layer thermal diffusivity [SI]
    14 real, allocatable, dimension(:),     save :: coefq_PEM            ! (FC) q_{k+1/2} coefficients [SI]
    15 real, allocatable, dimension(:,:),   save :: coefd_PEM            ! (FC) d_k coefficients [SI]
    16 real, allocatable, dimension(:,:),   save :: alph_PEM             ! (FC) alpha_k coefficients [SI]
    17 real, allocatable, dimension(:,:),   save :: beta_PEM             ! beta_k coefficients [SI]
    18 real,                                save :: mu_PEM               ! mu coefficient [SI]
    19 real,                                save :: fluxgeo              ! Geothermal flux [W/m^2]
    20 real,                                save :: depth_breccia        ! Depth at which we have breccia [m]
    21 real,                                save :: depth_bedrock        ! Depth at which we have bedrock [m]
    22 integer,                             save :: index_breccia        ! last index of the depth grid before having breccia
    23 integer,                             save :: index_bedrock        ! last index of the depth grid before having bedrock
    24 logical                                   :: soil_pem             ! True by default, to run with the subsurface physic. Read in pem.def
    25 real,                                save :: ini_h2o_bigreservoir ! Initial value for the large reservoir of water [kg/m^2]
    26 logical,                             save :: reg_thprop_dependp   ! thermal properites of the regolith vary with the surface pressure
     11real, allocatable, dimension(:,:,:), save :: tsoil_PEM          ! sub-surface temperatures [K]
     12real, allocatable, dimension(:,:),   save :: mthermdiff_PEM     ! (FC) mid-layer thermal diffusivity [SI]
     13real, allocatable, dimension(:,:),   save :: thermdiff_PEM      ! (FC) inter-layer thermal diffusivity [SI]
     14real, allocatable, dimension(:),     save :: coefq_PEM          ! (FC) q_{k+1/2} coefficients [SI]
     15real, allocatable, dimension(:,:),   save :: coefd_PEM          ! (FC) d_k coefficients [SI]
     16real, allocatable, dimension(:,:),   save :: alph_PEM           ! (FC) alpha_k coefficients [SI]
     17real, allocatable, dimension(:,:),   save :: beta_PEM           ! beta_k coefficients [SI]
     18real,                                save :: mu_PEM             ! mu coefficient [SI]
     19real,                                save :: fluxgeo            ! Geothermal flux [W/m^2]
     20real,                                save :: depth_breccia      ! Depth at which we have breccia [m]
     21real,                                save :: depth_bedrock      ! Depth at which we have bedrock [m]
     22integer,                             save :: index_breccia      ! last index of the depth grid before having breccia
     23integer,                             save :: index_bedrock      ! last index of the depth grid before having bedrock
     24logical                                   :: soil_pem           ! True by default, to run with the subsurface physic. Read in pem.def
     25real,                                save :: ini_huge_h2oice    ! Initial value for huge reservoirs of H2O ice [kg/m^2]
     26logical,                             save :: reg_thprop_dependp ! thermal properites of the regolith vary with the surface pressure
    2727
    2828!=======================================================================
  • trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90

    r3159 r3161  
    2525use time_evol_mod,         only: year_bp_ini, dt_pem, h2o_ice_crit, co2_ice_crit, ps_criterion, Max_iter_pem, &
    2626                          evol_orbit_pem, var_obl, var_ecc, var_lsp, convert_years
    27 use comsoil_h_pem,         only: soil_pem, fluxgeo, ini_h2o_bigreservoir, depth_breccia, depth_bedrock, reg_thprop_dependp
     27use comsoil_h_pem,         only: soil_pem, fluxgeo, ini_huge_h2oice, depth_breccia, depth_bedrock, reg_thprop_dependp
    2828use adsorption_mod,        only: adsorption_pem
    29 use glaciers_mod,          only: co2_ice_flow, h2o_ice_flow
     29use glaciers_mod,          only: h2oice_flow, co2ice_flow, inf_h2oice_threshold, metam_co2ice_threshold, metam_h2oice_threshold, metam_h2oice, metam_co2ice
    3030use ice_table_mod,         only: icetable_equilibrium, icetable_dynamic
    3131use time_phylmdz_mod,      only: ecritphy
    32 use constants_marspem_mod, only: inf_h2oice_threshold, metam_co2ice_threshold, metam_h2oice_threshold
    3332
    3433implicit none
     
    101100call getin('adsorption_pem',adsorption_pem)
    102101
    103 co2_ice_flow = .true.
    104 call getin('co2_ice_flow',co2_ice_flow)
    105 
    106 h2o_ice_flow = .true.
    107 call getin('h2o_ice_flow',h2o_ice_flow)
    108 
    109102reg_thprop_dependp = .false.
    110103call getin('reg_thprop_dependp',reg_thprop_dependp)
     
    158151
    159152!#---------- Ice management parameters ----------#
    160 ini_h2o_bigreservoir = 1.e4
    161 call getin('ini_h2o_bigreservoir',ini_h2o_bigreservoir)
     153ini_huge_h2oice = 1.e4
     154call getin('ini_huge_h2oice',ini_huge_h2oice)
    162155
    163156inf_h2oice_threshold = 2.e3
    164157call getin('inf_h2oice_threshold',inf_h2oice_threshold)
    165158
     159metam_h2oice = .false.
     160call getin('metam_h2oice',metam_h2oice)
     161
    166162metam_h2oice_threshold = 5.e-2
    167163call getin('metam_h2oice_threshold',metam_h2oice_threshold)
     164
     165h2oice_flow = .true.
     166call getin('h2oice_flow',h2oice_flow)
     167
     168metam_co2ice = .false.
     169call getin('metam_co2ice',metam_co2ice)
    168170
    169171metam_co2ice_threshold = 16.e3
    170172call getin('metam_co2ice_threshold',metam_co2ice_threshold)
    171173
     174co2ice_flow = .true.
     175call getin('co2ice_flow',co2ice_flow)
     176
    172177END SUBROUTINE conf_pem
    173178
  • trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90

    r3159 r3161  
    3737real, parameter :: Lco2 =  5.71e5 ! Pilorget and Forget 2016
    3838
    39 ! Thresholds for ice management
    40 real, save :: inf_h2oice_threshold   ! To consider the amount of H2O ice as an infinite reservoir
    41 real, save :: metam_h2oice_threshold ! To consider frost is becoming perennial H2O ice
    42 real, save :: metam_co2ice_threshold ! To consider frost is becoming perennial CO2 ice
    43 
    4439END MODULE constants_marspem_mod
  • trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def

    r3159 r3161  
    2828#  Max_iter_pem=100000000
    2929
    30 # Acceptance rate of sublimating H2o ice surface change? Default = 0.2
     30# Acceptance rate of sublimating H2O ice surface change? Default = 0.2
    3131# h2o_ice_crit=0.2
    3232
     
    4747# adsorption_pem=.true.
    4848
    49 # Do you want the glaciers to flow along subslope inside a mesh? Default = .true.
    50 # co2glaciersflow=.true.
    51 
    5249# Do you want to modify the soil thermal properties with the pressure? Default = .false.
    5350# reg_thprop_dependp=.false.
    5451
    5552#---------- Layering parameters ----------#
    56 # Value of the geothermal flux? Default= 0.
     53# Value of the geothermal flux? Default = 0.
    5754# fluxgeo=0.
    5855
     
    7067
    7168#---------- Ice management parameters ----------#
    72 # Amount of h2o ice to initialize the big reservoir if the variable is not present in "startfi_PEM.nc"? Default = 1.e4
    73 # ini_h2o_bigreservoir=1.e4
     69# Amount of H2O ice to initialize the huge reservoir if the variable is not present in "startfi_PEM.nc"? Default = 1.e4
     70# ini_huge_h2oice=1.e4
    7471
    7572# Threshold to consider the amount of H2O ice as an infinite reservoir? Default = 2.e3
    7673# inf_h2oice_threshold=2.e3
    7774
     75# Do you want H2O frost to transform into perennial H2O ice? Default = .false.
     76# metam_h2oice=.false.
     77
    7878# Threshold to consider frost is becoming perennial H2O ice? Default = 5.e-2
    7979# metam_h2oice_threshold=5.e-2
     80
     81# Do you want the H2O ice to flow along subslope inside a cell? Default = .true.
     82# h2oice_flow=.true.
     83
     84# Do you want CO2 frost to transform into perennial CO2 ice? Default = .false.
     85# metam_co2ice=.false.
    8086
    8187# Threshold to consider frost is becoming perennial CO2 ice? Default = 16.e3
    8288# metam_co2ice_threshold=16.e3
    8389
     90# Do you want the CO2 ice to flow along subslope inside a cell? Default = .true.
     91# co2ice_flow=.true.
     92
    8493# Some definitions for the physics, in file 'callphys.def'
    8594INCLUDEDEF=callphys.def
  • trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90

    r3149 r3161  
    33implicit none
    44
    5 logical :: co2_ice_flow ! True by default, to compute co2 ice flow. Read in run_PEM.def
    6 logical :: h2o_ice_flow ! True by default, to compute h2o ice flow. Read in run_PEM.def
     5
     6! Flags for ice management
     7logical :: h2oice_flow ! True by default, to compute H2O ice flow. Read in "run_PEM.def"
     8logical :: co2ice_flow ! True by default, to compute CO2 ice flow. Read in "run_PEM.def"
     9logical :: metam_h2oice ! False by default, to compute H2O ice metamorphism. Read in "run_PEM.def"
     10logical :: metam_co2ice ! False by default, to compute CO2 ice metamorphism. Read in "run_PEM.def"
     11
     12! Thresholds for ice management
     13real, save :: inf_h2oice_threshold   ! To consider the amount of H2O ice as an infinite reservoir
     14real, save :: metam_h2oice_threshold ! To consider frost is becoming perennial H2O ice
     15real, save :: metam_co2ice_threshold ! To consider frost is becoming perennial CO2 ice
    716
    817!=======================================================================
  • trunk/LMDZ.COMMON/libf/evolution/pem.F90

    r3159 r3161  
    11!------------------------
    22! I   Initialization
    3 !    I_a Read run.def
    4 !    I_b Read of start_evol.nc and starfi_evol.nc
     3!    I_a Read the "run.def"
     4!    I_b Read the "start_evol.nc" and "startfi_evol.nc"
    55!    I_c Subslope parametrisation
    6 !    I_d Read PCM data and convert to the physical grid
     6!    I_d Read the PCM data and convert them to the physical grid
    77!    I_e Initialization of the PEM variable and soil
    88!    I_f Compute tendencies
    99!    I_g Save initial PCM situation
    10 !    I_h Read the startpem.nc
     10!    I_h Read the "startpem.nc"
    1111!    I_i Compute orbit criterion
    1212
     
    2222! III Output
    2323!    III_a Update surface value for the PCM start files
    24 !    III_b Write restart_evol.nc and restartfi_evol.nc
    25 !    III_c Write restartpem.nc
     24!    III_b Write the "restart_evol.nc" and "restartfi_evol.nc"
     25!    III_c Write the "restartpem.nc"
    2626!------------------------
    2727
    2828PROGRAM pem
    2929
    30 use phyetat0_mod,                only: phyetat0
    31 use phyredem,                    only: physdem0, physdem1
    32 use netcdf,                      only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
    33 use turb_mod,                    only: q2, wstar
    34 use comslope_mod,                only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, major_slope, ini_comslope_h
    35 use logic_mod,                   only: iflag_phys
    36 use mod_const_mpi,               only: COMM_LMDZ
     30use phyetat0_mod,               only: phyetat0
     31use phyredem,                   only: physdem0, physdem1
     32use netcdf,                     only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
     33use turb_mod,                   only: q2, wstar
     34use comslope_mod,               only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, major_slope, ini_comslope_h
     35use logic_mod,                  only: iflag_phys
     36use mod_const_mpi,              only: COMM_LMDZ
    3737use infotrac
    38 use geometry_mod,                only: latitude_deg
    39 use conf_pem_mod,                only: conf_pem
    40 use pemredem,                    only: pemdem0, pemdem1
    41 use glaciers_mod,                only: flow_co2glaciers, flow_h2oglaciers, co2_ice_flow, h2o_ice_flow
    42 use stopping_crit_mod,           only: stopping_crit_h2o_ice, stopping_crit_co2
    43 use constants_marspem_mod,       only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, &
    44                                        m_noco2, inf_h2oice_threshold, metam_co2ice_threshold, metam_h2oice_threshold
    45 use evol_ice_mod,                only: evol_co2_ice, evol_h2o_ice
    46 use comsoil_h_PEM,               only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
    47                                        TI_PEM, inertiedat_PEM,           & ! soil thermal inertia
    48                                        tsoil_PEM, mlayer_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
    49                                        fluxgeo                             ! Geothermal flux for the PEM and PCM
    50 use adsorption_mod,              only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
    51                                        ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
    52                                        co2_adsorbded_phys, h2o_adsorbded_phys        ! Mass of co2 and h2O adsorbded
    53 use time_evol_mod,               only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
    54 use orbit_param_criterion_mod,   only: orbit_param_criterion
    55 use recomp_orb_param_mod,        only: recomp_orb_param
    56 use ice_table_mod,               only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
    57                                        ini_ice_table_porefilling, computeice_table_equilibrium,compute_massh2o_exchange_ssi
    58 use soil_thermalproperties_mod,  only: update_soil_thermalproperties
    59 use time_phylmdz_mod,            only: daysec, dtphys, day_end
    60 use abort_pem_mod,               only: abort_pem
    61 use soil_settings_PEM_mod,       only: soil_settings_PEM
    62 use compute_tend_mod,            only: compute_tend
    63 use info_PEM_mod,                only: info_PEM
    64 use interpol_TI_PEM2PCM_mod, only: interpol_TI_PEM2PCM
    65 use nb_time_step_PCM_mod,        only: nb_time_step_PCM
    66 use pemetat0_mod,                only: pemetat0
    67 use read_data_PCM_mod,           only: read_data_PCM
    68 use recomp_tend_co2_slope_mod,   only: recomp_tend_co2_slope
    69 use soil_pem_compute_mod,        only: soil_pem_compute
    70 use writediagpem_mod,            only: writediagpem
     38use geometry_mod,               only: latitude_deg
     39use conf_pem_mod,               only: conf_pem
     40use pemredem,                   only: pemdem0, pemdem1
     41use glaciers_mod,               only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, inf_h2oice_threshold, &
     42                                      metam_h2oice_threshold, metam_co2ice_threshold, metam_h2oice, metam_co2ice
     43use stopping_crit_mod,          only: stopping_crit_h2o_ice, stopping_crit_co2
     44use constants_marspem_mod,      only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2
     45use evol_ice_mod,               only: evol_co2_ice, evol_h2o_ice
     46use comsoil_h_PEM,              only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
     47                                      TI_PEM, inertiedat_PEM,           & ! soil thermal inertia
     48                                      tsoil_PEM, mlayer_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
     49                                      fluxgeo                             ! Geothermal flux for the PEM and PCM
     50use adsorption_mod,             only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
     51                                      ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
     52                                      co2_adsorbded_phys, h2o_adsorbded_phys        ! Mass of co2 and h2O adsorbded
     53use time_evol_mod,              only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
     54use orbit_param_criterion_mod,  only: orbit_param_criterion
     55use recomp_orb_param_mod,       only: recomp_orb_param
     56use ice_table_mod,              only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
     57                                      ini_ice_table_porefilling, computeice_table_equilibrium,compute_massh2o_exchange_ssi
     58use soil_thermalproperties_mod, only: update_soil_thermalproperties
     59use time_phylmdz_mod,           only: daysec, dtphys, day_end
     60use abort_pem_mod,              only: abort_pem
     61use soil_settings_PEM_mod,      only: soil_settings_PEM
     62use compute_tend_mod,           only: compute_tend
     63use info_PEM_mod,               only: info_PEM
     64use interpol_TI_PEM2PCM_mod,    only: interpol_TI_PEM2PCM
     65use nb_time_step_PCM_mod,       only: nb_time_step_PCM
     66use pemetat0_mod,               only: pemetat0
     67use read_data_PCM_mod,          only: read_data_PCM
     68use recomp_tend_co2_slope_mod,  only: recomp_tend_co2_slope
     69use soil_pem_compute_mod,       only: soil_pem_compute
     70use writediagpem_mod,           only: writediagpem
    7171
    7272#ifndef CPP_STD
     
    291291!------------------------
    292292! I   Initialization
    293 !    I_a Read run.def
     293!    I_a Read the "run.def"
    294294!------------------------
    295295#ifndef CPP_1D
     
    328328!------------------------
    329329! I   Initialization
    330 !    I_b Read of start_evol.nc and starfi_evol.nc
    331 !------------------------
    332 ! I_b.1 Read start_evol.nc
     330!    I_b Read of the "start_evol.nc" and starfi_evol.nc
     331!------------------------
     332! I_b.1 Read "start_evol.nc"
    333333allocate(ps_start_PCM(ngrid))
    334334#ifndef CPP_1D
     
    355355
    356356! In the PCM, these values are given to the physic by the dynamic.
    357 ! Here we simply read them in the startfi_evol.nc file
     357! Here we simply read them in the "startfi_evol.nc" file
    358358status = nf90_open(startfi_name, NF90_NOWRITE, ncid)
    359359
     
    372372status = nf90_close(ncid)
    373373
    374 ! I_b.2 Read startfi_evol.nc
     374! I_b.2 Read the "startfi_evol.nc"
    375375! First we read the initial state (starfi.nc)
    376376#ifndef CPP_STD
     
    468468!------------------------
    469469! I   Initialization
    470 !    I_d Read PCM data and convert to the physical grid
     470!    I_d Read the PCM data and convert them to the physical grid
    471471!------------------------
    472472! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the PCM run, saving only the minimum value
     
    579579!------------------------
    580580! I   Initialization
    581 !    I_h Read the startpem.nc
     581!    I_h Read the "startpem.nc"
    582582!------------------------
    583583allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope))
     
    655655    write(*,*) 'Global average pressure new time step',global_avg_press_new
    656656
    657 ! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
     657! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consumption)
    658658    allocate(zplev_old_timeseries(ngrid,nlayer + 1,timelen))
    659659    write(*,*) "Recomputing the old pressure levels timeserie adapted to the old pressure..."
     
    725725!    II_c Flow of glaciers
    726726!------------------------
    727     if (co2_ice_flow) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys,ps_timeseries, &
     727    if (co2ice_flow) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys,ps_timeseries, &
    728728                                            global_avg_press_PCM,global_avg_press_new,co2_ice,flag_co2flow,flag_co2flow_mesh)
    729     if (h2o_ice_flow) call flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,h2o_ice,flag_h2oflow,flag_h2oflow_mesh)
     729    if (h2oice_flow) call flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,h2o_ice,flag_h2oflow,flag_h2oflow_mesh)
    730730
    731731!------------------------
     
    944944do ig = 1,ngrid
    945945    ! H2O ice metamorphism
    946     if (sum(qsurf(ig,igcm_h2o_ice,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_h2oice_threshold) then
     946    if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_h2oice_threshold) then
    947947        h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
    948948        qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
     
    961961
    962962    ! CO2 ice metamorphism
    963     if (sum(qsurf(ig,igcm_co2,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_co2ice_threshold) then
     963    if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_co2ice_threshold) then
    964964        perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
    965965        qsurf(ig,igcm_co2,:) = metam_co2ice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
     
    10221022!------------------------
    10231023! III Output
    1024 !    III_b Write restart_evol.nc and restartfi_evol.nc
    1025 !------------------------
    1026 ! III_b.1 Write restart_evol.nc
     1024!    III_b Write "restart_evol.nc" and "restartfi_evol.nc"
     1025!------------------------
     1026! III_b.1 Write "restart_evol.nc"
    10271027ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
    10281028pday = day_ini
     
    10411041#endif
    10421042
    1043 ! III_b.2 Write restartfi_evol.nc
     1043! III_b.2 Write the "restartfi_evol.nc"
    10441044#ifndef CPP_STD
    10451045    call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
     
    10631063!------------------------
    10641064! III Output
    1065 !    III_c Write restartpem.nc
     1065!    III_c Write the "restartpem.nc"
    10661066!------------------------
    10671067call pemdem0("restartpem.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, &
  • trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90

    r3159 r3161  
    1313
    1414use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var
    15 use comsoil_h_PEM,              only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, ini_h2o_bigreservoir, depth_breccia, depth_bedrock, index_breccia, index_bedrock
     15use comsoil_h_PEM,              only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, ini_huge_h2oice, depth_breccia, depth_bedrock, index_breccia, index_bedrock
    1616use comsoil_h,                  only: volcapa, inertiedat
    1717use adsorption_mod,             only: regolith_adsorption, adsorption_pem
     
    2323use soil_pem_compute_mod,       only: soil_pem_compute
    2424use soil_pem_ini_mod,           only: soil_pem_ini
     25use comslope_mod,               only: def_slope_mean, subslope_dist
    2526
    2627#ifndef CPP_STD
    27     use comcstfi_h,   only: r, mugaz
    28     use surfdat_h,    only: watercaptag, perennial_co2ice
     28    use comcstfi_h,   only: r, mugaz, pi
     29    use surfdat_h,    only: watercaptag, watercap, perennial_co2ice
    2930#else
    30     use comcstfi_mod, only: r, mugaz
     31    use comcstfi_mod, only: r, mugaz, pi
    3132#endif
    3233
     
    120121        write(*,*)'Pemetat0: failed loading <h2o_ice>'
    121122        write(*,*)'will reconstruct the values from watercaptag'
    122         write(*,*)'with default value ''ini_h2o_bigreservoir'''
     123        write(*,*)'with default value ''ini_huge_h2oice'''
    123124        do ig = 1,ngrid
    124             if (watercaptag(ig)) h2o_ice(ig,:) = ini_h2o_bigreservoir
    125         enddo
     125            if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
     126        enddo
     127    else
     128        ! The variations of infinite reservoirs during the PCM years are taken into account
     129        h2o_ice = h2o_ice + watercap
    126130    endif
    127131
     
    319323    ! h2o ice
    320324    h2o_ice = 0.
    321     write(*,*)'There is no "startpem.nc" so ''h2o_ice'' is reconstructed from watercaptag with default value ''ini_h2o_bigreservoir''.'
     325    write(*,*)'There is no "startpem.nc" so ''h2o_ice'' is reconstructed from watercaptag with default value ''ini_huge_h2oice''.'
    322326        do ig = 1,ngrid
    323             if (watercaptag(ig)) h2o_ice(ig,:) = ini_h2o_bigreservoir
     327            if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
    324328        enddo
    325329
Note: See TracChangeset for help on using the changeset viewer.