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

Last change on this file since 3032 was 3032, checked in by llange, 15 months ago

MARS PEM
Fix bug with albedo (wrong albedo for surface where co2 ice is) and add the porosity in the computation commited previously (missing)
LL
File size: 17.8 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   SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness)
42!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
43!!!
44!!! Purpose: Compute the ice table depth (pore-filling) knowing the yearly average water
45!!! density at the surface and at depth.
46!!! Computations are made following the methods in Schorgofer et al., 2005
47!!! This subroutine only gives the ice table at equilibrium and does not consider exchange with the atmosphere
48!!!
49!!! Author: LL
50!!!
51!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
52    use math_mod,only: findroot
53    USE comsoil_h_PEM, only: mlayer_PEM,layer_PEM                  ! Depth of the vertical grid
54    USE soil_thermalproperties_mod, only: ice_thermal_properties
55    implicit none
56!   inputs
57! -----------
58    integer,intent(in) :: ngrid,nslope,nsoil_PEM                   ! Size of the physical grid, number of subslope, number of soil layer in the PEM
59    logical,intent(in) :: watercaptag(ngrid)                       ! Boolean to check the presence of a perennial glacier
60    real,intent(in) :: rhowatersurf_ave(ngrid,nslope)              ! Water density at the surface, yearly averaged [kg/m^3]
61    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]
62    real,intent(in) :: regolith_inertia(ngrid,nslope)              ! Thermal inertia of the regolith layer [SI]
63!   Ouputs
64! -----------
65    real,intent(out) :: ice_table_beg(ngrid,nslope)               ! ice table depth [m]
66    real,intent(out) :: ice_table_thickness(ngrid,nslope)         ! ice table thickness [m]
67!   Local
68! -----------
69    real :: z1,z2                                                  ! intermediate variables used when doing a linear interpolation between two depths to find the root
70    integer ig, islope,isoil,isoilend                              ! loop variables
71    real :: diff_rho(nsoil_PEM)                                    ! difference of water vapor density between the surface and at depth [kg/m^3]
72    real  :: ice_table_end                                         ! depth of the end of the ice table  [m]
73    real :: previous_icetable_depth(ngrid,nslope)                  ! Ice table computed at previous ice depth [m]
74    real :: stretch                                                ! stretch factor to improve the convergence of the ice table
75    real :: wice_inertia                                           ! Water Ice thermal Inertia [USI]   
76!   Code
77! -----------
78
79     previous_icetable_depth(:,:) = ice_table_beg(:,:)
80     do ig = 1,ngrid
81      if(watercaptag(ig)) then
82         ice_table_beg(ig,:) = 0.
83         ice_table_thickness(ig,:) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
84      else
85        do islope = 1,nslope
86           ice_table_beg(ig,islope) = -1.
87           ice_table_thickness(ig,islope) = 0.
88         do isoil = 1,nsoil_PEM
89           diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope)
90         enddo
91         if(diff_rho(1) > 0) then                                   ! ice is at the surface
92           ice_table_beg(ig,islope) = 0.
93               do isoilend = 2,nsoil_PEM-1
94                 if((diff_rho(isoilend).gt.0).and.(diff_rho(isoilend+1).lt.0.)) then
95                  call findroot(diff_rho(isoilend),diff_rho(isoilend+1),mlayer_PEM(isoilend),mlayer_PEM(isoilend+1),ice_table_end)
96                  ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope)
97                  exit
98                 endif
99               enddo
100         else
101           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
102             if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then
103               call findroot(diff_rho(isoil),diff_rho(isoil+1),mlayer_PEM(isoil),mlayer_PEM(isoil+1),ice_table_beg(ig,islope))
104               ! Now let's find the end of the ice table
105               ice_table_thickness(ig,islope) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
106               do isoilend = isoil+1,nsoil_PEM-1
107                 if((diff_rho(isoilend).gt.0).and.(diff_rho(isoilend+1).lt.0.)) then
108                  call findroot(diff_rho(isoilend),diff_rho(isoilend+1),mlayer_PEM(isoilend),mlayer_PEM(isoilend+1),ice_table_end)
109                  ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope)
110                  exit
111                 endif
112               enddo
113             endif !diff_rho(z) <0 & diff_rho(z+1) > 0
114           enddo
115          endif ! diff_rho(1) > 0
116        enddo
117       endif ! watercaptag
118      enddo
119
120! Small trick to speed up the convergence, Oded's idea.
121      do islope = 1,nslope
122        do ig = 1,ngrid
123         if((ice_table_beg(ig,islope).gt.previous_icetable_depth(ig,islope)).and.(previous_icetable_depth(ig,islope).ge.0)) then
124            call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),wice_inertia)
125            stretch = (regolith_inertia(ig,islope)/wice_inertia)**2
126
127            ice_table_thickness(ig,islope) = ice_table_thickness(ig,islope) + (ice_table_beg(ig,islope) - &
128            previous_icetable_depth(ig,islope)+(ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch)
129            ice_table_beg(ig,islope) = previous_icetable_depth(ig,islope)+(ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch
130         endif
131        enddo
132      enddo
133
134      RETURN
135      END
136
137!!! --------------------------------------
138
139
140   SUBROUTINE compute_massh2o_exchange_ssi(ngrid,nslope,nsoil_PEM,former_ice_table_thickness,new_ice_table_thickness,ice_table_depth,tsoil,delta_m_h2o)
141!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
142!!!
143!!! Purpose: Compute the mass of H2O that has sublimated from the ice table / condensed
144!!!
145!!! Author: LL
146!!!
147!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
148    use comsoil_h_PEM, only: mlayer_PEM   
149    use comslope_mod, only: subslope_dist,def_slope_mean           
150    use comconst_mod,only: pi       
151    use constants_marspem_mod,only:porosity
152    implicit none
153!   inputs
154! -----------
155    integer,intent(in) :: ngrid,nslope,nsoil_PEM                   ! Size of the physical grid, number of subslope, number of soil layer in the PEM
156    real,intent(in) :: former_ice_table_thickness(ngrid,nslope)    ! ice table thickness at the former iteration [m]
157    real,intent(in) :: new_ice_table_thickness(ngrid,nslope)       ! ice table thickness at the current iteration [m]
158    real,intent(out) :: ice_table_depth(ngrid,nslope)              ! ice table depth [m]
159    real,intent(in) :: tsoil(ngrid,nsoil_PEM,nslope)               ! Soil temperature [K]
160!   outputs
161! -----------
162    real,intent(out) :: delta_m_h2o(ngrid)                         ! Mass of H2O ice that has been condensed on the ice table / sublimates from the ice table [kg/m^2]
163
164! local
165    real :: rho(ngrid,nslope)                                      ! density of water ice [kg/m^3]
166    integer :: ig,islope,ilay,iref                                 ! loop index 
167    real :: Tice(ngrid,nslope)                                     ! ice temperature [k]
168!   Code
169! -----------
170   rho(:,:) = 0.
171   Tice(:,:) = 0.
172!1. First let's compute Tice using a linear interpolation between the mlayer level
173   do ig = 1,ngrid
174      do islope = 1,nslope
175         if(ice_table_depth(ig,islope).gt.0.) then
176            if (ice_table_depth(ig,islope).lt. 1e-10) then
177               Tice(ig,islope) = tsoil(ig,1,islope)
178            else
179               iref=0 ! initialize iref
180               do ilay=1,nsoil_PEM ! loop on layers to find the beginning of the ice table
181                  if (ice_table_depth(ig,islope).ge.mlayer_PEM(ilay)) then
182                     iref=ilay ! pure regolith layer up to here
183                  else
184                  ! correct iref was obtained in previous cycle
185                     exit
186                  endif
187               enddo
188               if(iref.eq.nsoil_PEM) then
189                  Tice(ig,islope) = tsoil(ig,nsoil_PEM,islope)
190               else
191                  Tice(ig,islope) = (tsoil(ig,iref+1,islope) - tsoil(ig,iref,islope))/(mlayer_PEM(iref+1) - mlayer_PEM(iref))* &
192                                    (ice_table_depth(ig,islope) - mlayer_PEM(iref)) +  tsoil(ig,iref,islope)           
193               endif         
194               rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+   0.0351* Tice(ig,islope) +  933.5030 ! Rottgers, 2012
195            endif
196         endif
197      enddo
198   enddo     
199
200         
201!2. Let's compute the amount of ice that has sublimated in each subslope
202   do ig = 1,ngrid
203      do islope = 1,nslope
204        delta_m_h2o(ig)  = delta_m_h2o(ig) +  porosity*rho(ig,islope)*(new_ice_table_thickness(ig,islope) - former_ice_table_thickness(ig,islope)) & ! convention > 0. <=> it condenses 
205                           *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
206      enddo
207   enddo
208
209   return
210end subroutine
211
212!!! --------------------------------------
213
214   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)
215!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
216!!!
217!!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM
218!!!
219!!! Author: LL
220!!!
221!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
222use comsoil_h_PEM,only: nsoilmx_PEM,mlayer_PEM
223use math_mod, only: deriv1,deriv1_onesided,colint,findroot,deriv2_simple
224implicit none
225! inputs
226! ------
227
228real,intent(in) :: porefill(nsoilmx_PEM) ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice
229real,intent(in) :: psat_soil(nsoilmx_PEM) ! Soil water pressure at saturation, yearly averaged [Pa]
230real,intent(in) :: psat_surf ! surface water pressure at saturation, yearly averaged [Pa]
231real,intent(in) :: pwat_surf ! Water vapor pressure  at the surface, not necesseraly at saturation, yearly averaged [Pa]
232real,intent(in) :: psat_bottom ! Boundary conditions for soil vapor pressure [Pa]
233real,intent(in) :: B ! constant (Eq. 8 from  Schorgofer, Icarus (2010).)
234integer, intent(in) :: index_IS ! index of the soil layer where the ice sheet begins [1]
235real, intent(inout) :: depth_filling ! depth where pore filling begins [m]
236
237! outputs
238! -------
239integer,intent(out) :: index_filling ! index where the pore filling begins [1]
240integer, intent(out) :: index_geothermal ! index where the ice table stops [1]
241real, intent(out) :: depth_geothermal ! depth where the ice table stops [m]
242real, intent(out) :: dz_etadz_rho(nsoilmx_PEM) ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase
243
244! local
245
246real :: eta(nsoilmx_PEM)  ! constriction
247integer :: ilay ! index for loop
248real :: old_depth_filling ! depth_filling saved
249real :: dz_psat(nsoilmx_PEM) ! first derivative of the vapor pressure at saturationn
250integer :: index_tmp ! for loop
251real :: Jdry ! flux trought the dry layer
252real :: Jsat ! flux trought the ice layer
253real :: Jdry_prevlay,Jsat_prevlay ! same but for the previous ice layer
254integer :: index_firstice ! first index where ice appears (i.e., f > 0)
255real :: dz_eta(nsoilmx_PEM) ! \partial z \eta
256real :: dz_eta_low ! same but evaluated at the interface for ice
257real :: dzz_psat(nsoilmx_PEM) ! \partial \partial psat
258real :: massfillabove,massfillafter ! h2O mass above and after index_geothermal
259
260! constant
261real :: pvap2rho = 18.e-3/8.314
262!
263 
264! 0. Compute constriction over the layer
265! Within the ice sheet, constriction is set to 0. Elsewhere, constriction =  (1-porefilling)**2
266if (index_IS.lt.0) then
267   index_tmp = nsoilmx_PEM
268   do ilay = 1,nsoilmx_PEM
269       call constriction(porefill(ilay),eta(ilay))
270   enddo
271else
272   index_tmp = index_IS
273   do ilay = 1,index_IS-1
274       call constriction(porefill(ilay),eta(ilay))
275   enddo
276   do ilay = index_IS,nsoilmx_PEM
277       eta(ilay) = 0.
278   enddo
279endif
280
281! 1. Depth at which pore filling occurs. We solve Eq. 9 from  Schorgofer, Icarus  (2010)   
282
283old_depth_filling = depth_filling
284
285call deriv1(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dz_psat) 
286
287do ilay = 1,index_tmp
288  Jdry = (psat_soil(ilay) - pwat_surf)/mlayer_PEM(ilay) ! left member of Eq. 9 from Schorgofer, Icarus  (2010)
289  Jsat =  eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010)
290  if((Jdry - Jsat).le.0) then
291     index_filling = ilay
292     exit
293   endif
294enddo
295
296if(index_filling.eq.1)  depth_filling = mlayer_PEM(1)
297
298if(index_filling.gt.1) then
299    Jdry_prevlay = (psat_soil(index_filling-1) - pwat_surf)/mlayer_PEM(index_filling-1)
300    Jsat_prevlay = eta(index_filling-1)*dz_psat(index_filling-1)
301   call findroot(Jdry-Jsat,Jdry_prevlay-Jsat_prevlay,mlayer_PEM(index_filling),mlayer_PEM(index_filling-1),depth_filling)
302endif
303
304
305! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010))
306
307
308! 2.0 preliminary: depth to shallowest  ice (discontinuity at interface)
309
310index_firstice = -1
311do ilay = 1,nsoilmx_PEM
312   if (porefill(ilay).le.0.) then
313        index_firstice = ilay  ! first point with ice
314        exit
315    endif
316enddo
317
318! 2.1: now we can computeCompute d_z (eta* d_z(rho))
319
320call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM-1),dz_eta)
321
322if ((index_firstice.gt.0).and.(index_firstice.lt.nsoilmx_PEM-2)) then
323     call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice))
324endif
325
326call deriv2_simple(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dzz_psat)
327dz_etadz_rho(:) = pvap2rho*(dz_eta(:)*dz_psat(:) + eta(:)*dzz_psat(:))
328
329! 3. Ice table boundary due to geothermal heating
330
331if(index_IS.gt.0) index_geothermal = -1
332if(index_geothermal.lt.0) depth_geothermal = -1.
333if((index_geothermal.gt.0).and.(index_IS.lt.0)) then ! Eq. 21 from Schorfoger, Icarus (2010)
334   index_geothermal = -1
335   do ilay=2,nsoilmx_PEM
336        if (dz_psat(ilay).gt.0.) then  ! first point with reversed flux
337           index_geothermal=ilay
338           call findroot(dz_psat(ilay-1),dz_psat(ilay),mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal)
339           exit
340        endif
341     enddo
342  else
343     index_geothermal = -1
344  endif
345 if ((index_geothermal.gt.0).and.(index_IS.lt.0)) then !  Eq. 24 from Schorgofer, Icarus (2010)
346     call  colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,index_geothermal-1,nsoilmx_PEM,massfillabove)
347     index_tmp = -1
348     do ilay=index_geothermal,nsoilmx_PEM
349        if (minval(eta(ilay:nsoilmx_PEM)).le.0.) cycle  ! eta=0 means completely full
350         call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter)
351        if (massfillafter<dz_psat(ilay)*pvap2rho*B) then  ! usually executes on i=typeG
352           if (ilay>index_geothermal) then
353!              write(34,*) '# adjustment to geotherm depth by',ilay-index_geothermal
354              call findroot(dz_psat(ilay-1)*pvap2rho*B-massfillabove, & 
355                   dz_psat(ilay)*pvap2rho*B-massfillafter,mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal)
356!               if (depth_geothermal.gt.mlayer_PEM(ilay) .or. depth_geothermal.lt.<mlayer_PEM(ilay-1)) write(34,*)
357!                     '# WARNING: zdepthG interpolation failed',ilay,mlayer_PEM(ilay-1),depth_geothermal,mlayer_PEM(ilay)
358              index_tmp=ilay
359           endif
360           ! otherwise leave depth_geothermal unchanged
361           exit
362        endif
363        massfillabove = massfillafter
364     enddo
365     if (index_tmp.gt.0) index_geothermal = index_tmp
366  end if
367  return
368end subroutine
369
370!!! --------------------------------------
371
372SUBROUTINE constriction(porefill,eta)
373
374        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
375        !!!
376        !!! Purpose: Compute the  constriction of vapor flux by pore ice
377        !!!
378        !!! Author: LL
379        !!!
380        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
381        implicit none
382        real,intent(in) :: porefill ! pore filling fraction
383        real,intent(out) :: eta ! constriction
384
385        !!!
386
387        if (porefill.le.0.) eta = 1.
388        if ((porefill.gt.0.) .and.(porefill.lt.1.)) then
389        eta = (1-porefill)**2  ! Hudson et al., JGR, 2009
390        endif
391        if (porefill.le.1.) eta = 0.
392        return
393        end subroutine
394
395end module
Note: See TracBrowser for help on using the repository browser.