Changeset 3525

Timestamp:

Nov 19, 2024, 5:44:27 PM (4 days ago)

Author:

jbclement

Message:

PEM:
Computation of <soil thermal intertia> and <H2O mass subsurface/surface exchange> according to the presence of subsurface ice provided by the (Norbert's) dynamic method + few cleanings.
JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• compute_soiltemp_mod.F90 (modified) (3 diffs)
• ice_table_mod.F90 (modified) (14 diffs)
• pem.F90 (modified) (5 diffs)
• pemetat0.F90 (modified) (6 diffs)
• soil_thermalproperties_mod.F90 (modified) (8 diffs)

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

@@ -489 +489 @@
 == 14/11/2024 == JBC
 Committing the right correction in previous commit (r3516) would have been better...
+
+== 19/11/2024 == JBC
+Computation of <soil thermal intertia> and <H2O mass subsurface/surface exchange> according to the presence of subsurface ice provided by the (Norbert's) dynamic method + few cleanings.

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

@@ -90 +90 @@
 endif ! of if (firstcall)
 
-IF (.not. firstcall) THEN
+if (.not. firstcall) THEN
 ! 2. Compute soil temperatures
 ! First layer:
@@ -120 +120 @@
 END SUBROUTINE compute_tsoil_pem
 
-
+!=======================================================================
 
 SUBROUTINE ini_tsoil_pem(ngrid,nsoil,therm_i,tsurf,tsoil)
@@ -221 +221 @@
 END SUBROUTINE ini_tsoil_pem
 
-
-
 END MODULE compute_soiltemp_mod

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

@@ -18 +18 @@
 !-----------------------------------------------------------------------
 contains
+
 !-----------------------------------------------------------------------
 SUBROUTINE ini_ice_table(ngrid,nslope,nsoil)
@@ -45 +46 @@
 
 !------------------------------------------------------------------------------------------------------
-SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness)
+SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
@@ -51 +52 @@
 !!! density at the surface and at depth.
 !!! Computations are made following the methods in Schorgofer et al., 2005
-!!! This SUBROUTINE only gives the ice table at equilibrium and does not consider exchange with the atmosphere
+!!! This subroutine only gives the ice table at equilibrium and does not consider exchange with the atmosphere
 !!!
 !!! Author: LL
@@ -64 +65 @@
 ! Inputs
 ! ------
-integer,                                 intent(in) :: ngrid, nslope, nsoil_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM
-logical, dimension(ngrid),               intent(in) :: watercaptag              ! Boolean to check the presence of a perennial glacier
-real, dimension(ngrid,nslope),           intent(in) :: rhowatersurf_ave         ! Water density at the surface, yearly averaged [kg/m^3]
-real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: rhowatersoil_ave         ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged  [kg/m^3]
-real, dimension(ngrid,nslope),           intent(in) :: regolith_inertia         ! Thermal inertia of the regolith layer [SI]
+integer,                             intent(in) :: ngrid, nslope, nsoil ! Size of the physical grid, number of subslope, number of soil layer in the PEM
+logical, dimension(ngrid),           intent(in) :: watercaptag          ! Boolean to check the presence of a perennial glacier
+real, dimension(ngrid,nslope),       intent(in) :: rhowatersurf_ave     ! Water density at the surface, yearly averaged [kg/m^3]
+real, dimension(ngrid,nsoil,nslope), intent(in) :: rhowatersoil_ave     ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged [kg/m^3]
+real, dimension(ngrid,nslope),       intent(in) :: regolith_inertia     ! Thermal inertia of the regolith layer [SI]
 ! Ouputs
 ! ------
@@ -76 +77 @@
 ! ------
 integer                       :: ig, islope, isoil, isoilend ! loop variables
-real, dimension(nsoil_PEM)    :: diff_rho                    ! difference of water vapor density between the surface and at depth [kg/m^3]
+real, dimension(nsoil)        :: diff_rho                    ! difference of water vapor density between the surface and at depth [kg/m^3]
 real                          :: ice_table_end               ! depth of the end of the ice table  [m]
 real, dimension(ngrid,nslope) :: previous_icetable_depth     ! Ice table computed at previous ice depth [m]
@@ -87 +88 @@
     if (watercaptag(ig)) then
         ice_table_beg(ig,:) = 0.
-        ice_table_thickness(ig,:) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
+        ice_table_thickness(ig,:) = layer_PEM(nsoil) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
     else
         do islope = 1,nslope
             ice_table_beg(ig,islope) = -1.
             ice_table_thickness(ig,islope) = 0.
-            do isoil = 1,nsoil_PEM
+            do isoil = 1,nsoil
                 diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope)
             enddo
             if (diff_rho(1) > 0) then ! ice is at the surface
                 ice_table_beg(ig,islope) = 0.
