Changeset 3159 for trunk/LMDZ.COMMON/libf/evolution

Timestamp:

Dec 14, 2023, 6:56:40 PM (12 months ago)

Author:

jbclement

Message:

PEM:

• A number of corrections related to r3149 for CO2 ice management:

  'co2_ice' was not transferred to 'perennial_co2ice' at the end of the PEM run;
  In the Mars PCM, 'perennial_co2ice' is now correctly handled regarding the frost in case of CO2 sublimation/condensation.
  

• 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.
• Thresholds value for ice management can now be set in the "run_PEM.def".

JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• conf_pem.F90 (modified) (6 diffs)
• constants_marspem_mod.F90 (modified) (1 diff)
• deftank/run_PEM.def (modified) (2 diffs)
• pem.F90 (modified) (7 diffs)
• pemetat0.F90 (modified) (2 diffs)
• stopping_crit_mod.F90 (modified) (6 diffs)
• time_evol_mod.F90 (modified) (1 diff)

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

@@ -172 +172 @@
 - The PEM deftank folder is now in LMDZ.COMMON/libf/evolution/ (not anymore in the Mars PCM deftank).
 - Small corrections related to r3149.
+
+== 14/12/2023 == JBC
+- A number of corrections related to r3149 for CO2 ice management:
+      > 'co2_ice' was not transferred to 'perennial_co2ice' at the end of the PEM run;
+      > In the Mars PCM, 'perennial_co2ice' is now correctly handled regarding the frost in case of CO2 sublimation/condensation.
+- 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.
+- Thresholds value for ice management can now be set in the "run_PEM.def".

trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90

@@ -23 +23 @@
 #endif
 
-use time_evol_mod,    only: year_bp_ini, dt_pem, water_ice_criterion, co2_ice_criterion, ps_criterion, Max_iter_pem, &
+use time_evol_mod,         only: year_bp_ini, dt_pem, h2o_ice_crit, co2_ice_crit, ps_criterion, Max_iter_pem, &
                           evol_orbit_pem, var_obl, var_ecc, var_lsp, convert_years
-use comsoil_h_pem,    only: soil_pem, fluxgeo, ini_h2o_bigreservoir, depth_breccia, depth_bedrock, reg_thprop_dependp
-use adsorption_mod,   only: adsorption_pem
-use glaciers_mod,     only: co2_ice_flow, h2o_ice_flow
-use ice_table_mod,    only: icetable_equilibrium, icetable_dynamic
-use time_phylmdz_mod, only: ecritphy
+use comsoil_h_pem,         only: soil_pem, fluxgeo, ini_h2o_bigreservoir, depth_breccia, depth_bedrock, reg_thprop_dependp
+use adsorption_mod,        only: adsorption_pem
+use glaciers_mod,          only: co2_ice_flow, h2o_ice_flow
+use ice_table_mod,         only: icetable_equilibrium, icetable_dynamic
+use time_phylmdz_mod,      only: ecritphy
+use constants_marspem_mod, only: inf_h2oice_threshold, metam_co2ice_threshold, metam_h2oice_threshold
 
 implicit none
@@ -52 +53 @@
 close(100)
 
-! #---------- Output parameters
+!#---------- Output parameters ----------#
 ! Frequency of outputs for the PEM
 ecritphy = 1 ! Default value: every year
 call getin('ecritpem',ecritphy)
 
-! #---------- Orbital parameters
+!#---------- Orbital parameters ----------#
 evol_orbit_pem = .false.
 call getin('evol_orbit_pem',evol_orbit_pem)
@@ -77 +78 @@
 write(*,*) 'Does Ls peri vary ?',var_lsp
 
-! #---------- Stopping criteria parameters
+!#---------- Stopping criteria parameters ----------#
 Max_iter_pem = 100000000
 call getin('Max_iter_pem',Max_iter_pem)
 
-water_ice_criterion = 0.2
-call getin('water_ice_criterion',water_ice_criterion)
+h2o_ice_crit = 0.2
+call getin('h2o_ice_crit',h2o_ice_crit)
 
