source: trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90 @ 2961

Last change on this file since 2961 was 2961, checked in by llange, 20 months ago

PEM

  • Move math functions (findroot, deriv, etc) in a specific module
  • Ice table is no longer infinite, the bottom of the ice table (zend) is set such at d(rho_vap) dz (zend) = 0. Future commit will transfer the loss of ice table to perenial glaciers

LL

File size: 14.3 KB
Line 
1module ice_table_mod
2
3  implicit none
4
5!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
6!!!
7!!! Purpose: Ice table (pore-filling) variables and methods to compute it (dynamic and static)
8!!! Author:  LL, 02/2023
9!!!
10!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
11
12  LOGICAL,SAVE :: icetable_equilibrium                    ! Boolean to say if the PEM needs to recompute the icetable depth when at  equilibrium
13  LOGICAL,SAVE :: icetable_dynamic                        ! Boolean to say if the PEM needs to recompute the icetable depth (dynamic method)
14  real,save,allocatable  :: porefillingice_depth(:,:)    ! ngrid x nslope: Depth of the ice table [m]
15  real,save,allocatable  :: porefillingice_thickness(:,:)    ! ngrid x nslope: Thickness of the ice table [m]
16
17  contains
18
19
20  subroutine ini_ice_table_porefilling(ngrid,nslope)
21 
22    implicit none
23    integer,intent(in) :: ngrid ! number of atmospheric columns
24    integer,intent(in) :: nslope ! number of slope within a mesh
25
26    allocate(porefillingice_depth(ngrid,nslope))
27    allocate(porefillingice_thickness(ngrid,nslope))
28
29  end subroutine ini_ice_table_porefilling
30
31  subroutine end_ice_table_porefilling
32
33    implicit none
34    if (allocated(porefillingice_depth)) deallocate(porefillingice_depth)
35    if (allocated(porefillingice_thickness)) deallocate(porefillingice_thickness)
36
37  end subroutine end_ice_table_porefilling
38
39
40
41
42   SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness)
43!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
44!!!
45!!! Purpose: Compute the ice table depth (pore-filling) knowing the yearly average water
46!!! density at the surface and at depth.
47!!! Computations are made following the methods in Schorgofer et al., 2005
48!!! This subroutine only gives the ice table at equilibrium and does not consider exchange with the atmosphere
49!!!
50!!! Author: LL
51!!!
52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
53    use math_mod,only: findroot
54#ifndef CPP_STD
55    USE comsoil_h_PEM, only: mlayer_PEM,layer_PEM                  ! Depth of the vertical grid
56    USE soil_thermalproperties_mod, only: ice_thermal_properties
57    implicit none
58!   inputs
59! -----------
60    integer,intent(in) :: ngrid,nslope,nsoil_PEM                   ! Size of the physical grid, number of subslope, number of soil layer in the PEM
61    logical,intent(in) :: watercaptag(ngrid)                       ! Boolean to check the presence of a perennial glacier
62    real,intent(in) :: rhowatersurf_ave(ngrid,nslope)              ! Water density at the surface, yearly averaged [kg/m^3]
63    real,intent(in) :: rhowatersoil_ave(ngrid,nsoil_PEM,nslope)    ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged  [kg/m^3]
64    real,intent(in) :: regolith_inertia(ngrid,nslope)              ! Thermal inertia of the regolith layer [SI]
65!   Ouputs
66! -----------
67    real,intent(out) :: ice_table_beg(ngrid,nslope)               ! ice table depth [m]
68    real,intent(out) :: ice_table_thickness(ngrid,nslope)         ! ice table thickness [m]
69!   Local
70! -----------
71    real :: z1,z2                                                  ! intermediate variables used when doing a linear interpolation between two depths to find the root
72    integer ig, islope,isoil,isoilend                              ! loop variables
73    real :: diff_rho(nsoil_PEM)                                    ! difference of water vapor density between the surface and at depth [kg/m^3]
74    real  :: ice_table_end                                         ! depth of the end of the ice table  [m]
75    real :: previous_icetable_depth(ngrid,nslope)                  ! Ice table computed at previous ice depth [m]
76    real :: stretch                                                ! stretch factor to improve the convergence of the ice table
77    real :: wice_inertia                                           ! Water Ice thermal Inertia [USI]   
78!   Code
79! -----------
80
81     previous_icetable_depth(:,:) = ice_table_beg(:,:)
82     do ig = 1,ngrid
83      if(watercaptag(ig)) then
84         ice_table_beg(ig,:) = 0.
85         ice_table_thickness(ig,:) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
86      else
87        do islope = 1,nslope
88           ice_table_beg(ig,islope) = -1.
89           ice_table_thickness(ig,islope) = 0.
90         do isoil = 1,nsoil_PEM
91           diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope)
92         enddo
93         if(diff_rho(1) > 0) then                                   ! ice is at the surface
94           ice_table_beg(ig,islope) = 0.
95               do isoilend = 2,nsoil_PEM-1
96                 if((diff_rho(isoilend).gt.0).and.(diff_rho(isoilend+1).lt.0.)) then
97                  call findroot(diff_rho(isoilend),diff_rho(isoilend+1),mlayer_PEM(isoilend),mlayer_PEM(isoilend+1),ice_table_end)
98                  ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope)
99                  exit
100                 endif
101               enddo
102         else
103           do isoil = 1,nsoil_PEM -1                                ! general case, we find the ice table depth by doing a linear approximation between the two depth, and then solve the first degree equation to find the root
104             if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then
105               call findroot(diff_rho(isoil),diff_rho(isoil+1),mlayer_PEM(isoil),mlayer_PEM(isoil+1),ice_table_beg(ig,islope))
106               ! Now let's find the end of the ice table
107               ice_table_thickness(ig,islope) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
108               do isoilend = isoil+1,nsoil_PEM-1
109                 if((diff_rho(isoilend).gt.0).and.(diff_rho(isoilend+1).lt.0.)) then
110                  call findroot(diff_rho(isoilend),diff_rho(isoilend+1),mlayer_PEM(isoilend),mlayer_PEM(isoilend+1),ice_table_end)
111                  ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope)
112                  exit
113                 endif
114               enddo
115             endif !diff_rho(z) <0 & diff_rho(z+1) > 0
116           enddo
117          endif ! diff_rho(1) > 0
118        enddo
119       endif ! watercaptag
120      enddo
121
122! Small trick to speed up the convergence, Oded's idea.
123      do islope = 1,nslope
124        do ig = 1,ngrid
125         if((ice_table_beg(ig,islope).gt.previous_icetable_depth(ig,islope)).and.(previous_icetable_depth(ig,islope).ge.0)) then
126            call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),wice_inertia)
127            stretch = (regolith_inertia(ig,islope)/wice_inertia)**2
128
129            ice_table_thickness(ig,islope) = ice_table_thickness(ig,islope) + (ice_table_beg(ig,islope) - &
130            previous_icetable_depth(ig,islope)+(ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch)
131            ice_table_beg(ig,islope) = previous_icetable_depth(ig,islope)+(ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch
132         endif
133        enddo
134      enddo
135
136!=======================================================================
137      RETURN
138#endif
139      END
140
141
142   SUBROUTINE find_layering_icetable(porefill,psat_soil,psat_surf,pwat_surf,psat_bottom,B,index_IS,depth_filling,index_filling,index_geothermal,depth_geothermal,dz_etadz_rho)
143!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
144!!!
145!!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM
146!!!
147!!! Author: LL
148!!!
149!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
150use comsoil_h_PEM,only: nsoilmx_PEM,mlayer_PEM
151use math_mod, only: deriv1,deriv1_onesided,colint,findroot,deriv2_simple
152implicit none
153! inputs
154! ------
155
156real,intent(in) :: porefill(nsoilmx_PEM) ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice
157real,intent(in) :: psat_soil(nsoilmx_PEM) ! Soil water pressure at saturation, yearly averaged [Pa]
158real,intent(in) :: psat_surf ! surface water pressure at saturation, yearly averaged [Pa]
159real,intent(in) :: pwat_surf ! Water vapor pressure  at the surface, not necesseraly at saturation, yearly averaged [Pa]
160real,intent(in) :: psat_bottom ! Boundary conditions for soil vapor pressure [Pa]
161real,intent(in) :: B ! constant (Eq. 8 from  Schorgofer, Icarus (2010).)
162integer, intent(in) :: index_IS ! index of the soil layer where the ice sheet begins [1]
163real, intent(inout) :: depth_filling ! depth where pore filling begins [m]
164
165! outputs
166! -------
167integer,intent(out) :: index_filling ! index where the pore filling begins [1]
168integer, intent(out) :: index_geothermal ! index where the ice table stops [1]
169real, intent(out) :: depth_geothermal ! depth where the ice table stops [m]
170real, intent(out) :: dz_etadz_rho(nsoilmx_PEM) ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase
171
172! local
173
174real :: eta(nsoilmx_PEM)  ! constriction
175integer :: ilay ! index for loop
176real :: old_depth_filling ! depth_filling saved
177real :: dz_psat(nsoilmx_PEM) ! first derivative of the vapor pressure at saturationn
178integer :: index_tmp ! for loop
179real :: Jdry ! flux trought the dry layer
180real :: Jsat ! flux trought the ice layer
181real :: Jdry_prevlay,Jsat_prevlay ! same but for the previous ice layer
182integer :: index_firstice ! first index where ice appears (i.e., f > 0)
183real :: dz_eta(nsoilmx_PEM) ! \partial z \eta
184real :: dz_eta_low ! same but evaluated at the interface for ice
185real :: dzz_psat(nsoilmx_PEM) ! \partial \partial psat
186real :: massfillabove,massfillafter ! h2O mass above and after index_geothermal
187
188! constant
189real :: pvap2rho = 18.e-3/8.314
190!
191
192
193
194
195 
196! 0. Compute constriction over the layer
197! Within the ice sheet, constriction is set to 0. Elsewhere, constriction =  (1-porefilling)**2
198if (index_IS.lt.0) then
199   index_tmp = nsoilmx_PEM
200   do ilay = 1,nsoilmx_PEM
201       call constriction(porefill(ilay),eta(ilay))
202   enddo
203else
204   index_tmp = index_IS
205   do ilay = 1,index_IS-1
206       call constriction(porefill(ilay),eta(ilay))
207   enddo
208   do ilay = index_IS,nsoilmx_PEM
209       eta(ilay) = 0.
210   enddo
211endif
212
213! 1. Depth at which pore filling occurs. We solve Eq. 9 from  Schorgofer, Icarus  (2010)   
214
215old_depth_filling = depth_filling
216
217call deriv1(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dz_psat) 
218
219do ilay = 1,index_tmp
220  Jdry = (psat_soil(ilay) - pwat_surf)/mlayer_PEM(ilay) ! left member of Eq. 9 from Schorgofer, Icarus  (2010)
221  Jsat =  eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010)
222  if((Jdry - Jsat).le.0) then
223     index_filling = ilay
224     exit
225   endif
226enddo
227
228if(index_filling.eq.1)  depth_filling = mlayer_PEM(1)
229
230if(index_filling.gt.1) then
231    Jdry_prevlay = (psat_soil(index_filling-1) - pwat_surf)/mlayer_PEM(index_filling-1)
232    Jsat_prevlay = eta(index_filling-1)*dz_psat(index_filling-1)
233   call findroot(Jdry-Jsat,Jdry_prevlay-Jsat_prevlay,mlayer_PEM(index_filling),mlayer_PEM(index_filling-1),depth_filling)
234endif
235
236
237! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010))
238
239
240! 2.0 preliminary: depth to shallowest  ice (discontinuity at interface)
241
242index_firstice = -1
243do ilay = 1,nsoilmx_PEM
244   if (porefill(ilay).le.0.) then
245        index_firstice = ilay  ! first point with ice
246        exit
247    endif
248enddo
249
250! 2.1: now we can computeCompute d_z (eta* d_z(rho))
251
252call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM-1),dz_eta)
253
254if ((index_firstice.gt.0).and.(index_firstice.lt.nsoilmx_PEM-2)) then
255     call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice))
256endif
257
258call deriv2_simple(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dzz_psat)
259dz_etadz_rho(:) = pvap2rho*(dz_eta(:)*dz_psat(:) + eta(:)*dzz_psat(:))
260
261! 3. Ice table boundary due to geothermal heating
262
263if(index_IS.gt.0) index_geothermal = -1
264if(index_geothermal.lt.0) depth_geothermal = -1.
265if((index_geothermal.gt.0).and.(index_IS.lt.0)) then ! Eq. 21 from Schorfoger, Icarus (2010)
266   index_geothermal = -1
267   do ilay=2,nsoilmx_PEM
268        if (dz_psat(ilay).gt.0.) then  ! first point with reversed flux
269           index_geothermal=ilay
270           call findroot(dz_psat(ilay-1),dz_psat(ilay),mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal)
271           exit
272        endif
273     enddo
274  else
275     index_geothermal = -1
276  endif
277 if ((index_geothermal.gt.0).and.(index_IS.lt.0)) then !  Eq. 24 from Schorgofer, Icarus (2010)
278     call  colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,index_geothermal-1,nsoilmx_PEM,massfillabove)
279     index_tmp = -1
280     do ilay=index_geothermal,nsoilmx_PEM
281        if (minval(eta(ilay:nsoilmx_PEM)).le.0.) cycle  ! eta=0 means completely full
282         call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter)
283        if (massfillafter<dz_psat(ilay)*pvap2rho*B) then  ! usually executes on i=typeG
284           if (ilay>index_geothermal) then
285!              write(34,*) '# adjustment to geotherm depth by',ilay-index_geothermal
286              call findroot(dz_psat(ilay-1)*pvap2rho*B-massfillabove, & 
287                   dz_psat(ilay)*pvap2rho*B-massfillafter,mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal)
288!               if (depth_geothermal.gt.mlayer_PEM(ilay) .or. depth_geothermal.lt.<mlayer_PEM(ilay-1)) write(34,*)
289!                     '# WARNING: zdepthG interpolation failed',ilay,mlayer_PEM(ilay-1),depth_geothermal,mlayer_PEM(ilay)
290              index_tmp=ilay
291           endif
292           ! otherwise leave depth_geothermal unchanged
293           exit
294        endif
295        massfillabove = massfillafter
296     enddo
297     if (index_tmp.gt.0) index_geothermal = index_tmp
298  end if
299  return
300end subroutine
301
302
303
304
305SUBROUTINE constriction(porefill,eta)
306
307        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
308        !!!
309        !!! Purpose: Compute the  constriction of vapor flux by pore ice
310        !!!
311        !!! Author: LL
312        !!!
313        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
314        implicit none
315        real,intent(in) :: porefill ! pore filling fraction
316        real,intent(out) :: eta ! constriction
317
318        !!!
319
320
321        if (porefill.le.0.) eta = 1.
322        if ((porefill.gt.0.) .and.(porefill.lt.1.)) then
323        eta = (1-porefill)**2  ! Hudson et al., JGR, 2009
324        endif
325        if (porefill.le.1.) eta = 0.
326        return
327        end subroutine
328
329
330end module
Note: See TracBrowser for help on using the repository browser.