-                do isoilend = 2,nsoil_PEM - 1
+                do isoilend = 2,nsoil - 1
                     if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend + 1) < 0.) then
                         call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend),mlayer_PEM(isoilend + 1),ice_table_end)
@@ -105 +106 @@
                 enddo
             else
-                do isoil = 1,nsoil_PEM - 1 ! general case, we find the ice table depth by doing a linear approximation between the two depth, and then solve the first degree equation to find the root
+                do isoil = 1,nsoil - 1 ! general case, we find the ice table depth by doing a linear approximation between the two depth, and then solve the first degree equation to find the root
                     if (diff_rho(isoil) < 0 .and. diff_rho(isoil + 1) > 0.) then
                         call findroot(diff_rho(isoil),diff_rho(isoil + 1),mlayer_PEM(isoil),mlayer_PEM(isoil + 1),ice_table_beg(ig,islope))
                         ! Now let's find the end of the ice table
-                        ice_table_thickness(ig,islope) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
-                        do isoilend = isoil + 1,nsoil_PEM - 1
+                        ice_table_thickness(ig,islope) = layer_PEM(nsoil) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
+                        do isoilend = isoil + 1,nsoil - 1
                             if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend + 1) < 0.) then
                                 call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend),mlayer_PEM(isoilend + 1),ice_table_end)
@@ -140 +141 @@
 
 !-----------------------------------------------------------------------
-SUBROUTINE compute_massh2o_exchange_ssi(ngrid,nslope,nsoil_PEM,former_ice_table_thickness,new_ice_table_thickness,ice_table_depth,tsurf,tsoil,delta_m_h2o)
+SUBROUTINE compute_massh2o_exchange_ssi(ngrid,nslope,nsoil,icetable_thickness_old,ice_porefilling_old,tsurf,tsoil,delta_m_h2o)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
@@ -161 +162 @@
 ! Inputs
 ! ------
-integer,                                 intent(in) :: ngrid, nslope, nsoil_PEM   ! Size of the physical grid, number of subslope, number of soil layer in the PEM
-real, dimension(ngrid,nslope),           intent(in) :: former_ice_table_thickness ! ice table thickness at the former iteration [m]
-real, dimension(ngrid,nslope),           intent(in) :: new_ice_table_thickness    ! ice table thickness at the current iteration [m]
-real, dimension(ngrid,nslope),           intent(in) :: ice_table_depth            ! ice table depth [m]
-real, dimension(ngrid,nslope),           intent(in) :: tsurf                      ! Surface temperature [K]
-real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil                      ! Soil temperature [K]
+integer,                             intent(in) :: ngrid, nslope, nsoil   ! Size of the physical grid, number of subslope, number of soil layer in the PEM
+real, dimension(ngrid,nslope),       intent(in) :: icetable_thickness_old ! Ice table thickness at the previous iteration [m]
+real, dimension(ngrid,nsoil,nslope), intent(in) :: ice_porefilling_old    ! Ice pore filling at the previous iteration [m]
+real, dimension(ngrid,nslope),       intent(in) :: tsurf                  ! Surface temperature [K]
+real, dimension(ngrid,nsoil,nslope), intent(in) :: tsoil                  ! Soil temperature [K]
 ! Outputs
 ! -------
