source: trunk/LMDZ.COMMON/libf/evolution/soil_thermalproperties_mod.F90 @ 3525

Last change on this file since 3525 was 3525, checked in by jbclement, 6 days ago

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
File size: 12.8 KB
Line 
1MODULE soil_thermalproperties_mod
2
3implicit none
4
5!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
6!!!
7!!! Purpose: Compute the soil thermal properties
8!!! Author:  LL, 04/2023
9!!!
10!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
11
12! Constants:
13real, 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]
14real, parameter :: reg_inertie_minvalue = 50.  ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
15real, parameter :: reg_inertie_maxvalue = 370. ! Maximum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
16real, parameter :: P610 = 610.0                ! current average pressure on Mars [Pa]
17real, parameter :: C = 0.0015                  ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [unitless]
18real, parameter :: K = 8.1*1e4                 ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [Torr, or 133.3Pa]
19real, parameter :: Pa2Torr = 1./133.3          ! Conversion from Pa to tor [Pa/Torr]
20
21!=======================================================================
22contains
23!=======================================================================
24
25SUBROUTINE ice_thermal_properties(ispureice,pore_filling,surf_thermalinertia,ice_thermalinertia)
26!=======================================================================
27!   subject: Compute ice thermal properties
28!   --------
29!
30!   author: LL, 04/2023
31!   -------
32!
33!=======================================================================
34
35use constants_marspem_mod, only: porosity
36
37implicit none
38
39!-----------------------------------------------------------------------
40!=======================================================================
41!    Declarations :
42!=======================================================================
43!
44!    Input/Output
45!    ------------
46logical, intent(in) :: ispureice           ! Boolean to check if ice is massive or just pore filling
47real,    intent(in) :: pore_filling        ! ice pore filling in each layer (1)
48real,    intent(in) :: surf_thermalinertia ! surface thermal inertia (J/m^2/K/s^1/2)
49real,    intent(out) :: ice_thermalinertia ! Thermal inertia of ice when present in the pore (J/m^2/K/s^1/2)
50
51!    Local Variables
52!    --------------
53REAL :: inertie_purewaterice = 2100 ! 2050 is better according to my computations with the formula from Siegler et al., 2012, but let's follow Mellon et al. (2004)
54!=======================================================================
55!   Beginning of the code
56!=======================================================================
57
58if (ispureice) then
59    ice_thermalinertia = inertie_purewaterice
60else
61    ice_thermalinertia = sqrt(surf_thermalinertia**2 + porosity*pore_filling*inertie_purewaterice**2) ! Siegler et al., 2012
62endif
63
64END SUBROUTINE
65!=======================================================================
66
67SUBROUTINE update_soil_thermalproperties(ngrid,nslope,nsoil,tendencies_waterice,waterice,p_avg_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
68
69use comsoil_h,             only: volcapa
70use comsoil_h_PEM,         only: layer_PEM, inertiedat_PEM, depth_breccia, depth_bedrock, index_breccia, index_bedrock, reg_thprop_dependp
71use constants_marspem_mod, only: TI_breccia, TI_bedrock, TI_regolith_avg
72
73implicit none
74
75! Input:
76integer,                             intent(in) :: ngrid, nslope, nsoil ! Shape of the arrays: physical grid, number of sub-grid slopes, number of layer in the soil
77real,                                intent(in) :: p_avg_new            ! Global average surface pressure [Pa]
78real, dimension(ngrid,nslope),       intent(in) :: tendencies_waterice  ! Tendencies on the water ice [kg/m^2/year]
79real, dimension(ngrid,nslope),       intent(in) :: waterice             ! Surface Water ice [kg/m^2]
80real, dimension(ngrid,nslope),       intent(in) :: icetable_depth       ! Ice table depth [m]
81real, dimension(ngrid,nslope),       intent(in) :: icetable_thickness   ! Ice table thickness [m]
82real, dimension(ngrid,nsoil,nslope), intent(in) :: ice_porefilling      ! Ice porefilling [m^3/m^3]
83logical,                             intent(in) :: icetable_equilibrium, icetable_dynamic ! Computing method for ice table
84
85! Outputs:
86real, dimension(ngrid,nsoil,nslope), intent(inout) :: TI_PEM ! Soil Thermal Inertia [J/m^2/K/s^1/2]
87
88! Local variables:
89integer :: ig, islope, isoil, iref, iend          ! Loop variables
90real, dimension(ngrid,nslope) :: regolith_inertia ! Thermal inertia of the regolith (first layer in the GCM) [J/m^2/K/s^1/2]
91real    :: delta                                  ! Difference of depth to compute the  thermal inertia in a mixed layer [m]
92real    :: ice_bottom_depth                       ! Bottom depth of the subsurface ice [m]
93real    :: d_part                                 ! Regolith particle size [micrometer]
94real    :: ice_inertia                            ! Inertia of water ice [SI]
95
96write(*,*) "Update soil properties"
97
98! 1. Modification of the regolith thermal inertia.
99do islope = 1,nslope
100    regolith_inertia(:,islope) = inertiedat_PEM(:,1)
101    do ig = 1,ngrid
102        if (tendencies_waterice(ig,islope) < -1.e-5 .and. waterice(ig,islope) == 0) regolith_inertia(ig,islope) = TI_regolith_avg
103        if (reg_thprop_dependp) then
104            if (TI_PEM(ig,1,islope) < reg_inertie_thresold) then
105                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
106                regolith_inertia(ig,islope) = sqrt(volcapa*C*(p_avg_new*Pa2Torr)**(0.6)*d_part**(-0.11*log10(p_avg_new*Pa2Torr/K)))
107                if (regolith_inertia(ig,islope) > reg_inertie_maxvalue) regolith_inertia(ig,islope) = reg_inertie_maxvalue
108                if (regolith_inertia(ig,islope) < reg_inertie_minvalue) regolith_inertia(ig,islope) = reg_inertie_minvalue
109            endif
110        endif
111    enddo
112enddo
113
114! 2. Build new Thermal Inertia
115do islope = 1,nslope
116    do ig = 1,ngrid
117        do isoil = 1,index_breccia
118            TI_PEM(ig,isoil,islope) = regolith_inertia(ig,islope)
119        enddo
120        if (regolith_inertia(ig,islope) < TI_breccia) then
121!!! transition
122            delta = depth_breccia
123            TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/              &
124                                                  (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
125                                                  ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
126            do isoil = index_breccia + 2,index_bedrock
127                TI_PEM(ig,isoil,islope) = TI_breccia
128            enddo
129        else ! we keep the high ti values
130            do isoil = index_breccia + 1,index_bedrock
131                TI_PEM(ig,isoil,islope) = TI_PEM(ig,index_breccia,islope)
132            enddo
133        endif ! TI PEM and breccia comparison
134!!! transition
135        delta = depth_bedrock
136        TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/              &
137                                              (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
138                                              ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
139        do isoil = index_bedrock + 2,nsoil
140            TI_PEM(ig,isoil,islope) = TI_bedrock
141        enddo
142    enddo ! ig
143enddo ! islope
144
145! 3. Build new TI in case of ice table
146do ig = 1,ngrid
147    do islope = 1,nslope
148        if (icetable_depth(ig,islope) > -1.e-6) then
149        ! 3.0 Case where it is perennial ice
150            if (icetable_depth(ig,islope) < 1.e-10) then
151                call ice_thermal_properties(.true.,1.,0.,ice_inertia)
152                do isoil = 1,nsoil
153                    TI_PEM(ig,isoil,islope) = ice_inertia
154                enddo
155            else
156                if (icetable_equilibrium) then
157                    call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),ice_inertia)
158                    ! 3.1.1 find the index of the mixed layer
159                    iref = 0 ! initialize iref
160                    do isoil = 1,nsoil ! loop on layers to find the beginning of the ice table
161                        if (icetable_depth(ig,islope) >= layer_PEM(isoil)) then
162                            iref = isoil ! pure regolith layer up to here
163                        else
164                            exit ! correct iref was obtained in previous cycle
165                        endif
166                    enddo
167                    ! 3.1.2 find the index of the end of the ice table
168                    iend = 0 ! initialize iend
169                    ice_bottom_depth = icetable_depth(ig,islope) + icetable_thickness(ig,islope)
170                    do isoil = 1,nsoil ! loop on layers to find the end of the ice table
171                        if (ice_bottom_depth >= layer_PEM(isoil)) then
172                            iend = isoil ! pure regolith layer up to here
173                        else
174                            exit ! correct iref was obtained in previous cycle
175                        endif
176                    enddo
177                    ! 3.2 Build the new ti
178                    if (iref < nsoil) then
179                        if (iref == iend) then
180                            ! Ice table begins and end in the same mixture with three components
181                            if (iref /= 0) then ! mixed layer
182                                TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/                                  &
183                                                             (((icetable_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + &
184                                                             ((ice_bottom_depth - icetable_depth(ig,islope))/(ice_inertia**2)) +            &
185                                                             ((layer_PEM(iref + 1) - ice_bottom_depth)/(TI_PEM(ig,iref + 1,islope)**2))))
186                            else ! first layer is already a mixed layer
187                                ! (ie: take layer(iref=0)=0)
188                                TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/                                                &
189                                                      (((icetable_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) +           &
190                                                      ((ice_bottom_depth - icetable_depth(ig,islope))/(ice_inertia**2)) + &
191                                                      ((layer_PEM(2) - ice_bottom_depth)/(TI_PEM(ig,2,islope)**2))))
192                            endif ! of if (iref /= 0)
193                        else
194                            if (iref /= 0) then ! mixed layer
195                                TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/                                  &
196                                                             (((icetable_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + &
197                                                             ((layer_PEM(iref + 1) - icetable_depth(ig,islope))/(ice_inertia**2))))
198                            else ! first layer is already a mixed layer
199                                ! (ie: take layer(iref=0)=0)
200                                TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/                                      &
201                                                      (((icetable_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) + &
202                                                      ((layer_PEM(1) - icetable_depth(ig,islope))/(ice_inertia**2))))
203                            endif ! of if (iref /= 0)
204                        endif ! iref == iend
205
206                        TI_PEM(ig,iref + 2:iend,islope) = ice_inertia
207                        if (iend < nsoil) then
208                            TI_PEM(ig,iend + 1,islope) = sqrt((layer_PEM(iend + 1) - layer_PEM(iend))/                         &
209                                                         (((ice_bottom_depth - layer_PEM(iend))/(TI_PEM(ig,iend,islope)**2)) + &
210                                                         ((layer_PEM(iend + 1) - ice_bottom_depth)/(ice_inertia**2))))
211                        endif
212                    endif ! of if (iref < nsoilmx)
213                else if (icetable_dynamic) then
214                    do  isoil = 1,nsoil
215                        call ice_thermal_properties(.false.,ice_porefilling(ig,isoil,islope),regolith_inertia(ig,islope),TI_PEM(ig,isoil,islope))
216                    enddo
217                endif ! of if icetable_equilibrium
218            endif ! permanent glaciers
219        endif ! icetable_depth(ig,islope) > 0.
220    enddo !islope
221enddo !ig
222
223END SUBROUTINE update_soil_thermalproperties
224
225END MODULE soil_thermalproperties_mod
Note: See TracBrowser for help on using the repository browser.