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

Last change on this file since 3026 was 2985, checked in by romain.vande, 18 months ago

Generic PEM :

Adapt PEM to run with the generic model.
(CPP_STD keyword to exclude some part of the code at the compilation)

RV
File size: 9.4 KB
Line 
1module soil_thermalproperties_mod
2
3  implicit none
4
5!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
6!!!
7!!! Purpose: Compute the soil thermal properties
8!!! Author:  LL, 04/2023
9!!!
10!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
11
12 contains
13
14 subroutine ice_thermal_properties(ispureice,pore_filling,surf_thermalinertia,ice_thermalinertia)
15
16         
17
18!=======================================================================
19!   subject: Compute ice thermal properties
20!   --------
21!
22!   author: LL, 04/2023
23!   ------
24!         
25!=======================================================================
26
27 USE constants_marspem_mod,only: porosity
28
29      IMPLICIT NONE
30
31!-----------------------------------------------------------------------
32!=======================================================================
33!    Declarations :
34!=======================================================================
35!
36!    Input/Output
37!    ------------
38     LOGICAL,INTENT(IN) :: ispureice                         ! Boolean to check if ice is massive or just pore filling
39     REAL,INTENT(IN) ::    pore_filling                         ! ice pore filling in each layer (1)
40     REAL,INTENT(IN) ::    surf_thermalinertia                  ! surface thermal inertia (J/m^2/K/s^1/2) 
41     REAL,INTENT(OUT) ::   ice_thermalinertia                   ! Thermal inertia of ice when present in the pore (J/m^2/K/s^1/2)
42
43!    Local Variables
44!    --------------
45     REAL :: 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)
46!=======================================================================
47!   Beginning of the code
48!=======================================================================
49
50   if(ispureice) then
51           ice_thermalinertia = inertie_purewaterice
52   else
53           ice_thermalinertia = sqrt(surf_thermalinertia**2 + porosity*pore_filling*inertie_purewaterice**2)  ! Siegler et al.,2012
54   endif
55
56 end subroutine
57!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
58
59  SUBROUTINE update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tendencies_waterice,waterice,p_avg_new,ice_depth,ice_thickness,TI_PEM)
60
61 USE comsoil_h, only:  inertiedat, volcapa
62 USE comsoil_h_PEM, only: layer_PEM,inertiedat_PEM,depth_breccia,depth_bedrock,index_breccia,index_bedrock,reg_thprop_dependp
63 USE vertical_layers_mod, ONLY: ap,bp
64 USE constants_marspem_mod,only: TI_breccia,TI_bedrock, TI_regolith_avg
65 implicit none
66! Input:
67 INTEGER,INTENT(IN) ::  ngrid, nslope, nsoil_PEM           ! Shape of the arrays: physical grid, number of sub-grid slopes, number of layer in the soil
68 REAL,INTENT(IN) :: p_avg_new                              ! Global average surface pressure [Pa]
69 REAL,INTENT(IN) :: tendencies_waterice(ngrid,nslope)      ! Tendencies on the water ice [kg/m^2/year]
70 REAL,INTENT(IN) :: waterice(ngrid,nslope)                 ! Surface Water ice [kg/m^2]
71 REAL,INTENT(in) :: ice_depth(ngrid,nslope)                ! Ice table depth [m]
72 REAL,INTENT(in) :: ice_thickness(ngrid,nslope)            ! Ice table thickness [m]
73! Outputs:
74
75 REAL,INTENT(INOUT) :: TI_PEM(ngrid,nsoil_PEM,nslope)      ! Soil Thermal Inertia [J/m^2/K/s^1/2]
76 
77! Constants:
78
79 REAL ::  inertie_thresold = 800. ! Above this thermal inertia, it's  ice [SI]
80 REAL ::  ice_inertia        ! Inertia of water ice [SI]
81 REAL ::  P610 = 610.0       ! current average pressure on Mars [Pa]
82
83! Local variables:
84
85 INTEGER :: ig,islope,iloop,iref,k,iend
86 REAL :: regolith_inertia(ngrid,nslope) ! TI of the regolith
87 REAL :: delta
88 REAL :: TI_breccia_new
89 REAL :: ice_bottom_depth
90
91! 1. Modification of the regolith thermal inertia.
92
93
94
95 do islope = 1,nslope
96   regolith_inertia(:,islope) = inertiedat_PEM(:,1)
97   do ig = 1,ngrid
98      if((tendencies_waterice(ig,islope).lt.-1e-5).and.(waterice(ig,islope).eq.0)) then
99              regolith_inertia(ig,islope) = TI_regolith_avg
100       endif
101      if(reg_thprop_dependp) then
102          if(TI_PEM(ig,1,islope).lt.inertie_thresold) then
103            regolith_inertia(ig,islope) = regolith_inertia(ig,islope)*(p_avg_new/P610)**0.3
104          endif
105         TI_breccia_new = TI_breccia*(p_avg_new/P610)**0.3
106      else
107         TI_breccia_new = TI_breccia
108      endif
109   enddo
110 enddo
111
112
113! 2. Build new Thermal Inertia
114
115do  islope=1,nslope
116   do ig = 1,ngrid
117     do iloop = 1,index_breccia
118         TI_PEM(ig,iloop,islope) = regolith_inertia(ig,islope)
119     enddo
120     if(regolith_inertia(ig,islope).lt.TI_breccia_new) 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_new**2))))
126            do iloop=index_breccia+2,index_bedrock
127              TI_PEM(ig,iloop,islope) = TI_breccia_new
128            enddo     
129      else ! we keep the high ti values
130            do iloop=index_breccia+1,index_bedrock
131              TI_PEM(ig,iloop,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 iloop=index_bedrock+2,nsoil_PEM
140            TI_PEM(ig,iloop,islope) = TI_bedrock
141       enddo   
142      enddo ! ig
143ENDDO ! islope
144
145!  3. Build new TI in case of ice table
146  do ig=1,ngrid
147    do islope=1,nslope
148     if (ice_depth(ig,islope).gt.-1e-6) then
149! 3.0 Case where it is perenial ice
150      if (ice_depth(ig,islope).lt. 1e-10) then
151       call ice_thermal_properties(.true.,1.,0.,ice_inertia)
152       do iloop = 1,nsoil_PEM
153         TI_PEM(ig,iloop,islope)=ice_inertia
154       enddo
155      else
156
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 k=1,nsoil_PEM ! loop on layers to find the beginning of the ice table
161          if (ice_depth(ig,islope).ge.layer_PEM(k)) then
162            iref=k ! pure regolith layer up to here
163          else
164         ! correct iref was obtained in previous cycle
165            exit
166          endif
167        enddo
168        ! 3.1.2 find the index of the end of the ice table
169        iend=0 ! initialize iend
170        ice_bottom_depth = ice_depth(ig,islope)+ice_thickness(ig,islope)
171        do k=1,nsoil_PEM ! loop on layers to find the end of the ice table
172          if (ice_bottom_depth.ge.layer_PEM(k)) then
173            iend=k ! pure regolith layer up to here
174          else
175         ! correct iref was obtained in previous cycle
176            exit
177          endif
178        enddo
179      ! 3.2 Build the new ti
180        do iloop=1,iref
181          TI_PEM(ig,iloop,islope) = TI_PEM(ig,1,islope)
182        enddo
183        if (iref.lt.nsoil_PEM) then
184           if(iref.eq.iend) then
185            ! Ice table begins and end in the same mixture Mixtures with three components
186            if (iref.ne.0) then !
187            ! mixed layer
188              TI_PEM(ig,iref+1,islope)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ &
189                      (((ice_depth(ig,islope)-layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2))+ &
190                      ((ice_bottom_depth-ice_depth(ig,islope))/(ice_inertia**2))+ &
191                      ((layer_PEM(iref+1)-ice_bottom_depth)/(TI_PEM(ig,iref+1,islope)**2))))
192            else ! first layer is already a mixed layer
193              ! (ie: take layer(iref=0)=0)
194              TI_PEM(ig,1,islope)=sqrt((layer_PEM(1))/ &
195                      (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2))+ &
196                      ((ice_bottom_depth-ice_depth(ig,islope))/(ice_inertia**2))+ &
197                      ((layer_PEM(2)-ice_bottom_depth)/(TI_PEM(ig,2,islope)**2))))
198            endif ! of if (iref.ne.0)
199           else
200            if (iref.ne.0) then
201            ! mixed layer
202              TI_PEM(ig,iref+1,islope)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ &
203                      (((ice_depth(ig,islope)-layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2))+ &
204                      ((layer_PEM(iref+1)-ice_depth(ig,islope))/(ice_inertia**2))))
205            else ! first layer is already a mixed layer
206              ! (ie: take layer(iref=0)=0)
207              TI_PEM(ig,1,islope)=sqrt((layer_PEM(1))/ &
208                      (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2))+ &
209                      ((layer_PEM(1)-ice_depth(ig,islope))/(ice_inertia**2))))
210            endif ! of if (iref.ne.0)       
211           endif ! iref.eq.iend
212      ! 3.3 Build the new ti
213          do iloop=iref+2,iend
214            TI_PEM(ig,iloop,islope)=ice_inertia
215          enddo
216
217
218          if(iend.lt.nsoil_PEM) then
219             TI_PEM(ig,iend+1,islope)=sqrt((layer_PEM(iend+1)-layer_PEM(iend))/ &
220                      (((ice_bottom_depth-layer_PEM(iend))/(TI_PEM(ig,iend,islope)**2)) + &
221                      ((layer_PEM(iend+1)-ice_bottom_depth)/(ice_inertia**2))))
222          endif
223
224
225        endif ! of if (iref.lt.(nsoilmx))
226      endif ! permanent glaciers
227     endif !ice_depth(ig,islope) > 0.
228!     write(*,*) 'TI=', TI_PEM(ig,:,islope)
229    enddo !islope
230   enddo !ig
231
232!=======================================================================
233      RETURN
234      END subroutine update_soil_thermalproperties
235
236end module soil_thermalproperties_mod
237
Note: See TracBrowser for help on using the repository browser.