@@ -172 +172 @@
 ! Locals
 !-------
-integer                       :: ig, islope, ilay, iref ! loop index
-real, dimension(ngrid,nslope) :: rho                    ! density of water ice [kg/m^3]
-real, dimension(ngrid,nslope) :: Tice                   ! ice temperature [k]
+integer :: ig, islope, isoil ! loop index
+real    :: rho               ! density of water ice [kg/m^3]
+real    :: Tice              ! ice temperature [k]
 ! Code
 ! ----
-rho = 0.
-Tice = 0.
-!1. First let's compute Tice using a linear interpolation between the mlayer level
-do ig = 1,ngrid
-    do islope = 1,nslope
-        call compute_Tice_pem(nsoil_PEM,tsoil(ig,:,islope),tsurf(ig,islope),ice_table_depth(ig,islope),Tice(ig,islope))
-        rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2 + 0.0351* Tice(ig,islope) + 933.5030 ! Rottgers, 2012
-    enddo
-enddo
-
-!2. Let's compute the amount of ice that has sublimated in each subslope
-do ig = 1,ngrid
-    do islope = 1,nslope
-        delta_m_h2o(ig) = delta_m_h2o(ig) + porosity*rho(ig,islope)*(new_ice_table_thickness(ig,islope) - former_ice_table_thickness(ig,islope)) & ! convention > 0. <=> it condenses
-                          *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
-    enddo
-enddo
+
+! Let's compute the amount of ice that has sublimated in each subslope
+if (icetable_equilibrium) then
+    delta_m_h2o = 0.
+    do ig = 1,ngrid
+        do islope = 1,nslope
+            call compute_Tice_pem(nsoil,tsoil(ig,:,islope),tsurf(ig,islope),icetable_depth(ig,islope),Tice)
+            rho = -3.5353e-4*Tice**2 + 0.0351*Tice + 933.5030 ! Rottgers, 2012
+            delta_m_h2o(ig) = delta_m_h2o(ig) + porosity*rho*(icetable_thickness(ig,islope) - icetable_thickness_old(ig,islope)) & ! convention > 0. <=> it condenses
+                              *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
+        enddo
+    enddo
+else if (icetable_dynamic) then
+    delta_m_h2o = 0.
+    do ig = 1,ngrid
+        do islope = 1,nslope
+            do isoil = 1,nsoil
+                call compute_Tice_pem(nsoil,tsoil(ig,:,islope),tsurf(ig,islope),mlayer_PEM(isoil - 1),Tice)
+                rho = -3.5353e-4*Tice**2 + 0.0351*Tice + 933.5030 ! Rottgers, 2012
+                delta_m_h2o(ig) = delta_m_h2o(ig) + rho*(ice_porefilling(ig,isoil,islope) - ice_porefilling_old(ig,isoil,islope)) & ! convention > 0. <=> it condenses
+                                  *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
+            enddo
+        enddo
+    enddo
+endif
 
 END SUBROUTINE compute_massh2o_exchange_ssi
@@ -285 +293 @@
 
 ! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010))
-! 2.0 preliminary: depth to shallowest  ice (discontinuity at interface)
+! 2.0 preliminary: depth to shallowest ice (discontinuity at interface)
 index_firstice = -1
 do ilay = 1,nsoilmx_PEM
@@ -294 +302 @@
 enddo
 
-! 2.1: now we can computeCompute d_z (eta* d_z(rho))
+! 2.1: now we can compute
 call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM - 1),dz_eta)
 
@@ -368 +376 @@
 
 !-----------------------------------------------------------------------
