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

Last change on this file since 3516 was 3347, checked in by llange, 6 months ago

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