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

Last change on this file since 3320 was 3320, checked in by jbclement, 10 months ago

PEM:
Small correction of a bug in the reading of "run_PEM.def" + some cleanings.
JBC
File size: 17.5 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!----------------------------------------
18  contains
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
40 END 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 :: diff_rho(nsoil_PEM)                                    ! 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 :: previous_icetable_depth(ngrid,nslope)                  ! 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! ------
81     previous_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
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
149use vdifc_mod,             only: compute_Tice
150#ifndef CPP_STD
151    use comcstfi_h,   only: pi
152#else
153    use comcstfi_mod, only: pi
154#endif
155
156implicit none
157
158!   Inputs
159! --------
160integer,                                 intent(in) :: ngrid, nslope, nsoil_PEM   ! Size of the physical grid, number of subslope, number of soil layer in the PEM
161real, dimension(ngrid,nslope),           intent(in) :: former_ice_table_thickness ! ice table thickness at the former iteration [m]
162real, dimension(ngrid,nslope),           intent(in) :: new_ice_table_thickness    ! ice table thickness at the current iteration [m]
163real, dimension(ngrid,nslope),           intent(in) :: ice_table_depth            ! ice table depth [m]
164real, dimension(ngrid,nslope),           intent(in) :: tsurf                      ! Surface temperature [K]
165real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil                      ! Soil temperature [K]
166!   Outputs
167! ---------
168real, 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]
169!   Locals
170!---------
171integer :: ig, islope, ilay, iref     ! loop index
172real, dimension(ngrid,nslope) :: rho  ! density of water ice [kg/m^3]
173real, dimension(ngrid,nslope) :: Tice ! ice temperature [k]
174!   Code
175! ------
176rho = 0.
177Tice = 0.
178!1. First let's compute Tice using a linear interpolation between the mlayer level
179do ig = 1,ngrid
180    do islope = 1,nslope
181        call compute_Tice(nsoil_PEM,tsoil(ig,:,islope),tsurf(ig,islope),ice_table_depth(ig,islope),Tice(ig,islope))
182        rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2 + 0.0351* Tice(ig,islope) + 933.5030 ! Rottgers, 2012
183    enddo
184enddo
185
186!2. Let's compute the amount of ice that has sublimated in each subslope
187do ig = 1,ngrid
188    do islope = 1,nslope
189        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
190                          *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
191    enddo
192enddo
193
194END SUBROUTINE
195
196!----------------------------------------
197
198SUBROUTINE 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)
199!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
200!!!
201!!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM
202!!!
203!!! Author: LL
204!!!
205!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
206use comsoil_h_PEM, only: nsoilmx_PEM, mlayer_PEM
207use math_mod,      only: deriv1, deriv1_onesided, colint, findroot, deriv2_simple
208
209implicit none
210
211! Inputs
212! ------
213real, dimension(nsoilmx_PEM), intent(in) :: porefill    ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice
214real, dimension(nsoilmx_PEM), intent(in) :: psat_soil   ! Soil water pressure at saturation, yearly averaged [Pa]
215real,                         intent(in) :: psat_surf   ! surface water pressure at saturation, yearly averaged [Pa]
216real,                         intent(in) :: pwat_surf   ! Water vapor pressure  at the surface, not necesseraly at saturation, yearly averaged [Pa]
217real,                         intent(in) :: psat_bottom ! Boundary conditions for soil vapor pressure [Pa]
218real,                         intent(in) :: B           ! constant (Eq. 8 from  Schorgofer, Icarus (2010).)
219integer,                      intent(in) :: index_IS    ! index of the soil layer where the ice sheet begins [1]
220real, intent(inout) :: depth_filling ! depth where pore filling begins [m]
221! Outputs
222! -------
223integer,                      intent(out) :: index_filling    ! index where the pore filling begins [1]
224integer,                      intent(out) :: index_geothermal ! index where the ice table stops [1]
225real,                         intent(out) :: depth_geothermal ! depth where the ice table stops [m]
226real, dimension(nsoilmx_PEM), intent(out) :: dz_etadz_rho     ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase
227!   Locals
228!---------
229real, dimension(nsoilmx_PEM) :: eta                          ! constriction
230real, dimension(nsoilmx_PEM) :: dz_psat                      ! first derivative of the vapor pressure at saturation
231real, dimension(nsoilmx_PEM) :: dz_eta                       ! \partial z \eta
232real, dimension(nsoilmx_PEM) :: dzz_psat                     ! \partial \partial psat
233integer                      :: ilay, index_tmp              ! index for loop
234real                         :: old_depth_filling            ! depth_filling saved
235real                         :: Jdry                         ! flux trought the dry layer
236real                         :: Jsat                         ! flux trought the ice layer
237real                         :: Jdry_prevlay, Jsat_prevlay   ! same but for the previous ice layer
238integer                      :: index_firstice               ! first index where ice appears (i.e., f > 0)
239real                         :: massfillabove, massfillafter ! h2O mass above and after index_geothermal
240!   Constant
241!-----------
242real, parameter :: pvap2rho = 18.e-3/8.314
243!   Code
244!-------
245! 0. Compute constriction over the layer
246! Within the ice sheet, constriction is set to 0. Elsewhere, constriction =  (1-porefilling)**2
247if (index_IS.lt.0) then
248    index_tmp = nsoilmx_PEM
249    do ilay = 1,nsoilmx_PEM
250        call constriction(porefill(ilay),eta(ilay))
251    enddo
252else
253    index_tmp = index_IS
254    do ilay = 1,index_IS-1
255        call constriction(porefill(ilay),eta(ilay))
256    enddo
257    do ilay = index_IS,nsoilmx_PEM
258        eta(ilay) = 0.
259    enddo
260endif
261
262! 1. Depth at which pore filling occurs. We solve Eq. 9 from  Schorgofer, Icarus  (2010)
263old_depth_filling = depth_filling
264
265call deriv1(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dz_psat)
266
267do ilay = 1,index_tmp
268    Jdry = (psat_soil(ilay) - pwat_surf)/mlayer_PEM(ilay) ! left member of Eq. 9 from Schorgofer, Icarus  (2010)
269    Jsat =  eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010)
270    if (Jdry - Jsat <= 0) then
271        index_filling = ilay
272        exit
273    endif
274enddo
275
276if (index_filling == 1) then
277    depth_filling = mlayer_PEM(1)
278else if (index_filling > 1) then
279    Jdry_prevlay = (psat_soil(index_filling - 1) - pwat_surf)/mlayer_PEM(index_filling - 1)
280    Jsat_prevlay = eta(index_filling - 1)*dz_psat(index_filling - 1)
281    call findroot(Jdry - Jsat,Jdry_prevlay - Jsat_prevlay,mlayer_PEM(index_filling),mlayer_PEM(index_filling - 1),depth_filling)
282endif
283
284! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010))
285! 2.0 preliminary: depth to shallowest  ice (discontinuity at interface)
286index_firstice = -1
287do ilay = 1,nsoilmx_PEM
288    if (porefill(ilay) <= 0.) then
289        index_firstice = ilay  ! first point with ice
290        exit
291    endif
292enddo
293
294! 2.1: now we can computeCompute d_z (eta* d_z(rho))
295call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM - 1),dz_eta)
296
297if (index_firstice > 0 .and. index_firstice < nsoilmx_PEM - 2) call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice))
298
299call deriv2_simple(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dzz_psat)
300dz_etadz_rho = pvap2rho*(dz_eta*dz_psat + eta*dzz_psat)
301
302! 3. Ice table boundary due to geothermal heating
303if (index_IS > 0) index_geothermal = -1
304if (index_geothermal < 0) depth_geothermal = -1.
305if ((index_geothermal > 0).and.(index_IS < 0)) then ! Eq. 21 from Schorfoger, Icarus (2010)
306    index_geothermal = -1
307    do ilay = 2,nsoilmx_PEM
308        if (dz_psat(ilay) > 0.) then  ! first point with reversed flux
309            index_geothermal = ilay
310            call findroot(dz_psat(ilay - 1),dz_psat(ilay),mlayer_PEM(ilay - 1),mlayer_PEM(ilay),depth_geothermal)
311            exit
312        endif
313    enddo
314else
315    index_geothermal = -1
316endif
317if (index_geothermal > 0 .and. index_IS < 0) then ! Eq. 24 from Schorgofer, Icarus (2010)
318    call  colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,index_geothermal-1,nsoilmx_PEM,massfillabove)
319    index_tmp = -1
320    do ilay = index_geothermal,nsoilmx_PEM
321        if (minval(eta(ilay:nsoilmx_PEM)).le.0.) cycle ! eta=0 means completely full
322        call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter)
323        if (massfillafter < dz_psat(ilay)*pvap2rho*B) then ! usually executes on i=typeG
324            if (ilay > index_geothermal) then
325!                write(34,*) '# adjustment to geotherm depth by',ilay-index_geothermal
326                call findroot(dz_psat(ilay - 1)*pvap2rho*B - massfillabove, &
327                              dz_psat(ilay)*pvap2rho*B - massfillafter,mlayer_PEM(ilay - 1),mlayer_PEM(ilay),depth_geothermal)
328!                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)
329              index_tmp = ilay
330           endif
331           ! otherwise leave depth_geothermal unchanged
332           exit
333        endif
334        massfillabove = massfillafter
335    enddo
336    if (index_tmp > 0) index_geothermal = index_tmp
337end if
338
339END SUBROUTINE
340
341!----------------------------------------
342SUBROUTINE constriction(porefill,eta)
343
344!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
345!!!
346!!! Purpose: Compute the constriction of vapor flux by pore ice
347!!!
348!!! Author: LL
349!!!
350!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
351
352implicit none
353
354real, intent(in) :: porefill ! pore filling fraction
355real, intent(out) :: eta ! constriction
356
357if (porefill <= 0.) then
358    eta = 1.
359else if (0 < porefill .and. porefill < 1.) then
360    eta = (1-porefill)**2 ! Hudson et al., JGR, 2009
361else
362    eta = 0.
363endif
364
365END SUBROUTINE
366
367END MODULE
Note: See TracBrowser for help on using the repository browser.