-co2_ice_criterion = 0.2
-call getin('co2_ice_criterion',co2_ice_criterion)
+co2_ice_crit = 0.2
+call getin('co2_ice_crit',co2_ice_crit)
 
 ps_criterion = 0.15
@@ -93 +94 @@
 call getin('dt_pem', dt_pem)
 
-! #---------- Subsurface parameters
+!#---------- Subsurface parameters ----------#
 soil_pem = .true.
 call getin('soil_pem',soil_pem)
@@ -110 +111 @@
 write(*,*)  'Thermal properties of the regolith vary with pressure ?', reg_thprop_dependp
 
-! #---------- Layering parameters
+!#---------- Layering parameters ----------#
 fluxgeo = 0.
 call getin('fluxgeo',fluxgeo)
@@ -156 +157 @@
 endif
 
+!#---------- Ice management parameters ----------#
 ini_h2o_bigreservoir = 1.e4
 call getin('ini_h2o_bigreservoir',ini_h2o_bigreservoir)
+
+inf_h2oice_threshold = 2.e3
+call getin('inf_h2oice_threshold',inf_h2oice_threshold)
+
+metam_h2oice_threshold = 5.e-2
+call getin('metam_h2oice_threshold',metam_h2oice_threshold)
+
+metam_co2ice_threshold = 16.e3
+call getin('metam_co2ice_threshold',metam_co2ice_threshold)
 
 END SUBROUTINE conf_pem

trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90

@@ -37 +37 @@
 real, parameter :: Lco2 =  5.71e5 ! Pilorget and Forget 2016
 
-! Threshold to consider the amount of H2O ice as an infinite reservoir
-real, parameter :: inf_h2oice_threshold = 2000.  !~ 1 m
+! Thresholds for ice management
+real, save :: inf_h2oice_threshold   ! To consider the amount of H2O ice as an infinite reservoir
+real, save :: metam_h2oice_threshold ! To consider frost is becoming perennial H2O ice
+real, save :: metam_co2ice_threshold ! To consider frost is becoming perennial CO2 ice
 
 END MODULE constants_marspem_mod

trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def

@@ -28 +28 @@
 #  Max_iter_pem=100000000
 
-# Acceptance rate of sublimating water ice surface change? Default = 0.2
-# water_ice_criterion=0.2
+# Acceptance rate of sublimating H2o ice surface change? Default = 0.2
+# h2o_ice_crit=0.2
 
 # Acceptance rate of sublimating CO2 ice surface change? Default = 0.2
-# co2_ice_criterion=0.2
+# co2_ice_crit=0.2
 
 # Acceptance rate of pressure surface change? Default = 0.15
@@ -69 +69 @@
 # icetable_dynamic=.false.
 
-#---------- H2O reservoir parameters ----------#
-# Amount of h2o ice to initialize the big reservoir if the variable is not present in startfi_PEM.nc? Default = 1.e4
-# ini_h2o_bigreservoir=1e4
+#---------- Ice management parameters ----------#
+# Amount of h2o ice to initialize the big reservoir if the variable is not present in "startfi_PEM.nc"? Default = 1.e4
+# ini_h2o_bigreservoir=1.e4
+
+# Threshold to consider the amount of H2O ice as an infinite reservoir? Default = 2.e3
+# inf_h2oice_threshold=2.e3
+
+# Threshold to consider frost is becoming perennial H2O ice? Default = 5.e-2
+# metam_h2oice_threshold=5.e-2
+
+# Threshold to consider frost is becoming perennial CO2 ice? Default = 16.e3
+# metam_co2ice_threshold=16.e3
 
 # Some definitions for the physics, in file 'callphys.def'

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

@@ -42 +42 @@
 use stopping_crit_mod,           only: stopping_crit_h2o_ice, stopping_crit_co2
 use constants_marspem_mod,       only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, &
-                                       m_noco2, inf_h2oice_threshold
+                                       m_noco2, inf_h2oice_threshold, metam_co2ice_threshold, metam_h2oice_threshold
 use evol_ice_mod,                only: evol_co2_ice, evol_h2o_ice
 use comsoil_h_PEM,               only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
