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

Last change on this file since 3154 was 3149, checked in by jbclement, 2 years ago

PEM:

  • Simplification of the algorithm managing the stopping criteria;
  • Complete rework of the ice management in the PEM (H2O & CO2);

    Subroutines to evolve the H2O and CO2 ice are now in the same module "evol_ice_mod.F90".
    Tendencies are computed from the variation of "ice + frost" between the 2 PCM runs.
    Evolving ice in the PEM is now called 'h2o_ice' or 'co2_ice' (not anymore in 'qsurf' and free of 'water_reservoir').
    Default value 'ini_h2o_bigreservoir' (= 10 m) initializes the H2O ice of the first PEM run where there is 'watercap'. For the next PEM runs, initialization is done with the value kept in "startpem.nc". CO2 ice is taken from 'perennial_co2ice' of the PCM (paleoclimate flag must be true).
    Simplification of the condition to compute the surface ice cover needed for the stopping criteria.
    Frost ('qsurf') is not evolved by the PEM and given back to the PCM.
    New default threshold value 'inf_h2oice_threshold' (= 2 m) to decide at the end of the PEM run if the H2O ice should be 'watercap' or not for the next PCM runs. If H2O ice cannot be 'watercap', then the remaining H2O ice is transferred to the frost ('qsurf').

  • Renaming of variables/subroutines for clarity;
  • Some cleanings throughout the code;
  • Small updates in files of the deftank.

JBC
File size: 9.5 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
91write(*,*) "Update soil propreties"
92       
93! 1. Modification of the regolith thermal inertia.
94 do islope = 1,nslope
95   regolith_inertia(:,islope) = inertiedat_PEM(:,1)
96   do ig = 1,ngrid
97      if((tendencies_waterice(ig,islope).lt.-1e-5).and.(waterice(ig,islope).eq.0)) then
98              regolith_inertia(ig,islope) = TI_regolith_avg
99       endif
100      if(reg_thprop_dependp) then
101          if(TI_PEM(ig,1,islope).lt.inertie_thresold) then
102            regolith_inertia(ig,islope) = regolith_inertia(ig,islope)*(p_avg_new/P610)**0.3
103          endif
104         TI_breccia_new = TI_breccia*(p_avg_new/P610)**0.3
105      else
106         TI_breccia_new = TI_breccia
107      endif
108   enddo
109 enddo
110
111
112! 2. Build new Thermal Inertia
113do  islope=1,nslope
114   do ig = 1,ngrid
115     do iloop = 1,index_breccia
116         TI_PEM(ig,iloop,islope) = regolith_inertia(ig,islope)
117     enddo
118     if(regolith_inertia(ig,islope).lt.TI_breccia_new) 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_new**2))))
124            do iloop=index_breccia+2,index_bedrock
125              TI_PEM(ig,iloop,islope) = TI_breccia_new
126            enddo     
127      else ! we keep the high ti values
128            do iloop=index_breccia+1,index_bedrock
129              TI_PEM(ig,iloop,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 iloop=index_bedrock+2,nsoil_PEM
138            TI_PEM(ig,iloop,islope) = TI_bedrock
139       enddo   
140      enddo ! ig
141ENDDO ! islope
142
143!  3. Build new TI in case of ice table
144  do ig=1,ngrid
145    do islope=1,nslope
146     if (ice_depth(ig,islope).gt.-1e-6) then
147! 3.0 Case where it is perennial ice
148      if (ice_depth(ig,islope).lt. 1e-10) then
149       call ice_thermal_properties(.true.,1.,0.,ice_inertia)
150       do iloop = 1,nsoil_PEM
151         TI_PEM(ig,iloop,islope)=ice_inertia
152       enddo
153      else
154
155        call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),ice_inertia)
156        ! 3.1.1 find the index of the mixed layer
157        iref=0 ! initialize iref
158        do k=1,nsoil_PEM ! loop on layers to find the beginning of the ice table
159          if (ice_depth(ig,islope).ge.layer_PEM(k)) then
160            iref=k ! pure regolith layer up to here
161          else
162         ! correct iref was obtained in previous cycle
163            exit
164          endif
165        enddo
166        ! 3.1.2 find the index of the end of the ice table
167        iend=0 ! initialize iend
168        ice_bottom_depth = ice_depth(ig,islope)+ice_thickness(ig,islope)
169        do k=1,nsoil_PEM ! loop on layers to find the end of the ice table
170          if (ice_bottom_depth.ge.layer_PEM(k)) then
171            iend=k ! pure regolith layer up to here
172          else
173         ! correct iref was obtained in previous cycle
174            exit
175          endif
176        enddo
177      ! 3.2 Build the new ti
178        do iloop=1,iref
179          TI_PEM(ig,iloop,islope) = TI_PEM(ig,1,islope)
180        enddo
181        if (iref.lt.nsoil_PEM) then
182           if(iref.eq.iend) then
183            ! Ice table begins and end in the same mixture Mixtures with three components
184            if (iref.ne.0) then !
185            ! mixed layer
186              TI_PEM(ig,iref+1,islope)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ &
187                      (((ice_depth(ig,islope)-layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2))+ &
188                      ((ice_bottom_depth-ice_depth(ig,islope))/(ice_inertia**2))+ &
189                      ((layer_PEM(iref+1)-ice_bottom_depth)/(TI_PEM(ig,iref+1,islope)**2))))
190            else ! first layer is already a mixed layer
191              ! (ie: take layer(iref=0)=0)
192              TI_PEM(ig,1,islope)=sqrt((layer_PEM(1))/ &
193                      (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2))+ &
194                      ((ice_bottom_depth-ice_depth(ig,islope))/(ice_inertia**2))+ &
195                      ((layer_PEM(2)-ice_bottom_depth)/(TI_PEM(ig,2,islope)**2))))
196            endif ! of if (iref.ne.0)
197           else
198            if (iref.ne.0) then
199            ! mixed layer
200              TI_PEM(ig,iref+1,islope)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ &
201                      (((ice_depth(ig,islope)-layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2))+ &
202                      ((layer_PEM(iref+1)-ice_depth(ig,islope))/(ice_inertia**2))))
203            else ! first layer is already a mixed layer
204              ! (ie: take layer(iref=0)=0)
205              TI_PEM(ig,1,islope)=sqrt((layer_PEM(1))/ &
206                      (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2))+ &
207                      ((layer_PEM(1)-ice_depth(ig,islope))/(ice_inertia**2))))
208            endif ! of if (iref.ne.0)       
209           endif ! iref.eq.iend
210      ! 3.3 Build the new ti
211          do iloop=iref+2,iend
212            TI_PEM(ig,iloop,islope)=ice_inertia
213          enddo
214
215
216          if(iend.lt.nsoil_PEM) then
217             TI_PEM(ig,iend+1,islope)=sqrt((layer_PEM(iend+1)-layer_PEM(iend))/ &
218                      (((ice_bottom_depth-layer_PEM(iend))/(TI_PEM(ig,iend,islope)**2)) + &
219                      ((layer_PEM(iend+1)-ice_bottom_depth)/(ice_inertia**2))))
220          endif
221
222
223        endif ! of if (iref.lt.(nsoilmx))
224      endif ! permanent glaciers
225     endif !ice_depth(ig,islope) > 0.
226!     write(*,*) 'TI=', TI_PEM(ig,:,islope)
227    enddo !islope
228   enddo !ig
229
230!=======================================================================
231      RETURN
232      END subroutine update_soil_thermalproperties
233
234end module soil_thermalproperties_mod
235
Note: See TracBrowser for help on using the repository browser.