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

Last change on this file since 3296 was 3206, checked in by jbclement, 10 months ago

PEM:
Cleanings of unused variables/arguments and bad type conversions.
JBC
File size: 10.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, intent(in) :: tendencies_waterice(ngrid,nslope) ! Tendencies on the water ice [kg/m^2/year]
70real, intent(in) :: waterice(ngrid,nslope)            ! Surface Water ice [kg/m^2]
71real, intent(in) :: ice_depth(ngrid,nslope)           ! Ice table depth [m]
72real, intent(in) :: ice_thickness(ngrid,nslope)       ! Ice table thickness [m]
73! Outputs:
74
75real, intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope) ! Soil Thermal Inertia [J/m^2/K/s^1/2]
76 
77! Constants:
78
79 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]
80 REAL ::  reg_inertie_minvalue = 50.                       ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
81 REAL ::  reg_inertie_maxvalue = 370.                      ! Maximum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
82 REAL ::  ice_inertia                                      ! Inertia of water ice [SI]
83 REAL ::  P610 = 610.0                                     ! current average pressure on Mars [Pa]
84 REAL ::  C = 0.0015                                       ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [uniteless]
85 REAL ::  K = 8.1*1e4                                      ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [tor, or 133.3Pa]
86 REAL ::  Pa2Tor = 1./133.3                                ! Conversion from Pa to tor [Pa/tor]
87
88
89! Local variables:
90
91 INTEGER :: ig,islope,isoil,iref,iend                      ! Loop variables
92 REAL :: regolith_inertia(ngrid,nslope)                    ! Thermal inertia of the regolith (first layer in the GCM) [J/m^2/K/s^1/2]
93 REAL :: delta                                             ! Difference of depth to compute the  thermal inertia in a mixed layer [m]
94 REAL :: ice_bottom_depth                                  ! Bottom depth of the subsurface ice [m]
95 REAL :: d_part                                            ! Regolith particle size [micrometer]
96 
97
98write(*,*) "Update soil propreties"
99       
100! 1. Modification of the regolith thermal inertia.
101 do islope = 1,nslope
102   regolith_inertia(:,islope) = inertiedat_PEM(:,1)
103   do ig = 1,ngrid
104      if((tendencies_waterice(ig,islope).lt.-1e-5).and.(waterice(ig,islope).eq.0)) then
105              regolith_inertia(ig,islope) = TI_regolith_avg
106       endif
107      if(reg_thprop_dependp) then
108          if(TI_PEM(ig,1,islope).lt.reg_inertie_thresold) then
109            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
110            regolith_inertia(ig,islope) = sqrt(volcapa*C*(p_avg_new*Pa2Tor)**(0.6)*d_part**(-0.11*log10(p_avg_new*Pa2Tor/K)))
111            if(regolith_inertia(ig,islope).gt.reg_inertie_maxvalue) regolith_inertia(ig,islope) = reg_inertie_maxvalue
112            if(regolith_inertia(ig,islope).lt.reg_inertie_minvalue) regolith_inertia(ig,islope) = reg_inertie_minvalue
113          endif
114      endif
115   enddo
116 enddo
117
118! 2. Build new Thermal Inertia
119do  islope=1,nslope
120   do ig = 1,ngrid
121     do isoil = 1,index_breccia
122         TI_PEM(ig,isoil,islope) = regolith_inertia(ig,islope)
123     enddo
124     if(regolith_inertia(ig,islope).lt.TI_breccia) then
125!!! transition
126             delta = depth_breccia
127             TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
128            (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
129                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
130            do isoil=index_breccia+2,index_bedrock
131              TI_PEM(ig,isoil,islope) = TI_breccia
132            enddo     
133      else ! we keep the high ti values
134            do isoil=index_breccia+1,index_bedrock
135              TI_PEM(ig,isoil,islope) = TI_PEM(ig,index_breccia,islope)
136            enddo
137       endif ! TI PEM and breccia comparison
138!! transition
139       delta = depth_bedrock
140       TI_PEM(ig,index_bedrock+1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
141            (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
142            ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
143       do isoil=index_bedrock+2,nsoil_PEM
144            TI_PEM(ig,isoil,islope) = TI_bedrock
145       enddo   
146      enddo ! ig
147ENDDO ! islope
148
149!  3. Build new TI in case of ice table
150  do ig=1,ngrid
151    do islope=1,nslope
152     if (ice_depth(ig,islope).gt.-1e-6) then
153! 3.0 Case where it is perennial ice
154      if (ice_depth(ig,islope).lt. 1e-10) then
155       call ice_thermal_properties(.true.,1.,0.,ice_inertia)
156       do isoil = 1,nsoil_PEM
157         TI_PEM(ig,isoil,islope)=ice_inertia
158       enddo
159      else
160
161        call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),ice_inertia)
162        ! 3.1.1 find the index of the mixed layer
163        iref=0 ! initialize iref
164        do isoil=1,nsoil_PEM ! loop on layers to find the beginning of the ice table
165          if (ice_depth(ig,islope).ge.layer_PEM(isoil)) then
166            iref=isoil ! pure regolith layer up to here
167          else
168         ! correct iref was obtained in previous cycle
169            exit
170          endif
171        enddo
172        ! 3.1.2 find the index of the end of the ice table
173        iend=0 ! initialize iend
174        ice_bottom_depth = ice_depth(ig,islope)+ice_thickness(ig,islope)
175        do isoil=1,nsoil_PEM ! loop on layers to find the end of the ice table
176          if (ice_bottom_depth.ge.layer_PEM(isoil)) then
177            iend=isoil ! pure regolith layer up to here
178          else
179         ! correct iref was obtained in previous cycle
180            exit
181          endif
182        enddo
183      ! 3.2 Build the new ti
184        do isoil=1,iref
185          TI_PEM(ig,isoil,islope) = TI_PEM(ig,1,islope)
186        enddo
187        if (iref.lt.nsoil_PEM) then
188           if(iref.eq.iend) then
189            ! Ice table begins and end in the same mixture Mixtures with three components
190            if (iref.ne.0) then !
191            ! 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                      ((ice_bottom_depth-ice_depth(ig,islope))/(ice_inertia**2))+ &
195                      ((layer_PEM(iref+1)-ice_bottom_depth)/(TI_PEM(ig,iref+1,islope)**2))))
196            else ! first layer is already a mixed layer
197              ! (ie: take layer(iref=0)=0)
198              TI_PEM(ig,1,islope)=sqrt((layer_PEM(1))/ &
199                      (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2))+ &
200                      ((ice_bottom_depth-ice_depth(ig,islope))/(ice_inertia**2))+ &
201                      ((layer_PEM(2)-ice_bottom_depth)/(TI_PEM(ig,2,islope)**2))))
202            endif ! of if (iref.ne.0)
203           else
204            if (iref.ne.0) then
205            ! mixed layer
206              TI_PEM(ig,iref+1,islope)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ &
207                      (((ice_depth(ig,islope)-layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2))+ &
208                      ((layer_PEM(iref+1)-ice_depth(ig,islope))/(ice_inertia**2))))
209            else ! first layer is already a mixed layer
210              ! (ie: take layer(iref=0)=0)
211              TI_PEM(ig,1,islope)=sqrt((layer_PEM(1))/ &
212                      (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2))+ &
213                      ((layer_PEM(1)-ice_depth(ig,islope))/(ice_inertia**2))))
214            endif ! of if (iref.ne.0)       
215           endif ! iref.eq.iend
216      ! 3.3 Build the new ti
217          do isoil=iref+2,iend
218            TI_PEM(ig,isoil,islope)=ice_inertia
219          enddo
220
221
222          if(iend.lt.nsoil_PEM) then
223             TI_PEM(ig,iend+1,islope)=sqrt((layer_PEM(iend+1)-layer_PEM(iend))/ &
224                      (((ice_bottom_depth-layer_PEM(iend))/(TI_PEM(ig,iend,islope)**2)) + &
225                      ((layer_PEM(iend+1)-ice_bottom_depth)/(ice_inertia**2))))
226          endif
227
228
229        endif ! of if (iref.lt.(nsoilmx))
230      endif ! permanent glaciers
231     endif !ice_depth(ig,islope) > 0.
232!     write(*,*) 'TI=', TI_PEM(ig,:,islope)
233    enddo !islope
234   enddo !ig
235
236END SUBROUTINE update_soil_thermalproperties
237
238END MODULE soil_thermalproperties_mod
Note: See TracBrowser for help on using the repository browser.