@@ -189 +189 @@
 real, dimension(:,:,:), allocatable :: co2_ice_PCM        ! Physics x NSLOPE x Times field: co2 ice given by the PCM [kg/m^2]
 real, dimension(:,:),   allocatable :: co2_ice_ini_sublim ! physical point field: Logical array indicating sublimating point of co2 ice
-real, dimension(:,:),   allocatable :: initial_h2o_ice    ! physical point field: Logical array indicating if there is water ice at initial state
+real, dimension(:,:),   allocatable :: initial_h2o_ice    ! physical point field: Logical array indicating if there is h2o ice at initial state
 real, dimension(:,:),   allocatable :: tend_co2_ice       ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
 real, dimension(:,:),   allocatable :: tend_co2_ice_ini   ! physical point x slope field x nslope: Tendency of evolution of perennial co2 ice over a year in the PCM
@@ -550 +550 @@
 !    I_g Save initial PCM situation
 !------------------------
-! We save the places where water ice is sublimating
-! We compute the surface of water ice sublimating
+! We save the places where h2o ice is sublimating
+! We compute the surface of h2o ice sublimating
 co2_surf_ini = 0.
 h2o_surf_ini = 0.
@@ -848 +848 @@
 ! II_d.5 Update the mass of the regolith adsorbed
         if (adsorption_pem) then
-            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tend_h2o_ice,tend_co2_ice,       &
+            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tend_h2o_ice,tend_co2_ice,                   &
                                      h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
                                      h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
@@ -858 +858 @@
                     do l = 1,nsoilmx_PEM - 1
                         totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l + 1) - layer_PEM(l))* &
-                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
+                                               subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
                         totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l + 1) - layer_PEM(l))* &
-                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
+                                               subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
                     enddo
                 enddo
@@ -907 +907 @@
         select case (stopPEM)
             case(1)
-                write(*,*) "STOPPING because surface of water ice sublimating is too low:", stopPEM, "See message above."
+                write(*,*) "STOPPING because surface of h2o ice sublimating is too low:", stopPEM, "See message above."
             case(2)
-                write(*,*) "STOPPING because tendencies on water ice = 0:", stopPEM, "See message above."
+                write(*,*) "STOPPING because tendencies on h2o ice = 0:", stopPEM, "See message above."
             case(3)
                 write(*,*) "STOPPING because surface of co2 ice sublimating is too low:", stopPEM, "See message above."
@@ -940 +940 @@
 ! III_a.1 Ice update (for startfi)
 
-! H2O ice
 watercap = 0.
+perennial_co2ice = co2_ice
 do ig = 1,ngrid
