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

Last change on this file since 3383 was 3328, checked in by llange, 6 months ago

PEM

  • Remove a wrong line which modified the initial surface sublimating during

the main time loop

  • Cleaning the previous cleaning commit (rage pas JB)

LL
File size: 19.9 KB
Line 
1MODULE ice_table_mod
2
3implicit 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
12logical, save                           :: icetable_equilibrium     ! Boolean to say if the PEM needs to recompute the icetable depth when at  equilibrium
13logical, save                           :: icetable_dynamic         ! Boolean to say if the PEM needs to recompute the icetable depth (dynamic method)
14real, allocatable, dimension(:,:), save :: porefillingice_depth     ! ngrid x nslope: Depth of the ice table [m]
15real, allocatable, dimension(:,:), save :: porefillingice_thickness ! ngrid x nslope: Thickness of the ice table [m]
16
17!-----------------------------------------------------------------------
18contains
19!-----------------------------------------------------------------------
20SUBROUTINE ini_ice_table_porefilling(ngrid,nslope)
21
22implicit none
23
24integer, intent(in) :: ngrid  ! number of atmospheric columns
25integer, intent(in) :: nslope ! number of slope within a mesh
26
27allocate(porefillingice_depth(ngrid,nslope))
28allocate(porefillingice_thickness(ngrid,nslope))
29
30END SUBROUTINE ini_ice_table_porefilling
31
32!-----------------------------------------------------------------------
33SUBROUTINE end_ice_table_porefilling
34
35implicit none
36
37if (allocated(porefillingice_depth)) deallocate(porefillingice_depth)
38if (allocated(porefillingice_thickness)) deallocate(porefillingice_thickness)
39
40END SUBROUTINE end_ice_table_porefilling
41
42!------------------------------------------------------------------------------------------------------
43SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness)
44!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
45!!!
46!!! Purpose: Compute the ice table depth (pore-filling) knowing the yearly average water
47!!! density at the surface and at depth.
48!!! Computations are made following the methods in Schorgofer et al., 2005
49!!! This SUBROUTINE only gives the ice table at equilibrium and does not consider exchange with the atmosphere
50!!!
51!!! Author: LL
52!!!
53!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
54use math_mod,                   only: findroot
55use comsoil_h_PEM,              only: mlayer_PEM, layer_PEM ! Depth of the vertical grid
56use soil_thermalproperties_mod, only: ice_thermal_properties
57
58implicit none
59
60! Inputs
61! ------
62integer,                                 intent(in) :: ngrid, nslope, nsoil_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM
63logical, dimension(ngrid),               intent(in) :: watercaptag              ! Boolean to check the presence of a perennial glacier
64real, dimension(ngrid,nslope),           intent(in) :: rhowatersurf_ave         ! Water density at the surface, yearly averaged [kg/m^3]
65real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: rhowatersoil_ave         ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged  [kg/m^3]
66real, dimension(ngrid,nslope),           intent(in) :: regolith_inertia         ! Thermal inertia of the regolith layer [SI]
67! Ouputs
68! ------
69real, dimension(ngrid,nslope), intent(out) :: ice_table_beg       ! ice table depth [m]
70real, dimension(ngrid,nslope), intent(out) :: ice_table_thickness ! ice table thickness [m]
71! Locals
72! ------
73integer                       :: ig, islope, isoil, isoilend ! loop variables
74real, dimension(nsoil_PEM)    :: diff_rho                    ! difference of water vapor density between the surface and at depth [kg/m^3]
75real                          :: ice_table_end               ! depth of the end of the ice table  [m]
76real, dimension(ngrid,nslope) :: previous_icetable_depth     ! Ice table computed at previous ice depth [m]
77real                          :: stretch                     ! stretch factor to improve the convergence of the ice table
78real                          :: wice_inertia                ! Water Ice thermal Inertia [USI]
79! Code
80! ----
81previous_icetable_depth = ice_table_beg
82do 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) > 0 .and. diff_rho(isoilend + 1) < 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) < 0 .and. diff_rho(isoil + 1) > 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) > 0 .and. diff_rho(isoilend + 1) < 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
120enddo
121
122! Small trick to speed up the convergence, Oded's idea.
123do islope = 1,nslope
124    do ig = 1,ngrid
125        if (ice_table_beg(ig,islope) > previous_icetable_depth(ig,islope) .and. previous_icetable_depth(ig,islope) >= 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            ice_table_thickness(ig,islope) = ice_table_thickness(ig,islope) + (ice_table_beg(ig,islope) - &
129                                             previous_icetable_depth(ig,islope) + (ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch)
130            ice_table_beg(ig,islope) = previous_icetable_depth(ig,islope) + (ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch
131        endif
132    enddo
133enddo
134
135END SUBROUTINE computeice_table_equilibrium
136
137!-----------------------------------------------------------------------
138SUBROUTINE compute_massh2o_exchange_ssi(ngrid,nslope,nsoil_PEM,former_ice_table_thickness,new_ice_table_thickness,ice_table_depth,tsurf,tsoil,delta_m_h2o)
139!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
140!!!
141!!! Purpose: Compute the mass of H2O that has sublimated from the ice table / condensed
142!!!
143!!! Author: LL
144!!!
145!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
146use comsoil_h_PEM,         only: mlayer_PEM
147use comslope_mod,          only: subslope_dist, def_slope_mean
148use constants_marspem_mod, only: porosity
149#ifndef CPP_STD
150    use comcstfi_h,   only: pi
151#else
152    use comcstfi_mod, only: pi
153#endif
154
155implicit none
156
157! Inputs
158! ------
159integer,                                 intent(in) :: ngrid, nslope, nsoil_PEM   ! Size of the physical grid, number of subslope, number of soil layer in the PEM
160real, dimension(ngrid,nslope),           intent(in) :: former_ice_table_thickness ! ice table thickness at the former iteration [m]
161real, dimension(ngrid,nslope),           intent(in) :: new_ice_table_thickness    ! ice table thickness at the current iteration [m]
162real, dimension(ngrid,nslope),           intent(in) :: ice_table_depth            ! ice table depth [m]
163real, dimension(ngrid,nslope),           intent(in) :: tsurf                      ! Surface temperature [K]
164real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil                      ! Soil temperature [K]
165! Outputs
166! -------
167real, dimension(ngrid), intent(out) :: delta_m_h2o ! Mass of H2O ice that has been condensed on the ice table / sublimates from the ice table [kg/m^2]
168! Locals
169!-------
170integer                       :: ig, islope, ilay, iref ! loop index
171real, dimension(ngrid,nslope) :: rho                    ! density of water ice [kg/m^3]
172real, dimension(ngrid,nslope) :: Tice                   ! ice temperature [k]
173! Code
174! ----
175rho = 0.
176Tice = 0.
177!1. First let's compute Tice using a linear interpolation between the mlayer level
178do ig = 1,ngrid
179    do islope = 1,nslope
180        call compute_Tice_pem(nsoil_PEM,tsoil(ig,:,islope),tsurf(ig,islope),ice_table_depth(ig,islope),Tice(ig,islope))
181        rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2 + 0.0351* Tice(ig,islope) + 933.5030 ! Rottgers, 2012
182    enddo
183enddo
184
185!2. Let's compute the amount of ice that has sublimated in each subslope
186do ig = 1,ngrid
187    do islope = 1,nslope
188        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
189                          *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
190    enddo
191enddo
192
193END SUBROUTINE compute_massh2o_exchange_ssi
194
195!-----------------------------------------------------------------------
196SUBROUTINE 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)
197!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
198!!!
199!!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM
200!!!
201!!! Author: LL
202!!!
203!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
204use comsoil_h_PEM, only: nsoilmx_PEM, mlayer_PEM
205use math_mod,      only: deriv1, deriv1_onesided, colint, findroot, deriv2_simple
206
207implicit none
208
209! Inputs
210! ------
211real, dimension(nsoilmx_PEM), intent(in) :: porefill    ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice
212real, dimension(nsoilmx_PEM), intent(in) :: psat_soil   ! Soil water pressure at saturation, yearly averaged [Pa]
213real,                         intent(in) :: psat_surf   ! surface water pressure at saturation, yearly averaged [Pa]
214real,                         intent(in) :: pwat_surf   ! Water vapor pressure  at the surface, not necesseraly at saturation, yearly averaged [Pa]
215real,                         intent(in) :: psat_bottom ! Boundary conditions for soil vapor pressure [Pa]
216real,                         intent(in) :: B           ! constant (Eq. 8 from  Schorgofer, Icarus (2010).)
217integer,                      intent(in) :: index_IS    ! index of the soil layer where the ice sheet begins [1]
218real, intent(inout) :: depth_filling ! depth where pore filling begins [m]
219! Outputs
220! -------
221integer,                      intent(out) :: index_filling    ! index where the pore filling begins [1]
222integer,                      intent(out) :: index_geothermal ! index where the ice table stops [1]
223real,                         intent(out) :: depth_geothermal ! depth where the ice table stops [m]
224real, dimension(nsoilmx_PEM), intent(out) :: dz_etadz_rho     ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase
225! Locals
226!-------
227real, dimension(nsoilmx_PEM) :: eta                          ! constriction
228real, dimension(nsoilmx_PEM) :: dz_psat                      ! first derivative of the vapor pressure at saturation
229real, dimension(nsoilmx_PEM) :: dz_eta                       ! \partial z \eta
230real, dimension(nsoilmx_PEM) :: dzz_psat                     ! \partial \partial psat
231integer                      :: ilay, index_tmp              ! index for loop
232real                         :: old_depth_filling            ! depth_filling saved
233real                         :: Jdry                         ! flux trought the dry layer
234real                         :: Jsat                         ! flux trought the ice layer
235real                         :: Jdry_prevlay, Jsat_prevlay   ! same but for the previous ice layer
236integer                      :: index_firstice               ! first index where ice appears (i.e., f > 0)
237real                         :: massfillabove, massfillafter ! h2O mass above and after index_geothermal
238! Constant
239!---------
240real, parameter :: pvap2rho = 18.e-3/8.314
241! Code
242!-----
243! 0. Compute constriction over the layer
244! Within the ice sheet, constriction is set to 0. Elsewhere, constriction =  (1-porefilling)**2
245if (index_IS < 0) then
246    index_tmp = nsoilmx_PEM
247    do ilay = 1,nsoilmx_PEM
248        call constriction(porefill(ilay),eta(ilay))
249    enddo
250else
251    index_tmp = index_IS
252    do ilay = 1,index_IS - 1
253        call constriction(porefill(ilay),eta(ilay))
254    enddo
255    do ilay = index_IS,nsoilmx_PEM
256        eta(ilay) = 0.
257    enddo
258endif
259
260! 1. Depth at which pore filling occurs. We solve Eq. 9 from  Schorgofer, Icarus  (2010)
261old_depth_filling = depth_filling
262
263call deriv1(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dz_psat)
264
265do ilay = 1,index_tmp
266    Jdry = (psat_soil(ilay) - pwat_surf)/mlayer_PEM(ilay) ! left member of Eq. 9 from Schorgofer, Icarus  (2010)
267    Jsat =  eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010)
268    if (Jdry - Jsat <= 0) then
269        index_filling = ilay
270        exit
271    endif
272enddo
273
274if (index_filling == 1) then
275    depth_filling = mlayer_PEM(1)
276else if (index_filling > 1) then
277    Jdry_prevlay = (psat_soil(index_filling - 1) - pwat_surf)/mlayer_PEM(index_filling - 1)
278    Jsat_prevlay = eta(index_filling - 1)*dz_psat(index_filling - 1)
279    call findroot(Jdry - Jsat,Jdry_prevlay - Jsat_prevlay,mlayer_PEM(index_filling),mlayer_PEM(index_filling - 1),depth_filling)
280endif
281
282! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010))
283! 2.0 preliminary: depth to shallowest  ice (discontinuity at interface)
284index_firstice = -1
285do ilay = 1,nsoilmx_PEM
286    if (porefill(ilay) <= 0.) then
287        index_firstice = ilay  ! first point with ice
288        exit
289    endif
290enddo
291
292! 2.1: now we can computeCompute d_z (eta* d_z(rho))
293call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM - 1),dz_eta)
294
295if (index_firstice > 0 .and. index_firstice < nsoilmx_PEM - 2) call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice))
296
297call deriv2_simple(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dzz_psat)
298dz_etadz_rho = pvap2rho*(dz_eta*dz_psat + eta*dzz_psat)
299
300! 3. Ice table boundary due to geothermal heating
301if (index_IS > 0) index_geothermal = -1
302if (index_geothermal < 0) depth_geothermal = -1.
303if ((index_geothermal > 0).and.(index_IS < 0)) then ! Eq. 21 from Schorfoger, Icarus (2010)
304    index_geothermal = -1
305    do ilay = 2,nsoilmx_PEM
306        if (dz_psat(ilay) > 0.) then  ! first point with reversed flux
307            index_geothermal = ilay
308            call findroot(dz_psat(ilay - 1),dz_psat(ilay),mlayer_PEM(ilay - 1),mlayer_PEM(ilay),depth_geothermal)
309            exit
310        endif
311    enddo
312else
313    index_geothermal = -1
314endif
315if (index_geothermal > 0 .and. index_IS < 0) then ! Eq. 24 from Schorgofer, Icarus (2010)
316    call  colint(porefill/eta,mlayer_PEM,nsoilmx_PEM,index_geothermal - 1,nsoilmx_PEM,massfillabove)
317    index_tmp = -1
318    do ilay = index_geothermal,nsoilmx_PEM
319        if (minval(eta(ilay:nsoilmx_PEM)).le.0.) cycle ! eta=0 means completely full
320        call colint(porefill/eta,mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter)
321        if (massfillafter < dz_psat(ilay)*pvap2rho*B) then ! usually executes on i=typeG
322            if (ilay > index_geothermal) then
323!                write(34,*) '# adjustment to geotherm depth by',ilay-index_geothermal
324                call findroot(dz_psat(ilay - 1)*pvap2rho*B - massfillabove, &
325                              dz_psat(ilay)*pvap2rho*B - massfillafter,mlayer_PEM(ilay - 1),mlayer_PEM(ilay),depth_geothermal)
326!                if (depth_geothermal > mlayer_PEM(ilay) .or. depth_geothermal < mlayer_PEM(ilay - 1)) write(*,*) '# WARNING: zdepthG interpolation failed',ilay,mlayer_PEM(ilay - 1),depth_geothermal,mlayer_PEM(ilay)
327              index_tmp = ilay
328           endif
329           ! otherwise leave depth_geothermal unchanged
330           exit
331        endif
332        massfillabove = massfillafter
333    enddo
334    if (index_tmp > 0) index_geothermal = index_tmp
335end if
336
337END SUBROUTINE find_layering_icetable
338
339!-----------------------------------------------------------------------
340SUBROUTINE constriction(porefill,eta)
341
342!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
343!!!
344!!! Purpose: Compute the constriction of vapor flux by pore ice
345!!!
346!!! Author: LL
347!!!
348!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
349
350implicit none
351
352real, intent(in) :: porefill ! pore filling fraction
353real, intent(out) :: eta ! constriction
354
355if (porefill <= 0.) then
356    eta = 1.
357else if (0 < porefill .and. porefill < 1.) then
358    eta = (1-porefill)**2 ! Hudson et al., JGR, 2009
359else
360    eta = 0.
361endif
362
363END SUBROUTINE constriction
364
365!-----------------------------------------------------------------------
366SUBROUTINE compute_Tice_pem(nsoil, ptsoil, ptsurf, ice_depth, Tice)
367!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
368!!!
369!!! Purpose: Compute subsurface ice temperature by interpolating the temperatures between the two adjacent cells.
370!!!
371!!! Author: LL
372!!!
373!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
374
375use comsoil_h_PEM, only: layer_PEM, mlayer_PEM
376
377implicit none
378
379! Inputs
380! ------
381integer,                 intent(in) :: nsoil     ! Number of soil layers
382real, dimension(nsoil),  intent(in) :: ptsoil    ! Soil temperature (K)
383real,                    intent(in) :: ptsurf    ! Soil temperature (K)
384real,                    intent(in) :: ice_depth ! Ice depth (m)
385
386! Outputs
387! ------
388real, intent(out) :: Tice ! Ice temperatures (K)
389
390! Locals
391! ------
392integer :: ik       ! Loop variables
393integer :: indexice ! Index of the ice
394
395! Code
396!-----
397indexice = -1
398if (ice_depth >= mlayer_PEM(nsoil - 1)) then
399    Tice = ptsoil(nsoil)
400else
401    if(ice_depth < mlayer_PEM(0)) then
402        indexice = 0.
403    else     
404        do ik = 0,nsoil - 2 ! go through all the layers to find the ice locations
405            if (mlayer_PEM(ik) <= ice_depth .and. mlayer_PEM(ik + 1) > ice_depth) then
406                indexice = ik + 1
407                exit
408            endif
409        enddo
410    endif
411    if (indexice < 0) then
412        call abort_physic("compute_Tice_pem","subsurface ice is below the last soil layer",1)
413    else
414        if(indexice >= 1) then ! Linear inteprolation between soil temperature
415            Tice = (ptsoil(indexice) - ptsoil(indexice + 1))/(mlayer_PEM(indexice - 1) - mlayer_PEM(indexice))*(ice_depth - mlayer_PEM(indexice)) + ptsoil(indexice + 1)
416        else ! Linear inteprolation between the 1st soil temperature and the surface temperature
417            Tice = (ptsoil(1) - ptsurf)/mlayer_PEM(0)*(ice_depth - mlayer_PEM(0)) + ptsoil(1)
418        endif ! index ice >= 1
419    endif !indexice < 0
420endif ! icedepth > mlayer_PEM(nsoil - 1)
421
422END SUBROUTINE compute_Tice_pem
423
424END MODULE ice_table_mod
Note: See TracBrowser for help on using the repository browser.