-SUBROUTINE compute_Tice_pem(nsoil, ptsoil, ptsurf, ice_depth, Tice)
+SUBROUTINE compute_Tice_pem(nsoil,ptsoil,ptsurf,ice_depth,Tice)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
@@ -383 +391 @@
 ! Inputs
 ! ------
-integer,                 intent(in) :: nsoil     ! Number of soil layers
-real, dimension(nsoil),  intent(in) :: ptsoil    ! Soil temperature (K)
-real,                    intent(in) :: ptsurf    ! Soil temperature (K)
-real,                    intent(in) :: ice_depth ! Ice depth (m)
+integer,                intent(in) :: nsoil     ! Number of soil layers
+real, dimension(nsoil), intent(in) :: ptsoil    ! Soil temperature (K)
+real,                   intent(in) :: ptsurf    ! Surface temperature (K)
+real,                   intent(in) :: ice_depth ! Ice depth (m)
 
 ! Outputs
-! ------
+! -------
 real, intent(out) :: Tice ! Ice temperatures (K)
 

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

@@ -223 +223 @@
 logical                               :: bool_sublim                      ! logical to check if there is sublimation or not
 logical, dimension(:,:),  allocatable :: co2ice_disappeared               ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
-real, dimension(:,:),     allocatable :: icetable_thickness_prev_iter     ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
+real, dimension(:,:),     allocatable :: icetable_thickness_old           ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
+real, dimension(:,:,:),   allocatable :: ice_porefilling_old              ! ngrid x nslope: Ice pore filling at the previous iteration [m]
 real, dimension(:),       allocatable :: delta_h2o_icetablesublim         ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
 real, dimension(:),       allocatable :: porefill                         ! Pore filling (output) to compute the dynamic ice table
@@ -555 +556 @@
 allocate(delta_co2_adsorbded(ngrid))
 allocate(co2ice_disappeared(ngrid,nslope))
-allocate(icetable_thickness_prev_iter(ngrid,nslope))
+allocate(icetable_thickness_old(ngrid,nslope))
+allocate(ice_porefilling_old(ngrid,nsoilmx_PEM,nslope))
 allocate(delta_h2o_icetablesublim(ngrid))
 allocate(delta_h2o_adsorbded(ngrid))
@@ -927 +929 @@
         if (icetable_equilibrium) then
             write(*,*) "Compute ice table (equilibrium method)"
-            icetable_thickness_prev_iter = icetable_thickness
+            icetable_thickness_old = icetable_thickness
             call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
-            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_prev_iter,icetable_thickness,icetable_depth,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
+            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
         else if (icetable_dynamic) then
             write(*,*) "Compute ice table (dynamic method)"
+            ice_porefilling_old = ice_porefilling
             allocate(porefill(nsoilmx_PEM))
             do ig = 1,ngrid
@@ -941 +944 @@
             enddo
             deallocate(porefill)
-            !call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_prev_iter,icetable_thickness,icetable_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
+            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
         endif
 
 ! II_d.4 Update the soil thermal properties
-        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,d_h2oice,h2o_ice,global_avg_press_new,icetable_depth,icetable_thickness,TI_PEM)
+        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,d_h2oice,h2o_ice,global_avg_press_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
 
 ! II_d.5 Update the mass of the regolith adsorbed
@@ -1230 +1233 @@
 deallocate(watersurf_density_avg,watersoil_density_timeseries,Tsurfavg_before_saved)
 deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg)
-deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,icetable_thickness_prev_iter)
+deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim)
+deallocate(icetable_thickness_old,ice_porefilling_old)
 deallocate(is_co2ice_ini,co2ice_disappeared,ini_co2ice_sublim,ini_h2oice_sublim,stratif)
 !----------------------------- END OUTPUT ------------------------------

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

@@ -7 +7 @@
 !=======================================================================
 
-SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table_depth,ice_table_thickness,      &
-                    ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,d_h2oice,d_co2ice,co2_ice,h2o_ice, &
-                    global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,                         &
+SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table_depth,ice_table_thickness, &
+                    ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,d_h2oice,d_co2ice,co2_ice,h2o_ice,    &
+                    global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,                    &
                     m_h2o_regolith_phys,deltam_h2o_regolith_phys,stratif)
 
@@ -16 +16 @@
 use comsoil_h,                  only: volcapa, inertiedat
 use adsorption_mod,             only: regolith_adsorption, adsorption_pem
-use ice_table_mod,              only: computeice_table_equilibrium, icetable_equilibrium, icetable_dynamic
+use ice_table_mod,              only: computeice_table_equilibrium, icetable_depth, icetable_thickness, icetable_equilibrium, icetable_dynamic
 use constants_marspem_mod,      only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock
 use soil_thermalproperties_mod, only: update_soil_thermalproperties
@@ -315 +315 @@
                 write(*,*)'will reconstruct the values of the ice table given the current state'
                 call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
-                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
+                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
                 do islope = 1,nslope
                     call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
@@ -327 +327 @@
                 write(*,*)'will reconstruct the values of the ice table given the current state'
                 ice_table_depth = -9999.
-                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
+                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
                 do islope = 1,nslope
                     call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
@@ -504 +504 @@
         if (icetable_equilibrium) then
             call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
-            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
+            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
             do islope = 1,nslope
                 call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
@@ -512 +512 @@
             ice_porefilling = 1.
             ice_table_depth = -9999.
-            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
+            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
             do islope = 1,nslope
                 call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))

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

@@ -9 +9 @@
 !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+! Constants:
+real, parameter :: reg_inertie_thresold = 370. ! Above this thermal inertia, the regolith has too much cementation to be dependant on the pressure [J/m^2/K/s^1/2]
+real, parameter :: reg_inertie_minvalue = 50.  ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
+real, parameter :: reg_inertie_maxvalue = 370. ! Maximum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
+real, parameter :: P610 = 610.0                ! current average pressure on Mars [Pa]
+real, parameter :: C = 0.0015                  ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [unitless]
+real, parameter :: K = 8.1*1e4                 ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [Torr, or 133.3Pa]
+real, parameter :: Pa2Torr = 1./133.3          ! Conversion from Pa to tor [Pa/Torr]
 
 !=======================================================================
@@ -50 +59 @@
     ice_thermalinertia = inertie_purewaterice
 else
-    ice_thermalinertia = sqrt(surf_thermalinertia**2 + porosity*pore_filling*inertie_purewaterice**2)  ! Siegler et al., 2012
+    ice_thermalinertia = sqrt(surf_thermalinertia**2 + porosity*pore_filling*inertie_purewaterice**2) ! Siegler et al., 2012
 endif
 
@@ -56 +65 @@
 !=======================================================================
 
