source: trunk/LMDZ.COMMON/libf/evolution/soil_therm_inertia.F90 @ 4174

Last change on this file since 4174 was 4135, checked in by jbclement, 5 weeks ago

PEM:
Relocate Mars-specific parameters out of "planet" module and into the modules that own their physics domain.
JBC
File size: 13.6 KB
Line 
1MODULE soil_therm_inertia
2!-----------------------------------------------------------------------
3! NAME
4!     soil_therm_inertia
5!
6! DESCRIPTION
7!     Compute and update soil thermal properties based on ice content,
8!     pressure, and cementation state.
9!
10! AUTHORS & DATE
11!     L. Lange, 2023
12!     JB Clement, 2025
13!
14! NOTES
15!
16!-----------------------------------------------------------------------
17
18! DEPENDENCIES
19! ------------
20use numerics, only: dp, di, k4, minieps
21
22! DECLARATION
23! -----------
24implicit none
25
26! PARAMETERS
27! ----------
28real(dp), parameter          :: TI_regolith_avg = 250._dp      ! Average of observed thermal inertia for regolith [SI], Mellon et al. 2000
29real(dp), parameter          :: TI_breccia = 750._dp           ! Thermal inertia for breccia [SI], Mellon et al. 2000; Wood et al. 2008
30real(dp), parameter          :: TI_bedrock = 2300._dp          ! Thermal inertia for bedrock [SI], Mellon et al. 2000; Wood et al. 2008
31real(dp), parameter, private :: P610 = 610._dp                 ! Reference mean surface pressure [Pa]
32real(dp), parameter, private :: reg_inertie_thresold = 370._dp ! Above this thermal inertia, the regolith has too much cementation to be dependant on the pressure [J/m^2/K/s^1/2]
33real(dp), parameter, private :: reg_inertie_minvalue = 50._dp  ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
34real(dp), parameter, private :: reg_inertie_maxvalue = 370._dp ! Maximum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
35real(dp), parameter, private :: C = 0.0015_dp                  ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [unitless]
36real(dp), parameter, private :: K = 8.1*1e4_dp                 ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [Torr, or 133.3Pa]
37real(dp), parameter, private :: Pa2Torr = 1./133.3_dp          ! Conversion from Pa to tor [Pa/Torr]
38
39contains
40!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
41
42!=======================================================================
43SUBROUTINE get_ice_TI(ispureice,pore_filling,surf_thermalinertia,ice_thermalinertia)
44!-----------------------------------------------------------------------
45! NAME
46!     get_ice_TI
47!
48! DESCRIPTION
49!     Compute ice thermal properties based on pore filling and purity.
50!
51! AUTHORS & DATE
52!     L. Lange, 2023
53!     JB Clement, 2025
54!
55! NOTES
56!     Uses Siegler et al. (2012) formula for pore-filling ice;
57!     uses Mellon et al. (2004) value for pure water ice.
58!-----------------------------------------------------------------------
59
60! DEPENDENCIES
61! ------------
62use soil, only: regolith_porosity
63
64! DECLARATION
65! -----------
66implicit none
67
68! ARGUMENTS
69! ---------
70logical(k4), intent(in)  :: ispureice           ! Boolean to check if ice is massive or just pore filling
71real(dp),    intent(in)  :: pore_filling        ! ice pore filling in each layer (1)
72real(dp),    intent(in)  :: surf_thermalinertia ! surface thermal inertia (J/m^2/K/s^1/2)
73real(dp),    intent(out) :: ice_thermalinertia  ! Thermal inertia of ice when present in the pore (J/m^2/K/s^1/2)
74
75! LOCAL VARIABLES
76! ---------------
77real(dp) :: inertie_purewaterice = 2100._dp ! 2050 is better according to my computations with the formula from Siegler et al., 2012, but let's follow Mellon et al. (2004)
78
79! CODE
80! ----
81if (ispureice) then
82    ice_thermalinertia = inertie_purewaterice
83else
84    ice_thermalinertia = sqrt(surf_thermalinertia**2 + regolith_porosity*pore_filling*inertie_purewaterice**2) ! Siegler et al., 2012
85end if
86
87END SUBROUTINE get_ice_TI
88!=======================================================================
89
90!=======================================================================
91SUBROUTINE update_soil_TI(h2o_ice,p_avg_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI)
92!-----------------------------------------------------------------------
93! NAME
94!     update_soil_TI
95!
96! DESCRIPTION
97!     Update soil thermal inertia based on ice table, regolith properties,
98!     and pressure-dependent cementation.
99!
100! AUTHORS & DATE
101!     L. Lange, 2023
102!     JB Clement, 2025
103!
104! NOTES
105!
106!-----------------------------------------------------------------------
107
108! DEPENDENCIES
109! ------------
110use geometry, only: ngrid, nslope, nsoil
111use soil,     only: volcapa, layer, inertiedat, depth_breccia, depth_bedrock, index_breccia, index_bedrock, reg_thprop_dependp
112use display,  only: print_msg, LVL_NFO
113
114! DECLARATION
115! -----------
116implicit none
117
118! ARGUMENTS
119! ---------
120real(dp),                   intent(in)    :: p_avg_new          ! Global average surface pressure [Pa]
121real(dp), dimension(:,:),   intent(in)    :: h2o_ice            ! Surface Water ice [kg/m^2]
122real(dp), dimension(:,:),   intent(in)    :: icetable_depth     ! Ice table depth [m]
123real(dp), dimension(:,:),   intent(in)    :: icetable_thickness ! Ice table thickness [m]
124real(dp), dimension(:,:,:), intent(in)    :: ice_porefilling    ! Ice porefilling [m^3/m^3]
125logical(k4),                intent(in)    :: icetable_equilibrium, icetable_dynamic ! Computing method for ice table
126real(dp), dimension(:,:,:), intent(inout) :: TI                 ! Soil Thermal Inertia [J/m^2/K/s^1/2]
127
128! LOCAL VARIABLES
129! ---------------
130integer(di)                       :: ig, islope, isoil, iref, iend ! Loop variables
131real(dp), dimension(ngrid,nslope) :: regolith_inertia              ! Thermal inertia of the regolith (first layer in the GCM) [J/m^2/K/s^1/2]
132real(dp)                          :: delta                         ! Difference of depth to compute the  thermal inertia in a mixed layer [m]
133real(dp)                          :: ice_bottom_depth              ! Bottom depth of the subsurface ice [m]
134real(dp)                          :: d_part                        ! Regolith particle size [micrometer]
135real(dp)                          :: ice_inertia                   ! Inertia of water ice [SI]
136
137! CODE
138! ----
139call print_msg("> Updating soil properties",LVL_NFO)
140
141! 1. Modification of the regolith thermal inertia.
142do islope = 1,nslope
143    regolith_inertia(:,islope) = inertiedat(:,1)
144    do ig = 1,ngrid
145        if (abs(h2o_ice(ig,islope)) < minieps) regolith_inertia(ig,islope) = TI_regolith_avg
146        if (reg_thprop_dependp) then
147            if (TI(ig,1,islope) < reg_inertie_thresold) then
148                d_part = (regolith_inertia(ig,islope)**2/(volcapa*C*(P610*Pa2Torr)**0.6))**(-1./(0.11_dp*log10(P610*Pa2Torr/K))) ! compute particle size, in micrometer
149                regolith_inertia(ig,islope) = sqrt(volcapa*C*(p_avg_new*Pa2Torr)**0.6*d_part**(-0.11_dp*log10(p_avg_new*Pa2Torr/K)))
150                if (regolith_inertia(ig,islope) > reg_inertie_maxvalue) regolith_inertia(ig,islope) = reg_inertie_maxvalue
151                if (regolith_inertia(ig,islope) < reg_inertie_minvalue) regolith_inertia(ig,islope) = reg_inertie_minvalue
152            end if
153        end if
154    end do
155end do
156
157! 2. Build new thermal inertia
158do islope = 1,nslope
159    do ig = 1,ngrid
160        do isoil = 1,index_breccia
161            TI(ig,isoil,islope) = regolith_inertia(ig,islope)
162        end do
163        if (regolith_inertia(ig,islope) < TI_breccia) then
164!!! transition
165            delta = depth_breccia
166            TI(ig,index_breccia + 1,islope) = sqrt((layer(index_breccia + 1) - layer(index_breccia))/                  &
167                                                  (((delta - layer(index_breccia))/(TI(ig,index_breccia,islope)**2)) + &
168                                                  ((layer(index_breccia + 1) - delta)/(TI_breccia**2))))
169            do isoil = index_breccia + 2,index_bedrock
170                TI(ig,isoil,islope) = TI_breccia
171            end do
172        else ! we keep the high ti values
173            do isoil = index_breccia + 1,index_bedrock
174                TI(ig,isoil,islope) = TI(ig,index_breccia,islope)
175            end do
176        end if ! TI PEM and breccia comparison
177!!! transition
178        delta = depth_bedrock
179        TI(ig,index_bedrock + 1,islope) = sqrt((layer(index_bedrock + 1) - layer(index_bedrock))/                  &
180                                              (((delta - layer(index_bedrock))/(TI(ig,index_bedrock,islope)**2)) + &
181                                              ((layer(index_bedrock + 1) - delta)/(TI_bedrock**2))))
182        do isoil = index_bedrock + 2,nsoil
183            TI(ig,isoil,islope) = TI_bedrock
184        end do
185    end do ! ig
186end do ! islope
187
188! 3. Build new TI in case of ice table
189do ig = 1,ngrid
190    do islope = 1,nslope
191        if (icetable_depth(ig,islope) > -1.e-6_dp) then
192        ! 3.0 Case where it is perennial ice
193            if (icetable_depth(ig,islope) < minieps) then
194                call get_ice_TI(.true.,1._dp,0._dp,ice_inertia)
195                do isoil = 1,nsoil
196                    TI(ig,isoil,islope) = ice_inertia
197                end do
198            else
199                if (icetable_equilibrium) then
200                    call get_ice_TI(.false.,1._dp,regolith_inertia(ig,islope),ice_inertia)
201                    ! 3.1.1 find the index of the mixed layer
202                    iref = 0 ! initialize iref
203                    do isoil = 1,nsoil ! loop on layers to find the beginning of the ice table
204                        if (icetable_depth(ig,islope) >= layer(isoil)) then
205                            iref = isoil ! pure regolith layer up to here
206                        else
207                            exit ! correct iref was obtained in previous cycle
208                        end if
209                    end do
210                    ! 3.1.2 find the index of the end of the ice table
211                    iend = 0 ! initialize iend
212                    ice_bottom_depth = icetable_depth(ig,islope) + icetable_thickness(ig,islope)
213                    do isoil = 1,nsoil ! loop on layers to find the end of the ice table
214                        if (ice_bottom_depth >= layer(isoil)) then
215                            iend = isoil ! pure regolith layer up to here
216                        else
217                            exit ! correct iref was obtained in previous cycle
218                        end if
219                    end do
220                    ! 3.2 Build the new ti
221                    if (iref < nsoil) then
222                        if (iref == iend) then
223                            ! Ice table begins and end in the same mixture with three components
224                            if (iref /= 0) then ! mixed layer
225                                TI(ig,iref + 1,islope) = sqrt((layer(iref + 1) - layer(iref))/                                      &
226                                                             (((icetable_depth(ig,islope) - layer(iref))/(TI(ig,iref,islope)**2)) + &
227                                                             ((ice_bottom_depth - icetable_depth(ig,islope))/(ice_inertia**2)) +    &
228                                                             ((layer(iref + 1) - ice_bottom_depth)/(TI(ig,iref + 1,islope)**2))))
229                            else ! first layer is already a mixed layer
230                                ! (ie: take layer(iref=0)=0)
231                                TI(ig,1,islope) = sqrt((layer(1))/                                                        &
232                                                      (((icetable_depth(ig,islope))/(TI(ig,1,islope)**2)) +               &
233                                                      ((ice_bottom_depth - icetable_depth(ig,islope))/(ice_inertia**2)) + &
234                                                      ((layer(2) - ice_bottom_depth)/(TI(ig,2,islope)**2))))
235                            end if ! of if (iref /= 0)
236                        else
237                            if (iref /= 0) then ! mixed layer
238                                TI(ig,iref + 1,islope) = sqrt((layer(iref + 1) - layer(iref))/                                      &
239                                                             (((icetable_depth(ig,islope) - layer(iref))/(TI(ig,iref,islope)**2)) + &
240                                                             ((layer(iref + 1) - icetable_depth(ig,islope))/(ice_inertia**2))))
241                            else ! first layer is already a mixed layer
242                                ! (ie: take layer(iref=0)=0)
243                                TI(ig,1,islope) = sqrt((layer(1))/                                          &
244                                                      (((icetable_depth(ig,islope))/(TI(ig,1,islope)**2)) + &
245                                                      ((layer(1) - icetable_depth(ig,islope))/(ice_inertia**2))))
246                            end if ! of if (iref /= 0)
247                        end if ! iref == iend
248
249                        TI(ig,iref + 2:iend,islope) = ice_inertia
250                        if (iend < nsoil) then
251                            TI(ig,iend + 1,islope) = sqrt((layer(iend + 1) - layer(iend))/                             &
252                                                         (((ice_bottom_depth - layer(iend))/(TI(ig,iend,islope)**2)) + &
253                                                         ((layer(iend + 1) - ice_bottom_depth)/(ice_inertia**2))))
254                        end if
255                    end if ! of if (iref < nsoil)
256                else if (icetable_dynamic) then
257                    do  isoil = 1,nsoil
258                        call get_ice_TI(.false.,ice_porefilling(ig,isoil,islope),regolith_inertia(ig,islope),TI(ig,isoil,islope))
259                    end do
260                end if ! of if icetable_equilibrium
261            end if ! permanent glaciers
262        end if ! icetable_depth(ig,islope) > 0.
263    end do !islope
264end do !ig
265
266END SUBROUTINE update_soil_TI
267!=======================================================================
268
269END MODULE soil_therm_inertia
Note: See TracBrowser for help on using the repository browser.