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

Last change on this file since 2937 was 2902, checked in by llange, 3 years ago

PEM
First implementation of a dynamic ice table following Schorgofer, Icarus (2010) / MSIM
The dynamic of the ice is implemented here. Impact on the soil properties (thermal, inertia) is currently implemented and will be commit later
Small fix for the ice table at equilibrium
LL
File size: 14.4 KB
Line 
1module ice_table_mod
2
3  implicit none
4
5  contains
6
7
8!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
9!!!
10!!! Purpose: Compute the ice table in two ways: dynamic and at equilibrium
11!!! Author:  LL, 02/2023
12!!!
13!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
14
15
16
17   SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,ice_table)
18!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
19!!!
20!!! Purpose: Compute the ice table depth knowing the yearly average water
21!!! density at the surface and at depth.
22!!! Computations are made following the methods in Schorgofer et al., 2005
23!!! This subroutine only gives the ice table at equilibrium and does not consider exchange with the atmosphere
24!!!
25!!! Author: LL
26!!!
27!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
28
29#ifndef CPP_STD
30    USE comsoil_h_PEM, only: mlayer_PEM                             ! Depth of the vertical grid
31    implicit none
32
33    integer,intent(in) :: ngrid,nslope,nsoil_PEM                   ! Size of the physical grid, number of subslope, number of soil layer in the PEM
34    logical,intent(in) :: watercaptag(ngrid)                       ! Boolean to check the presence of a perennial glacier
35    real,intent(in) :: rhowatersurf_ave(ngrid,nslope)              ! Water density at the surface, yearly averaged [kg/m^3]
36    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]
37    real,intent(inout) :: ice_table(ngrid,nslope)                  ! ice table depth [m]
38    real :: z1,z2                                                  ! intermediate variables used when doing a linear interpolation between two depths to find the root
39    integer ig, islope,isoil                                       ! loop variables
40    real :: diff_rho(nsoil_PEM)                                    ! difference of water vapor density between the surface and at depth [kg/m^3]
41
42
43     do ig = 1,ngrid
44      if(watercaptag(ig)) then
45         ice_table(ig,:) = 0.
46      else
47        do islope = 1,nslope
48           ice_table(ig,islope) = -1.
49
50         do isoil = 1,nsoil_PEM
51           diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope)
52
53         enddo
54         if(diff_rho(1) > 0) then                                   ! ice is at the surface
55           ice_table(ig,islope) = 0.
56         else
57           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
58             if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then
59               call findroot(diff_rho(isoil),diff_rho(isoil+1),mlayer_PEM(isoil),mlayer_PEM(isoil+1),ice_table(ig,islope))
60               exit
61             endif !diff_rho(z) <0 & diff_rho(z+1) > 0
62           enddo
63          endif ! diff_rho(1) > 0
64        enddo
65       endif ! watercaptag
66      enddo
67!=======================================================================
68      RETURN
69#endif
70      END
71
72   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)
73!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
74!!!
75!!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM
76!!!
77!!! Author: LL
78!!!
79!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
80use comsoil_h_PEM,only: nsoilmx_PEM,mlayer_PEM
81implicit none
82! inputs
83! ------
84
85real,intent(in) :: porefill(nsoilmx_PEM) ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice
86real,intent(in) :: psat_soil(nsoilmx_PEM) ! Soil water pressure at saturation, yearly averaged [Pa]
87real,intent(in) :: psat_surf ! surface water pressure at saturation, yearly averaged [Pa]
88real,intent(in) :: pwat_surf ! Water vapor pressure  at the surface, not necesseraly at saturation, yearly averaged [Pa]
89real,intent(in) :: psat_bottom ! Boundary conditions for soil vapor pressure [Pa]
90real,intent(in) :: B ! constant (Eq. 8 from  Schorgofer, Icarus (2010).)
91integer, intent(in) :: index_IS ! index of the soil layer where the ice sheet begins [1]
92real, intent(inout) :: depth_filling ! depth where pore filling begins [m]
93
94! outputs
95! -------
96integer,intent(out) :: index_filling ! index where the pore filling begins [1]
97integer, intent(out) :: index_geothermal ! index where the ice table stops [1]
98real, intent(out) :: depth_geothermal ! depth where the ice table stops [m]
99real, intent(out) :: dz_etadz_rho(nsoilmx_PEM) ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase
100
101! local
102
103real :: eta(nsoilmx_PEM)  ! constriction
104integer :: ilay ! index for loop
105real :: old_depth_filling ! depth_filling saved
106real :: dz_psat(nsoilmx_PEM) ! first derivative of the vapor pressure at saturationn
107integer :: index_tmp ! for loop
108real :: Jdry ! flux trought the dry layer
109real :: Jsat ! flux trought the ice layer
110real :: Jdry_prevlay,Jsat_prevlay ! same but for the previous ice layer
111integer :: index_firstice ! first index where ice appears (i.e., f > 0)
112real :: dz_eta(nsoilmx_PEM) ! \partial z \eta
113real :: dz_eta_low ! same but evaluated at the interface for ice
114real :: dzz_psat(nsoilmx_PEM) ! \partial \partial psat
115real :: massfillabove,massfillafter ! h2O mass above and after index_geothermal
116
117! constant
118real :: pvap2rho = 18.e-3/8.314
119!
120
121
122
123
124 
125! 0. Compute constriction over the layer
126! Within the ice sheet, constriction is set to 0. Elsewhere, constriction =  (1-porefilling)**2
127if (index_IS.lt.0) then
128   index_tmp = nsoilmx_PEM
129   do ilay = 1,nsoilmx_PEM
130       call constriction(porefill(ilay),eta(ilay))
131   enddo
132else
133   index_tmp = index_IS
134   do ilay = 1,index_IS-1
135       call constriction(porefill(ilay),eta(ilay))
136   enddo
137   do ilay = index_IS,nsoilmx_PEM
138       eta(ilay) = 0.
139   enddo
140endif
141
142! 1. Depth at which pore filling occurs. We solve Eq. 9 from  Schorgofer, Icarus  (2010)   
143
144old_depth_filling = depth_filling
145
146call deriv1(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dz_psat) 
147
148do ilay = 1,index_tmp
149  Jdry = (psat_soil(ilay) - pwat_surf)/mlayer_PEM(ilay) ! left member of Eq. 9 from Schorgofer, Icarus  (2010)
150  Jsat =  eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010)
151  if((Jdry - Jsat).le.0) then
152     index_filling = ilay
153     exit
154   endif
155enddo
156
157if(index_filling.eq.1)  depth_filling = mlayer_PEM(1)
158
159if(index_filling.gt.1) then
160    Jdry_prevlay = (psat_soil(index_filling-1) - pwat_surf)/mlayer_PEM(index_filling-1)
161    Jsat_prevlay = eta(index_filling-1)*dz_psat(index_filling-1)
162   call findroot(Jdry-Jsat,Jdry_prevlay-Jsat_prevlay,mlayer_PEM(index_filling),mlayer_PEM(index_filling-1),depth_filling)
163endif
164
165
166! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010))
167
168
169! 2.0 preliminary: depth to shallowest  ice (discontinuity at interface)
170
171index_firstice = -1
172do ilay = 1,nsoilmx_PEM
173   if (porefill(ilay).le.0.) then
174        index_firstice = ilay  ! first point with ice
175        exit
176    endif
177enddo
178
179! 2.1: now we can computeCompute d_z (eta* d_z(rho))
180
181call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM-1),dz_eta)
182
183if ((index_firstice.gt.0).and.(index_firstice.lt.nsoilmx_PEM-2)) then
184     call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice))
185endif
186
187call deriv2_simple(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dzz_psat)
188dz_etadz_rho(:) = pvap2rho*(dz_eta(:)*dz_psat(:) + eta(:)*dzz_psat(:))
189
190! 3. Ice table boundary due to geothermal heating
191
192if(index_IS.gt.0) index_geothermal = -1
193if(index_geothermal.lt.0) depth_geothermal = -1.
194if((index_geothermal.gt.0).and.(index_IS.lt.0)) then ! Eq. 21 from Schorfoger, Icarus (2010)
195   index_geothermal = -1
196   do ilay=2,nsoilmx_PEM
197        if (dz_psat(ilay).gt.0.) then  ! first point with reversed flux
198           index_geothermal=ilay
199           call findroot(dz_psat(ilay-1),dz_psat(ilay),mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal)
200           exit
201        endif
202     enddo
203  else
204     index_geothermal = -1
205  endif
206 if ((index_geothermal.gt.0).and.(index_IS.lt.0)) then !  Eq. 24 from Schorgofer, Icarus (2010)
207     call  colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,index_geothermal-1,nsoilmx_PEM,massfillabove)
208     index_tmp = -1
209     do ilay=index_geothermal,nsoilmx_PEM
210        if (minval(eta(ilay:nsoilmx_PEM)).le.0.) cycle  ! eta=0 means completely full
211         call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter)
212        if (massfillafter<dz_psat(ilay)*pvap2rho*B) then  ! usually executes on i=typeG
213           if (ilay>index_geothermal) then
214!              write(34,*) '# adjustment to geotherm depth by',ilay-index_geothermal
215              call findroot(dz_psat(ilay-1)*pvap2rho*B-massfillabove, & 
216                   dz_psat(ilay)*pvap2rho*B-massfillafter,mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal)
217!               if (depth_geothermal.gt.mlayer_PEM(ilay) .or. depth_geothermal.lt.<mlayer_PEM(ilay-1)) write(34,*)
218!                     '# WARNING: zdepthG interpolation failed',ilay,mlayer_PEM(ilay-1),depth_geothermal,mlayer_PEM(ilay)
219              index_tmp=ilay
220           endif
221           ! otherwise leave depth_geothermal unchanged
222           exit
223        endif
224        massfillabove = massfillafter
225     enddo
226     if (index_tmp.gt.0) index_geothermal = index_tmp
227  end if
228  return
229end subroutine
230
231
232
233
234
235SUBROUTINE findroot(y1,y2,z1,z2,zr)
236        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
237        !!!
238        !!! Purpose: Compute the root zr, between two values y1 and y2 at depth z1,z2
239        !!!
240        !!! Author: LL
241        !!!
242        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
243
244        implicit none
245        real,intent(in) :: y1,y2 ! difference between surface water density and at depth  [kg/m^3]
246        real,intent(in) :: z1,z2 ! depth [m]
247        real,intent(out) :: zr ! depth at which we have zero
248        zr = (y1*z2 - y2*z1)/(y1-y2)
249        RETURN
250        end
251
252SUBROUTINE constriction(porefill,eta)
253
254        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
255        !!!
256        !!! Purpose: Compute the  constriction of vapor flux by pore ice
257        !!!
258        !!! Author: LL
259        !!!
260        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
261        implicit none
262        real,intent(in) :: porefill ! pore filling fraction
263        real,intent(out) :: eta ! constriction
264
265        !!!
266
267
268        if (porefill.le.0.) eta = 1.
269        if ((porefill.gt.0.) .and.(porefill.lt.1.)) then
270        eta = (1-porefill)**2  ! Hudson et al., JGR, 2009
271        endif
272        if (porefill.le.1.) eta = 0.
273        return
274        end
275
276
277
278subroutine deriv1(z,nz,y,y0,ybot,dzY)
279!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
280!!!
281!!! Purpose: Compute the first derivative of a function y(z) on an irregular grid
282!!!           
283!!! Author: From N.S (N.S, Icarus 2010), impletented here by LL
284!!!
285!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
286        implicit none
287
288! first derivative of a function y(z) on irregular grid
289! upper boundary conditions: y(0)=y0
290! lower boundary condition.: yp = ybottom
291  integer, intent(IN) :: nz ! number of layer
292  real, intent(IN) :: z(nz) ! depth layer
293  real, intent(IN) :: y(nz) ! function which needs to be derived
294  real, intent(IN) :: y0,ybot ! boundary conditions
295  real, intent(OUT) :: dzY(nz) ! derivative of y w.r.t depth
296! local
297  integer :: j
298  real :: hm,hp,c1,c2,c3
299
300  hp = z(2)-z(1)
301  c1 = z(1)/(hp*z(2))
302  c2 = 1/z(1) - 1/(z(2)-z(1))
303  c3 = -hp/(z(1)*z(2))
304  dzY(1) = c1*y(2) + c2*y(1) + c3*y0
305  do j=2,nz-1
306     hp = z(j+1)-z(j)
307     hm = z(j)-z(j-1)
308     c1 = +hm/(hp*(z(j+1)-z(j-1)))
309     c2 = 1/hm - 1/hp
310     c3 = -hp/(hm*(z(j+1)-z(j-1)))
311     dzY(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1)
312  enddo
313  dzY(nz) = (ybot - y(nz-1))/(2.*(z(nz)-z(nz-1)))
314 return
315end subroutine deriv1
316
317
318
319subroutine deriv2_simple(z,nz,y,y0,yNp1,yp2)
320!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
321!!!
322!!! Purpose: Compute the second derivative of a function y(z) on an irregular grid
323!!!           
324!!! Author: N.S (raw copy/paste from MSIM)
325!!!
326!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
327
328  ! second derivative y_zz on irregular grid
329  ! boundary conditions: y(0)=y0, y(nz+1)=yNp1
330  implicit none
331  integer, intent(IN) :: nz
332  real, intent(IN) :: z(nz),y(nz),y0,yNp1
333  real, intent(OUT) :: yp2(nz)
334  integer j
335  real hm,hp,c1,c2,c3
336
337  c1 = +2./((z(2)-z(1))*z(2))
338  c2 = -2./((z(2)-z(1))*z(1))
339  c3 = +2./(z(1)*z(2))
340  yp2(1) = c1*y(2) + c2*y(1) + c3*y0
341
342  do j=2,nz-1
343     hp = z(j+1)-z(j)
344     hm = z(j)-z(j-1)
345     c1 = +2./(hp*(z(j+1)-z(j-1)))
346     c2 = -2./(hp*hm)
347     c3 = +2./(hm*(z(j+1)-z(j-1)))
348     yp2(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1)
349  enddo
350  yp2(nz) = (yNp1 - 2*y(nz) + y(nz-1))/(z(nz)-z(nz-1))**2
351  return
352end subroutine deriv2_simple
353
354
355
356subroutine  deriv1_onesided(j,z,nz,y,dy_zj)
357!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
358!!!
359!!! Purpose:   First derivative of function y(z) at z(j)  one-sided derivative on irregular grid
360!!!           
361!!! Author: N.S (raw copy/paste from MSIM)
362!!!
363!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
364
365  implicit none
366  integer, intent(IN) :: nz,j
367  real, intent(IN) :: z(nz),y(nz)
368  real, intent(out) :: dy_zj
369  real h1,h2,c1,c2,c3
370  if (j<1 .or. j>nz-2) then
371     dy_zj = -1.
372  else
373     h1 = z(j+1)-z(j)
374     h2 = z(j+2)-z(j+1)
375     c1 = -(2*h1+h2)/(h1*(h1+h2))
376     c2 =  (h1+h2)/(h1*h2)
377     c3 = -h1/(h2*(h1+h2))
378     dy_zj = c1*y(j) + c2*y(j+1) + c3*y(j+2)
379  endif
380  return
381end subroutine deriv1_onesided
382
383
384subroutine  colint(y,z,nz,i1,i2,integral)
385!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
386!!!
387!!! Purpose:  Column integrates y on irregular grid
388!!!           
389!!! Author: N.S (raw copy/paste from MSIM)
390!!!
391!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
392 
393  implicit none
394  integer, intent(IN) :: nz, i1, i2
395  real, intent(IN) :: y(nz), z(nz)
396  real,intent(out) :: integral
397  integer i
398  real dz(nz)
399 
400  dz(1) = (z(2)-0.)/2
401  do i=2,nz-1
402     dz(i) = (z(i+1)-z(i-1))/2.
403  enddo
404  dz(nz) = z(nz)-z(nz-1)
405  integral = sum(y(i1:i2)*dz(i1:i2))
406end subroutine colint
407
408
409
410
411
412
413        end module
Note: See TracBrowser for help on using the repository browser.