+    ! H2O ice metamorphism
+    if (sum(qsurf(ig,igcm_h2o_ice,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_h2oice_threshold) then
+        h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
+        qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
+    endif
+
+    ! Is H2O ice still considered as an infinite reservoir for the PCM?
     if (sum(h2o_ice(ig,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then
-    ! If there is enough ice to be considered as an infinite reservoir
+        ! There is enough ice to be considered as an infinite reservoir
         watercaptag(ig) = .true.
     else
-    ! If there too little ice to be considered as an infinite reservoir
-    ! Ice is transferred to the frost
+        ! There too little ice to be considered as an infinite reservoir so ice is transferred to the frost
         watercaptag(ig) = .false.
         qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
         h2o_ice(ig,:) = 0.
+    endif
+
+    ! CO2 ice metamorphism
+    if (sum(qsurf(ig,igcm_co2,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_co2ice_threshold) then
+        perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
+        qsurf(ig,igcm_co2,:) = metam_co2ice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
     endif
 enddo

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

@@ -82 +82 @@
 #ifdef CPP_STD
     logical, dimension(ngrid) :: watercaptag
-    watercaptag(:) = .false.
+    watercaptag = .false.
 #endif
 
@@ -152 +152 @@
                             TI_PEM(ig,iloop,islope) = TI_breccia
                         enddo
-                    else ! we keep the high ti values
+                    else ! we keep the high TI values
                         do iloop = index_breccia + 1,index_bedrock
                             TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)

trunk/LMDZ.COMMON/libf/evolution/stopping_crit_mod.F90

@@ -16 +16 @@
 SUBROUTINE stopping_crit_h2o_ice(cell_area,surf_ini,h2o_ice,stopPEM,ngrid)
 
-use time_evol_mod, only: water_ice_criterion
+use time_evol_mod, only: h2o_ice_crit
 use comslope_mod,  only: subslope_dist, nslope
 
@@ -32 +32 @@
 integer,                       intent(in) :: ngrid     ! # of physical grid points
 real, dimension(ngrid),        intent(in) :: cell_area ! Area of the cells
-real, dimension(ngrid,nslope), intent(in) :: h2o_ice   ! Actual density of water ice
+real, dimension(ngrid,nslope), intent(in) :: h2o_ice   ! Actual density of h2o ice
 real,                          intent(in) :: surf_ini  ! Initial surface of h2o ice that was sublimating
 !   OUTPUT
@@ -51 +51 @@
 
 ! Check of the criterion
-if (surf_now < surf_ini*(1. - water_ice_criterion)) then
+if (surf_now < surf_ini*(1. - h2o_ice_crit)) then
     stopPEM = 1
-    write(*,*) "Reason of stopping: the surface of water ice sublimating reaches the threshold"
-    write(*,*) "surf_now < surf_ini*(1. - water_ice_criterion)", surf_now < surf_ini*(1. - water_ice_criterion)
-    write(*,*) "Current surface of water ice sublimating =", surf_now
-    write(*,*) "Initial surface of water ice sublimating =", surf_ini
-    write(*,*) "Percentage of change accepted =", water_ice_criterion*100
-else if (surf_now > surf_ini*(1. + water_ice_criterion)) then
+    write(*,*) "Reason of stopping: the surface of h2o ice sublimating reaches the threshold"
+    write(*,*) "surf_now < surf_ini*(1. - h2o_ice_crit)", surf_now < surf_ini*(1. - h2o_ice_crit)
+    write(*,*) "Current surface of h2o ice sublimating =", surf_now
+    write(*,*) "Initial surface of h2o ice sublimating =", surf_ini
+    write(*,*) "Percentage of change accepted =", h2o_ice_crit*100
+else if (surf_now > surf_ini*(1. + h2o_ice_crit)) then
     stopPEM = 1
-    write(*,*) "Reason of stopping: the surface of water ice sublimating reaches the threshold"
-    write(*,*) "surf_now > surf_ini*(1. + water_ice_criterion)", surf_now > surf_ini*(1. + water_ice_criterion)
-    write(*,*) "Current surface of water ice sublimating =", surf_now
-    write(*,*) "Initial surface of water ice sublimating =", surf_ini
-    write(*,*) "Percentage of change accepted =", water_ice_criterion*100
+    write(*,*) "Reason of stopping: the surface of h2o ice sublimating reaches the threshold"
+    write(*,*) "surf_now > surf_ini*(1. + h2o_ice_crit)", surf_now > surf_ini*(1. + h2o_ice_crit)
+    write(*,*) "Current surface of h2o ice sublimating =", surf_now
+    write(*,*) "Initial surface of h2o ice sublimating =", surf_ini
+    write(*,*) "Percentage of change accepted =", h2o_ice_crit*100
 endif
 
@@ -75 +75 @@
 SUBROUTINE stopping_crit_co2(cell_area,surf_ini,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)
 
-use time_evol_mod, only: co2_ice_criterion, ps_criterion
+use time_evol_mod, only: co2_ice_crit, ps_criterion
 use comslope_mod,  only: subslope_dist
 
@@ -91 +91 @@
 integer,                       intent(in) :: ngrid, nslope        ! # of grid physical grid points
 real, dimension(ngrid),        intent(in) :: cell_area            ! Area of the cells
-real, dimension(ngrid,nslope), intent(in) :: co2_ice              ! Actual density of water ice
+real, dimension(ngrid,nslope), intent(in) :: co2_ice              ! Actual density of h2o ice
 real,                          intent(in) :: surf_ini             ! Initial surface of co2 ice that was sublimating
 real,                          intent(in) :: global_avg_press_PCM ! Planet average pressure from the PCM start files
@@ -113 +113 @@
 
 ! Check of the criterion
-if (surf_now < surf_ini*(1. - co2_ice_criterion)) then
+if (surf_now < surf_ini*(1. - co2_ice_crit)) then
     stopPEM = 3
     write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
-    write(*,*) "surf_now < surf_ini*(1. - co2_ice_criterion)", surf_now < surf_ini*(1. - co2_ice_criterion)
+    write(*,*) "surf_now < surf_ini*(1. - co2_ice_crit)", surf_now < surf_ini*(1. - co2_ice_crit)
     write(*,*) "Current surface of co2 ice sublimating =", surf_now
     write(*,*) "Initial surface of co2 ice sublimating =", surf_ini
-    write(*,*) "Percentage of change accepted =", co2_ice_criterion*100.
-else if (surf_now > surf_ini*(1. + co2_ice_criterion)) then
+    write(*,*) "Percentage of change accepted =", co2_ice_crit*100.
+else if (surf_now > surf_ini*(1. + co2_ice_crit)) then
     stopPEM = 3
     write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
-    write(*,*) "surf_now > surf_ini*(1. + co2_ice_criterion)", surf_now > surf_ini*(1. + co2_ice_criterion)
+    write(*,*) "surf_now > surf_ini*(1. + co2_ice_crit)", surf_now > surf_ini*(1. + co2_ice_crit)
     write(*,*) "Current surface of co2 ice sublimating =", surf_now
     write(*,*) "Initial surface of co2 ice sublimating =", surf_ini
-    write(*,*) "Percentage of change accepted =", co2_ice_criterion*100.
+    write(*,*) "Percentage of change accepted =", co2_ice_crit*100.
 endif
 

trunk/LMDZ.COMMON/libf/evolution/time_evol_mod.F90

@@ -3 +3 @@
 implicit none
 
-integer :: year_bp_ini         ! Initial year (in Planetary years) of the simulation of the PEM defined in run.def (in Earth years)
-integer :: dt_pem              ! Time step used by the PEM in Planetary years
-real    :: convert_years       ! Conversion ratio from Planetary years to Earth years
-real    :: water_ice_criterion ! Percentage of change of the surface of water ice sublimating before stopping the PEM
-real    :: co2_ice_criterion   ! Percentage of change of the surface of co2 ice sublimating before stopping the PEM
-real    :: ps_criterion        ! Percentage of change of averaged surface pressure before stopping the PEM
-integer :: Max_iter_pem        ! Maximal number of iteration when converging to a steady state, read in evol.def
-logical :: evol_orbit_pem      ! True if we want to follow the orbital parameters of obl_ecc_lsp.asc, read in evol.def
-logical :: var_obl             ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for obliquity
-logical :: var_ecc             ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for eccenticity
-logical :: var_lsp             ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for ls perihelie
+integer :: year_bp_ini    ! Initial year (in Planetary years) of the simulation of the PEM defined in run.def (in Earth years)
+integer :: dt_pem         ! Time step used by the PEM in Planetary years
+real    :: convert_years  ! Conversion ratio from Planetary years to Earth years
+real    :: h2o_ice_crit   ! Percentage of change of the surface of h2o ice sublimating before stopping the PEM
+real    :: co2_ice_crit   ! Percentage of change of the surface of co2 ice sublimating before stopping the PEM
+real    :: ps_criterion   ! Percentage of change of averaged surface pressure before stopping the PEM
+integer :: Max_iter_pem   ! Maximal number of iteration when converging to a steady state, read in evol.def
+logical :: evol_orbit_pem ! True if we want to follow the orbital parameters of obl_ecc_lsp.asc, read in evol.def
+logical :: var_obl        ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for obliquity
+logical :: var_ecc        ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for eccenticity
+logical :: var_lsp        ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for ls perihelie
 
 END MODULE time_evol_mod