-SUBROUTINE update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tendencies_waterice,waterice,p_avg_new,ice_depth,ice_thickness,TI_PEM)
+SUBROUTINE update_soil_thermalproperties(ngrid,nslope,nsoil,tendencies_waterice,waterice,p_avg_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
 
 use comsoil_h,             only: volcapa
@@ -65 +74 @@
 
 ! Input:
-integer,                       intent(in) :: ngrid, nslope, nsoil_PEM ! Shape of the arrays: physical grid, number of sub-grid slopes, number of layer in the soil
-real,                          intent(in) :: p_avg_new                ! Global average surface pressure [Pa]
-real, dimension(ngrid,nslope), intent(in) :: tendencies_waterice      ! Tendencies on the water ice [kg/m^2/year]
-real, dimension(ngrid,nslope), intent(in) :: waterice                 ! Surface Water ice [kg/m^2]
-real, dimension(ngrid,nslope), intent(in) :: ice_depth                ! Ice table depth [m]
-real, dimension(ngrid,nslope), intent(in) :: ice_thickness            ! Ice table thickness [m]
+integer,                             intent(in) :: ngrid, nslope, nsoil ! Shape of the arrays: physical grid, number of sub-grid slopes, number of layer in the soil
+real,                                intent(in) :: p_avg_new            ! Global average surface pressure [Pa]
+real, dimension(ngrid,nslope),       intent(in) :: tendencies_waterice  ! Tendencies on the water ice [kg/m^2/year]
+real, dimension(ngrid,nslope),       intent(in) :: waterice             ! Surface Water ice [kg/m^2]
+real, dimension(ngrid,nslope),       intent(in) :: icetable_depth       ! Ice table depth [m]
+real, dimension(ngrid,nslope),       intent(in) :: icetable_thickness   ! Ice table thickness [m]
+real, dimension(ngrid,nsoil,nslope), intent(in) :: ice_porefilling      ! Ice porefilling [m^3/m^3]
+logical,                             intent(in) :: icetable_equilibrium, icetable_dynamic ! Computing method for ice table
 
 ! Outputs:
-real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! Soil Thermal Inertia [J/m^2/K/s^1/2]
-
-! Constants:
-real :: reg_inertie_thresold = 370. ! Above this thermal inertia, the regolith has too much cementation to be dependant on the pressure [J/m^2/K/s^1/2]
-real :: reg_inertie_minvalue = 50.  ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
-real :: reg_inertie_maxvalue = 370. ! Maximum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
-real :: ice_inertia                 ! Inertia of water ice [SI]
-real :: P610 = 610.0                ! current average pressure on Mars [Pa]
-real :: C = 0.0015                  ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [uniteless]
-real :: K = 8.1*1e4                 ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [tor, or 133.3Pa]
-real :: Pa2Tor = 1./133.3           ! Conversion from Pa to tor [Pa/tor]
+real, dimension(ngrid,nsoil,nslope), intent(inout) :: TI_PEM ! Soil Thermal Inertia [J/m^2/K/s^1/2]
 
 ! Local variables:
@@ -91 +92 @@
 real    :: ice_bottom_depth                       ! Bottom depth of the subsurface ice [m]
 real    :: d_part                                 ! Regolith particle size [micrometer]
+real    :: ice_inertia                            ! Inertia of water ice [SI]
 
 write(*,*) "Update soil properties"
@@ -101 +103 @@
         if (reg_thprop_dependp) then
             if (TI_PEM(ig,1,islope) < reg_inertie_thresold) then
-                d_part = (regolith_inertia(ig,islope)**2/(volcapa*C*(P610*Pa2Tor)**(0.6)))**(-1./(0.11*log10(P610*Pa2Tor/K))) ! compute particle size, in micrometer
-                regolith_inertia(ig,islope) = sqrt(volcapa*C*(p_avg_new*Pa2Tor)**(0.6)*d_part**(-0.11*log10(p_avg_new*Pa2Tor/K)))
+                d_part = (regolith_inertia(ig,islope)**2/(volcapa*C*(P610*Pa2Torr)**(0.6)))**(-1./(0.11*log10(P610*Pa2Torr/K))) ! compute particle size, in micrometer
+                regolith_inertia(ig,islope) = sqrt(volcapa*C*(p_avg_new*Pa2Torr)**(0.6)*d_part**(-0.11*log10(p_avg_new*Pa2Torr/K)))
                 if (regolith_inertia(ig,islope) > reg_inertie_maxvalue) regolith_inertia(ig,islope) = reg_inertie_maxvalue
                 if (regolith_inertia(ig,islope) < reg_inertie_minvalue) regolith_inertia(ig,islope) = reg_inertie_minvalue
@@ -135 +137 @@
                                               (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
                                               ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
-        do isoil = index_bedrock + 2,nsoil_PEM
+        do isoil = index_bedrock + 2,nsoil
             TI_PEM(ig,isoil,islope) = TI_bedrock
         enddo
@@ -141 +143 @@
 enddo ! islope
 
-!  3. Build new TI in case of ice table
+! 3. Build new TI in case of ice table
 do ig = 1,ngrid
     do islope = 1,nslope
-        if (ice_depth(ig,islope) > -1.e-6) then
+        if (icetable_depth(ig,islope) > -1.e-6) then
         ! 3.0 Case where it is perennial ice
-            if (ice_depth(ig,islope) < 1.e-10) then
+            if (icetable_depth(ig,islope) < 1.e-10) then
                 call ice_thermal_properties(.true.,1.,0.,ice_inertia)
-                do isoil = 1,nsoil_PEM
+                do isoil = 1,nsoil
                     TI_PEM(ig,isoil,islope) = ice_inertia
                 enddo
             else
-                call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),ice_inertia)
-        ! 3.1.1 find the index of the mixed layer
-                iref = 0 ! initialize iref
-                do isoil = 1,nsoil_PEM ! loop on layers to find the beginning of the ice table
-                    if (ice_depth(ig,islope) >= layer_PEM(isoil)) then
-                        iref = isoil ! pure regolith layer up to here
-                    else
-                        exit ! correct iref was obtained in previous cycle
-                    endif
-                enddo
-        ! 3.1.2 find the index of the end of the ice table
-                iend = 0 ! initialize iend
-                ice_bottom_depth = ice_depth(ig,islope) + ice_thickness(ig,islope)
-                do isoil = 1,nsoil_PEM ! loop on layers to find the end of the ice table
-                    if (ice_bottom_depth >= layer_PEM(isoil)) then
-                        iend = isoil ! pure regolith layer up to here
-                    else
-                        exit ! correct iref was obtained in previous cycle
-                    endif
-                enddo
-        ! 3.2 Build the new ti
-                if (iref < nsoil_PEM) then
-                    if (iref == iend) then
-                        ! Ice table begins and end in the same mixture Mixtures with three components
-                        if (iref /= 0) then ! mixed layer
-                            TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/                             &
-                                                         (((ice_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + &
-                                                         ((ice_bottom_depth - ice_depth(ig,islope))/(ice_inertia**2)) +            &
-                                                         ((layer_PEM(iref + 1) - ice_bottom_depth)/(TI_PEM(ig,iref + 1,islope)**2))))
-                        else ! first layer is already a mixed layer
-                            ! (ie: take layer(iref=0)=0)
-                            TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/ &
-                                                  (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) +           &
-                                                  ((ice_bottom_depth - ice_depth(ig,islope))/(ice_inertia**2)) + &
-                                                  ((layer_PEM(2) - ice_bottom_depth)/(TI_PEM(ig,2,islope)**2))))
-                        endif ! of if (iref /= 0)
-                    else
-                        if (iref /= 0) then ! mixed layer
-                            TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/                             &
-                                                         (((ice_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + &
-                                                         ((layer_PEM(iref + 1) - ice_depth(ig,islope))/(ice_inertia**2))))
-                        else ! first layer is already a mixed layer
-                            ! (ie: take layer(iref=0)=0)
-                            TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/                                 &
-                                                  (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) + &
-                                                  ((layer_PEM(1) - ice_depth(ig,islope))/(ice_inertia**2))))
-                        endif ! of if (iref /= 0)
-                    endif ! iref == iend
-        ! 3.3 Build the new ti
-                    do isoil = iref + 2,iend
-                        TI_PEM(ig,isoil,islope) = ice_inertia
+                if (icetable_equilibrium) then
+                    call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),ice_inertia)
+                    ! 3.1.1 find the index of the mixed layer
+                    iref = 0 ! initialize iref
+                    do isoil = 1,nsoil ! loop on layers to find the beginning of the ice table
+                        if (icetable_depth(ig,islope) >= layer_PEM(isoil)) then
+                            iref = isoil ! pure regolith layer up to here
+                        else
+                            exit ! correct iref was obtained in previous cycle
+                        endif
                     enddo
