source: trunk/LMDZ.COMMON/libf/evolution/soil_thermal_inertia.F90 @ 4003

Last change on this file since 4003 was 3991, checked in by jbclement, 2 months ago

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