-                    if (iend < nsoil_PEM) then
-                        TI_PEM(ig,iend + 1,islope) = sqrt((layer_PEM(iend + 1) - layer_PEM(iend))/                         &
-                                                     (((ice_bottom_depth - layer_PEM(iend))/(TI_PEM(ig,iend,islope)**2)) + &
-                                                     ((layer_PEM(iend + 1) - ice_bottom_depth)/(ice_inertia**2))))
-                    endif
-                endif ! of if (iref < nsoilmx)
+                    ! 3.1.2 find the index of the end of the ice table
+                    iend = 0 ! initialize iend
+                    ice_bottom_depth = icetable_depth(ig,islope) + icetable_thickness(ig,islope)
+                    do isoil = 1,nsoil ! loop on layers to find the end of the ice table
+                        if (ice_bottom_depth >= layer_PEM(isoil)) then
+                            iend = isoil ! pure regolith layer up to here
+                        else
+                            exit ! correct iref was obtained in previous cycle
+                        endif
+                    enddo
+                    ! 3.2 Build the new ti
+                    if (iref < nsoil) then
+                        if (iref == iend) then
+                            ! Ice table begins and end in the same mixture with three components
+                            if (iref /= 0) then ! mixed layer
+                                TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/                                  &
+                                                             (((icetable_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + &
+                                                             ((ice_bottom_depth - icetable_depth(ig,islope))/(ice_inertia**2)) +            &
+                                                             ((layer_PEM(iref + 1) - ice_bottom_depth)/(TI_PEM(ig,iref + 1,islope)**2))))
+                            else ! first layer is already a mixed layer
+                                ! (ie: take layer(iref=0)=0)
+                                TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/                                                &
+                                                      (((icetable_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) +           &
+                                                      ((ice_bottom_depth - icetable_depth(ig,islope))/(ice_inertia**2)) + &
+                                                      ((layer_PEM(2) - ice_bottom_depth)/(TI_PEM(ig,2,islope)**2))))
+                            endif ! of if (iref /= 0)
+                        else
+                            if (iref /= 0) then ! mixed layer
+                                TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/                                  &
+                                                             (((icetable_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + &
+                                                             ((layer_PEM(iref + 1) - icetable_depth(ig,islope))/(ice_inertia**2))))
+                            else ! first layer is already a mixed layer
+                                ! (ie: take layer(iref=0)=0)
+                                TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/                                      &
+                                                      (((icetable_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) + &
+                                                      ((layer_PEM(1) - icetable_depth(ig,islope))/(ice_inertia**2))))
+                            endif ! of if (iref /= 0)
+                        endif ! iref == iend
+
+                        TI_PEM(ig,iref + 2:iend,islope) = ice_inertia
+                        if (iend < nsoil) then
+                            TI_PEM(ig,iend + 1,islope) = sqrt((layer_PEM(iend + 1) - layer_PEM(iend))/                         &
+                                                         (((ice_bottom_depth - layer_PEM(iend))/(TI_PEM(ig,iend,islope)**2)) + &
+                                                         ((layer_PEM(iend + 1) - ice_bottom_depth)/(ice_inertia**2))))
+                        endif
+                    endif ! of if (iref < nsoilmx)
+                else if (icetable_dynamic) then
+                    do  isoil = 1,nsoil
+                        call ice_thermal_properties(.false.,ice_porefilling(ig,isoil,islope),regolith_inertia(ig,islope),TI_PEM(ig,isoil,islope))
+                    enddo
+                endif ! of if icetable_equilibrium
             endif ! permanent glaciers
-        endif ! ice_depth(ig,islope) > 0.
-!        write(*,*) 'TI=', TI_PEM(ig,:,islope)
+        endif ! icetable_depth(ig,islope) > 0.
     enddo !islope
 enddo !ig