source: trunk/LMDZ.COMMON/libf/evolution/ice_table.F90 @ 4135

Last change on this file since 4135 was 4135, checked in by jbclement, 37 hours ago

PEM:
Relocate Mars-specific parameters out of "planet" module and into the modules that own their physics domain.
JBC
File size: 26.2 KB
Line 
1MODULE ice_table
2!-----------------------------------------------------------------------
3! NAME
4!     ice_table
5!
6! DESCRIPTION
7!     Ice table variables and methods to compute it (dynamic and static).
8!
9! AUTHORS & DATE
10!     L. Lange, 02/2023
11!     JB Clement, 2023-2025
12!
13! NOTES
14!
15!-----------------------------------------------------------------------
16
17! DEPENDENCIES
18! ------------
19use numerics, only: dp, di, k4
20
21! DECLARATION
22! -----------
23implicit none
24
25! PARAMETERS
26! ----------
27logical(k4), protected :: icetable_equilibrium ! Flag to compute the icetable depth according to equilibrium
28logical(k4), protected :: icetable_dynamic     ! Flag to compute the icetable depth according to the dynamic method
29
30! VARIABLES
31! ---------
32real(dp), allocatable, dimension(:,:)   :: icetable_depth     ! Depth of the ice table [m]
33real(dp), allocatable, dimension(:,:)   :: icetable_thickness ! Thickness of the ice table [m]
34real(dp), allocatable, dimension(:,:,:) :: ice_porefilling    ! Amount of porefilling in each layer in each grid [m^3/m^3]
35
36contains
37!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
38
39!=======================================================================
40SUBROUTINE ini_ice_table()
41!-----------------------------------------------------------------------
42! NAME
43!     ini_ice_table
44!
45! DESCRIPTION
46!     Allocate module arrays for ice table fields.
47!
48! AUTHORS & DATE
49!     L. Lange
50!     JB Clement, 2023-2025
51!
52! NOTES
53!
54!-----------------------------------------------------------------------
55
56! DEPENDENCIES
57! ------------
58use geometry, only: ngrid, nslope, nsoil
59
60! DECLARATION
61! -----------
62implicit none
63
64! CODE
65! ----
66allocate(icetable_depth(ngrid,nslope))
67allocate(icetable_thickness(ngrid,nslope))
68allocate(ice_porefilling(ngrid,nsoil,nslope))
69icetable_depth(:,:) = 0._dp
70icetable_thickness(:,:) = 0._dp
71ice_porefilling(:,:,:) = 0._dp
72
73END SUBROUTINE ini_ice_table
74!=======================================================================
75
76!=======================================================================
77SUBROUTINE end_ice_table
78!-----------------------------------------------------------------------
79! NAME
80!     end_ice_table
81!
82! DESCRIPTION
83!     Deallocate ice table arrays.
84!
85! AUTHORS & DATE
86!     L. Lange
87!     JB Clement, 2023-2025
88!
89! NOTES
90!
91!-----------------------------------------------------------------------
92
93! DECLARATION
94! -----------
95implicit none
96
97! CODE
98! ----
99if (allocated(icetable_depth)) deallocate(icetable_depth)
100if (allocated(icetable_thickness)) deallocate(icetable_thickness)
101if (allocated(ice_porefilling)) deallocate(ice_porefilling)
102
103END SUBROUTINE end_ice_table
104!=======================================================================
105
106!=======================================================================
107SUBROUTINE set_ice_table_config(icetable_equilibrium_in,icetable_dynamic_in)
108!-----------------------------------------------------------------------
109! NAME
110!     set_ice_table_config
111!
112! DESCRIPTION
113!     Setter for 'ice_table' configuration parameters.
114!
115! AUTHORS & DATE
116!     JB Clement, 02/2026
117!
118! NOTES
119!
120!-----------------------------------------------------------------------
121
122! DEPENDENCIES
123! ------------
124use utility, only: bool2str
125use display, only: print_msg, LVL_NFO, LVL_WRN
126
127! DECLARATION
128! -----------
129implicit none
130
131! ARGUMENTS
132! ---------
133logical(k4), intent(in) :: icetable_equilibrium_in, icetable_dynamic_in
134
135! CODE
136! ----
137icetable_equilibrium = icetable_equilibrium_in
138icetable_dynamic = icetable_dynamic_in
139if (icetable_equilibrium .and. icetable_dynamic) then
140    call print_msg('Ice table is asked to be computed both by the equilibrium and dynamic method. Then, the dynamic method is chosen.',LVL_WRN)
141    icetable_equilibrium = .false.
142end if
143call print_msg('icetable_equilibrium = '//bool2str(icetable_equilibrium_in),LVL_NFO)
144call print_msg('icetable_dynamic     = '//bool2str(icetable_dynamic_in),LVL_NFO)
145
146END SUBROUTINE set_ice_table_config
147!=======================================================================
148
149!=======================================================================
150SUBROUTINE computeice_table_equilibrium(watercaptag,rhowatersurf_avg,rhowatersoil_avg,regolith_inertia,ice_table_depth,ice_table_thickness)
151!-----------------------------------------------------------------------
152! NAME
153!     computeice_table_equilibrium
154!
155! DESCRIPTION
156!     Compute the ice table depth knowing the yearly average water
157!     density at the surface and at depth. Computations are made following
158!     the methods in Schorghofer et al., 2005.
159!
160! AUTHORS & DATE
161!     L. Lange
162!     JB Clement, 12/12/2025
163!
164! NOTES
165!     This subroutine only gives the ice table at equilibrium and does
166!     not consider exchange with the atmosphere.
167!-----------------------------------------------------------------------
168
169! DEPENDENCIES
170! ------------
171use geometry,           only: ngrid, nslope, nsoil
172use maths,              only: findroot
173use soil,               only: mlayer, layer ! Depth of the vertical grid
174use soil_therm_inertia, only: get_ice_TI
175
176! DECLARATION
177! -----------
178implicit none
179
180! ARGUMENTS
181! ---------
182logical(k4), dimension(:),     intent(in)  :: watercaptag         ! Boolean to check the presence of a perennial glacier
183real(dp),    dimension(:,:),   intent(in)  :: rhowatersurf_avg    ! Water density at the surface, yearly averaged [kg/m^3]
184real(dp),    dimension(:,:,:), intent(in)  :: rhowatersoil_avg    ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged [kg/m^3]
185real(dp),    dimension(:,:),   intent(in)  :: regolith_inertia    ! Thermal inertia of the regolith layer [SI]
186real(dp),    dimension(:,:),   intent(out) :: ice_table_depth     ! Ice table depth [m]
187real(dp),    dimension(:,:),   intent(out) :: ice_table_thickness ! Ice table thickness [m]
188
189! LOCAL VARIABLES
190! ---------------
191integer(di)                       :: ig, islope, isoil, isoilend
192real(dp), dimension(nsoil)        :: diff_rho                ! Difference of water vapor density between the surface and at depth [kg/m^3]
193real(dp)                          :: ice_table_end           ! Depth of the end of the ice table  [m]
194real(dp), dimension(ngrid,nslope) :: previous_icetable_depth ! Ice table computed at previous ice depth [m]
195real(dp)                          :: stretch                 ! Stretch factor to improve the convergence of the ice table
196real(dp)                          :: wice_inertia            ! Water Ice thermal Inertia [USI]
197
198! CODE
199! ----
200previous_icetable_depth = ice_table_depth
201do ig = 1,ngrid
202    if (watercaptag(ig)) then
203        ice_table_depth(ig,:) = 0._dp
204        ice_table_thickness(ig,:) = layer(nsoil) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
205    else
206        do islope = 1,nslope
207            ice_table_depth(ig,islope) = -1._dp
208            ice_table_thickness(ig,islope) = 0._dp
209            do isoil = 1,nsoil
210                diff_rho(isoil) = rhowatersurf_avg(ig,islope) - rhowatersoil_avg(ig,isoil,islope)
211            end do
212            if (diff_rho(1) > 0.) then ! ice is at the surface
213                ice_table_depth(ig,islope) = 0._dp
214                do isoilend = 2,nsoil - 1
215                    if (diff_rho(isoilend) > 0._dp .and. diff_rho(isoilend + 1) < 0._dp) then
216                        call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer(isoilend - 1),mlayer(isoilend),ice_table_end)
217                        ice_table_thickness(ig,islope) = ice_table_end - ice_table_depth(ig,islope)
218                        exit
219                    end if
220                end do
221            else
222                do isoil = 1,nsoil - 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
223                    if (diff_rho(isoil) < 0._dp .and. diff_rho(isoil + 1) > 0._dp) then
224                        call findroot(diff_rho(isoil),diff_rho(isoil + 1),mlayer(isoil),mlayer(isoil + 1),ice_table_depth(ig,islope))
225                        ! Now let's find the end of the ice table
226                        ice_table_thickness(ig,islope) = layer(nsoil) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
227                        do isoilend = isoil + 1,nsoil - 1
228                            if (diff_rho(isoilend) > 0._dp .and. diff_rho(isoilend + 1) < 0._dp) then
229                                call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer(isoilend - 1),mlayer(isoilend),ice_table_end)
230                                ice_table_thickness(ig,islope) = ice_table_end - ice_table_depth(ig,islope)
231                                exit
232                            end if
233                        end do
234                    end if ! diff_rho(z) <0 & diff_rho(z+1) > 0
235                end do
236            end if ! diff_rho(1) > 0
237        end do
238    end if ! watercaptag
239end do
240
241! Small trick to speed up the convergence, Oded's idea.
242do islope = 1,nslope
243    do ig = 1,ngrid
244        if (ice_table_depth(ig,islope) > previous_icetable_depth(ig,islope) .and. previous_icetable_depth(ig,islope) >= 0._dp) then
245            call get_ice_TI(.false.,1._dp,regolith_inertia(ig,islope),wice_inertia)
246            stretch = (regolith_inertia(ig,islope)/wice_inertia)**2
247            ice_table_thickness(ig,islope) = ice_table_thickness(ig,islope) + (ice_table_depth(ig,islope) - &
248                                             previous_icetable_depth(ig,islope) + (ice_table_depth(ig,islope) - previous_icetable_depth(ig,islope))/stretch)
249            ice_table_depth(ig,islope) = previous_icetable_depth(ig,islope) + (ice_table_depth(ig,islope) - previous_icetable_depth(ig,islope))/stretch
250        end if
251    end do
252end do
253
254END SUBROUTINE computeice_table_equilibrium
255!=======================================================================
256
257!=======================================================================
258SUBROUTINE compute_deltam_icetable(icetable_thickness_old,ice_porefilling_old,tsurf,tsoil,delta_icetable)
259!-----------------------------------------------------------------------
260! NAME
261!     compute_deltam_icetable
262!
263! DESCRIPTION
264!     Compute the mass of H2O that has sublimated from the ice table /
265!     condensed.
266!
267! AUTHORS & DATE
268!     L. Lange
269!     JB Clement, 2023-2025
270!
271! NOTES
272!
273!-----------------------------------------------------------------------
274
275! DEPENDENCIES
276! ------------
277use geometry, only: ngrid, nslope, nsoil
278use soil,     only: mlayer, regolith_porosity
279use slopes,   only: subslope_dist, def_slope_mean
280use maths,    only: pi
281
282! DECLARATION
283! -----------
284implicit none
285
286! ARGUMENTS
287! ---------
288real(dp),    dimension(:,:),   intent(in)  :: icetable_thickness_old   ! Ice table thickness at the previous iteration [m]
289real(dp),    dimension(:,:,:), intent(in)  :: ice_porefilling_old      ! Ice pore filling at the previous iteration [m]
290real(dp),    dimension(:,:),   intent(in)  :: tsurf                    ! Surface temperature [K]
291real(dp),    dimension(:,:,:), intent(in)  :: tsoil                    ! Soil temperature [K]
292real(dp),    dimension(:),     intent(out) :: delta_icetable ! Mass of H2O ice that has been condensed on the ice table / sublimates from the ice table [kg/m^2]
293
294! LOCAL VARIABLES
295! ---------------
296integer(di) :: ig, islope, isoil
297real(dp)    :: Tice ! Ice temperature [k]
298
299! CODE
300! ----
301! Let's compute the amount of ice that has sublimated in each subslope
302if (icetable_equilibrium) then
303    delta_icetable = 0._dp
304    do ig = 1,ngrid
305        do islope = 1,nslope
306            call compute_Tice(tsoil(ig,:,islope),tsurf(ig,islope),icetable_depth(ig,islope),Tice)
307            delta_icetable(ig) = delta_icetable(ig) + regolith_porosity*rho_ice(Tice,'H2O')*(icetable_thickness(ig,islope) - icetable_thickness_old(ig,islope)) & ! convention > 0. <=> it condenses
308                              *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180._dp)
309        end do
310    end do
311else if (icetable_dynamic) then
312    delta_icetable = 0._dp
313    do ig = 1,ngrid
314        do islope = 1,nslope
315            do isoil = 1,nsoil
316                call compute_Tice(tsoil(ig,:,islope),tsurf(ig,islope),mlayer(isoil - 1),Tice)
317                delta_icetable(ig) = delta_icetable(ig) + rho_ice(Tice,'H2O')*(ice_porefilling(ig,isoil,islope) - ice_porefilling_old(ig,isoil,islope)) & ! convention > 0. <=> it condenses
318                                  *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180._dp)
319            end do
320        end do
321    end do
322end if
323
324END SUBROUTINE compute_deltam_icetable
325!=======================================================================
326
327!=======================================================================
328SUBROUTINE 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)
329!-----------------------------------------------------------------------
330! NAME
331!     find_layering_icetable
332!
333! DESCRIPTION
334!     Compute layering between dry soil, pore filling ice and ice
335!     sheet based on Schorghofer, Icarus (2010). Adapted from NS MSIM.
336!
337! AUTHORS & DATE
338!     L. Lange
339!     JB Clement, 2023-2025
340!
341! NOTES
342!
343!-----------------------------------------------------------------------
344
345! DEPENDENCIES
346! ------------
347use soil,           only: mlayer
348use geometry,       only: nsoil
349use maths,          only: deriv1, deriv1_onesided, colint, findroot, deriv2_simple
350use subsurface_ice, only: constriction
351
352! DECLARATION
353! -----------
354implicit none
355
356! ARGUMENTS
357! ---------
358real(dp), dimension(:), intent(in)    :: porefill         ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice
359real(dp), dimension(:), intent(in)    :: psat_soil        ! Soil water pressure at saturation, yearly averaged [Pa]
360real(dp),               intent(in)    :: psat_surf        ! Surface water pressure at saturation, yearly averaged [Pa]
361real(dp),               intent(in)    :: pwat_surf        ! Water vapor pressure  at the surface, not necesseraly at saturation, yearly averaged [Pa]
362real(dp),               intent(in)    :: psat_bottom      ! Boundary conditions for soil vapor pressure [Pa]
363real(dp),               intent(in)    :: B                ! Constant (Eq. 8 from  Schorgofer, Icarus (2010).)
364integer(di),            intent(in)    :: index_IS         ! Index of the soil layer where the ice sheet begins [1]
365real(dp),               intent(inout) :: depth_filling    ! Depth where pore filling begins [m]
366integer(di),            intent(out)   :: index_filling    ! Index where the pore filling begins [1]
367integer(di),            intent(out)   :: index_geothermal ! Index where the ice table stops [1]
368real(dp),               intent(out)   :: depth_geothermal ! Depth where the ice table stops [m]
369real(dp), dimension(:), intent(out)   :: dz_etadz_rho     ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase
370
371! LOCAL VARIABLES
372! ---------------
373real(dp), dimension(nsoil) :: eta                          ! Constriction
374real(dp), dimension(nsoil) :: dz_psat                      ! First derivative of the vapor pressure at saturation
375real(dp), dimension(nsoil) :: dz_eta                       ! \partial z \eta
376real(dp), dimension(nsoil) :: dzz_psat                     ! \partial \partial psat
377integer(di)                :: ilay, index_tmp              ! Index for loop
378real(dp)                   :: old_depth_filling            ! Depth_filling saved
379real(dp)                   :: Jdry                         ! Flux trought the dry layer
380real(dp)                   :: Jsat                         ! Flux trought the ice layer
381real(dp)                   :: Jdry_prevlay, Jsat_prevlay   ! Same but for the previous ice layer
382integer(di)                :: index_firstice               ! First index where ice appears (i.e., f > 0)
383real(dp)                   :: massfillabove, massfillafter ! H2O mass above and after index_geothermal
384real(dp), parameter        :: pvap2rho = 18.e-3/8.314
385
386! CODE
387! ----
388! 0. Compute constriction over the layer
389! Within the ice sheet, constriction is set to 0. Elsewhere, constriction =  (1-porefilling)**2
390if (index_IS < 0) then
391    index_tmp = nsoil
392    do ilay = 1,nsoil
393        eta(ilay) = constriction(porefill(ilay))
394    end do
395else
396    index_tmp = index_IS
397    do ilay = 1,index_IS - 1
398        eta(ilay) = constriction(porefill(ilay))
399    end do
400    do ilay = index_IS,nsoil
401        eta(ilay) = 0._dp
402    end do
403end if
404
405! 1. Depth at which pore filling occurs. We solve Eq. 9 from  Schorgofer, Icarus  (2010)
406old_depth_filling = depth_filling
407
408call deriv1(mlayer,nsoil,psat_soil,psat_surf,psat_bottom,dz_psat)
409
410do ilay = 1,index_tmp
411    Jdry = (psat_soil(ilay) - pwat_surf)/mlayer(ilay) ! left member of Eq. 9 from Schorgofer, Icarus  (2010)
412    Jsat =  eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010)
413    if (Jdry - Jsat <= 0) then
414        index_filling = ilay
415        exit
416    end if
417end do
418
419if (index_filling == 1) then
420    depth_filling = mlayer(1)
421else if (index_filling > 1) then
422    Jdry_prevlay = (psat_soil(index_filling - 1) - pwat_surf)/mlayer(index_filling - 1)
423    Jsat_prevlay = eta(index_filling - 1)*dz_psat(index_filling - 1)
424    call findroot(Jdry - Jsat,Jdry_prevlay - Jsat_prevlay,mlayer(index_filling),mlayer(index_filling - 1),depth_filling)
425end if
426
427! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010))
428! 2.0 preliminary: depth to shallowest ice (discontinuity at interface)
429index_firstice = -1
430do ilay = 1,nsoil
431    if (porefill(ilay) <= 0._dp) then
432        index_firstice = ilay  ! first point with ice
433        exit
434    end if
435end do
436
437! 2.1: now we can compute
438call deriv1(mlayer,nsoil,eta,1.,eta(nsoil - 1),dz_eta)
439
440if (index_firstice > 0 .and. index_firstice < nsoil - 2) call deriv1_onesided(index_firstice,mlayer,nsoil,eta,dz_eta(index_firstice))
441
442call deriv2_simple(mlayer,nsoil,psat_soil,psat_surf,psat_bottom,dzz_psat)
443dz_etadz_rho = pvap2rho*(dz_eta*dz_psat + eta*dzz_psat)
444
445! 3. Ice table boundary due to geothermal heating
446if (index_IS > 0) index_geothermal = -1
447if (index_geothermal < 0) depth_geothermal = -1._dp
448if ((index_geothermal > 0).and.(index_IS < 0)) then ! Eq. 21 from Schorfoger, Icarus (2010)
449    index_geothermal = -1
450    do ilay = 2,nsoil
451        if (dz_psat(ilay) > 0._dp) then  ! first point with reversed flux
452            index_geothermal = ilay
453            call findroot(dz_psat(ilay - 1),dz_psat(ilay),mlayer(ilay - 1),mlayer(ilay),depth_geothermal)
454            exit
455        end if
456    end do
457else
458    index_geothermal = -1
459end if
460if (index_geothermal > 0 .and. index_IS < 0) then ! Eq. 24 from Schorgofer, Icarus (2010)
461    call  colint(porefill/eta,mlayer,nsoil,index_geothermal - 1,nsoil,massfillabove)
462    index_tmp = -1
463    do ilay = index_geothermal,nsoil
464        if (minval(eta(ilay:nsoil)) <= 0._dp) cycle ! eta=0 means completely full
465        call colint(porefill/eta,mlayer,nsoil,ilay,nsoil,massfillafter)
466        if (massfillafter < dz_psat(ilay)*pvap2rho*B) then ! usually executes on i=typeG
467            if (ilay > index_geothermal) then
468                call findroot(dz_psat(ilay - 1)*pvap2rho*B - massfillabove, &
469                              dz_psat(ilay)*pvap2rho*B - massfillafter,mlayer(ilay - 1),mlayer(ilay),depth_geothermal)
470!                if (depth_geothermal > mlayer(ilay) .or. depth_geothermal < mlayer(ilay - 1)) write(*,*) '# WARNING: zdepthG interpolation failed',ilay,mlayer(ilay - 1),depth_geothermal,mlayer(ilay)
471              index_tmp = ilay
472           end if
473           ! otherwise leave depth_geothermal unchanged
474           exit
475        end if
476        massfillabove = massfillafter
477    end do
478    if (index_tmp > 0) index_geothermal = index_tmp
479end if
480
481END SUBROUTINE find_layering_icetable
482!=======================================================================
483
484!=======================================================================
485SUBROUTINE compute_Tice(ptsoil,ptsurf,ice_depth,Tice)
486!-----------------------------------------------------------------------
487! NAME
488!     compute_Tice
489!
490! DESCRIPTION
491!     Compute subsurface ice temperature by interpolating the
492!     temperatures between the two adjacent cells.
493!
494! AUTHORS & DATE
495!     L. Lange
496!     JB Clement, 12/12/2025
497!
498! NOTES
499!
500!-----------------------------------------------------------------------
501
502! DEPENDENCIES
503! ------------
504use geometry, only: nsoil
505use soil,     only: mlayer
506use stoppage, only: stop_clean
507
508! DECLARATION
509! -----------
510implicit none
511
512! ARGUMENTS
513! ---------
514real(dp), dimension(:), intent(in)  :: ptsoil    ! Soil temperature (K)
515real(dp),               intent(in)  :: ptsurf    ! Surface temperature (K)
516real(dp),               intent(in)  :: ice_depth ! Ice depth (m)
517real(dp),               intent(out) :: Tice      ! Ice temperatures (K)
518
519! LOCAL VARIABLES
520! ---------------
521integer(di) :: ik       ! Loop variables
522integer(di) :: indexice ! Index of the ice
523
524! CODE
525! ----
526indexice = -1
527if (ice_depth >= mlayer(nsoil - 1)) then
528    Tice = ptsoil(nsoil)
529else
530    if(ice_depth < mlayer(0)) then
531        indexice = 0
532    else
533        do ik = 0,nsoil - 2 ! go through all the layers to find the ice locations
534            if (mlayer(ik) <= ice_depth .and. mlayer(ik + 1) > ice_depth) then
535                indexice = ik + 1
536                exit
537            end if
538        end do
539    end if
540    if (indexice < 0) then
541        call stop_clean(__FILE__,__LINE__,"subsurface ice is below the last soil layer!",1)
542    else
543        if(indexice >= 1) then ! Linear inteprolation between soil temperature
544            Tice = (ptsoil(indexice) - ptsoil(indexice + 1))/(mlayer(indexice - 1) - mlayer(indexice))*(ice_depth - mlayer(indexice)) + ptsoil(indexice + 1)
545        else ! Linear inteprolation between the 1st soil temperature and the surface temperature
546            Tice = (ptsoil(1) - ptsurf)/mlayer(0)*(ice_depth - mlayer(0)) + ptsoil(1)
547        end if ! index ice >= 1
548    end if !indexice < 0
549end if ! icedepth > mlayer(nsoil - 1)
550
551END SUBROUTINE compute_Tice
552!=======================================================================
553
554!=======================================================================
555FUNCTION rho_ice(T,ice) RESULT(rho)
556!-----------------------------------------------------------------------
557! NAME
558!     rho_ice
559!
560! DESCRIPTION
561!     Compute subsurface ice density.
562!
563! AUTHORS & DATE
564!     JB Clement, 12/12/2025
565!
566! NOTES
567!
568!-----------------------------------------------------------------------
569
570! DEPENDENCIES
571! ------------
572use stoppage, only: stop_clean
573
574! DECLARATION
575! -----------
576implicit none
577
578! ARGUMENTS
579! ---------
580real(dp),     intent(in) :: T
581character(3), intent(in) :: ice
582
583! LOCAL VARIABLES
584! ---------------
585real(dp) :: rho
586
587! CODE
588! ----
589select case (trim(adjustl(ice)))
590    case ('H2O')
591        rho = -3.5353e-4_dp*T**2 + 0.0351_dp*T + 933.5030_dp ! Rottgers 2012
592    case ('CO2')
593        rho = (1.72391_dp - 2.53e-4_dp*T - 2.87_dp*1.e-7_dp*T**2)*1.e3_dp ! Mangan et al. 2017
594    case default
595        call stop_clean(__FILE__,__LINE__,"type of ice unknown!",1)
596end select
597
598END FUNCTION rho_ice
599!=======================================================================
600
601!=======================================================================
602SUBROUTINE evolve_ice_table(h2o_ice,h2o_surfdensity_avg,h2o_soildensity_avg,tsoil_avg,tsurf_avg,delta_icetable, &
603                            q_h2o_ts,ps_avg,icetable_depth,icetable_thickness,ice_porefilling,icetable_depth_old)
604!-----------------------------------------------------------------------
605! NAME
606!     evolve_ice_table
607!
608! DESCRIPTION
609!     Update the subsurface ice table using either the equilibrium or
610!     dynamic method and calculate the resulting H2O mass exchange.
611!
612! AUTHORS & DATE
613!     JB Clement, 12/2025
614!     C. Metz, 02/2026
615!-----------------------------------------------------------------------
616
617! DEPENDENCIES
618! ------------
619use geometry,       only: ngrid, nsoil, nslope, nday
620use slopes,         only: subslope_dist, def_slope_mean
621use surf_ice,       only: threshold_h2oice_cap
622use maths,          only: pi
623use soil,           only: TI
624use subsurface_ice, only: dyn_ss_ice_m
625use display,        only: print_msg, LVL_NFO
626
627! DECLARATION
628! -----------
629implicit none
630
631! ARGUMENTS
632! ---------
633real(dp), dimension(:),     intent(in)    :: ps_avg
634real(dp), dimension(:,:),   intent(in)    :: h2o_ice, h2o_surfdensity_avg, tsurf_avg, q_h2o_ts
635real(dp), dimension(:,:,:), intent(in)    :: h2o_soildensity_avg, tsoil_avg
636real(dp), dimension(:,:),   intent(out)   :: icetable_depth_old
637real(dp), dimension(:),     intent(out)   :: delta_icetable
638real(dp), dimension(:,:),   intent(inout) :: icetable_depth, icetable_thickness
639real(dp), dimension(:,:,:), intent(inout) :: ice_porefilling
640
641! LOCAL VARIABLES
642! ---------------
643integer(di)                                :: i, islope
644logical(k4), dimension(ngrid)              :: is_h2o_perice
645real(dp),    dimension(ngrid,nslope)       :: icetable_thickness_old
646real(dp),    dimension(ngrid,nsoil,nslope) :: ice_porefilling_old
647
648! CODE
649! ----
650do i = 1,ngrid
651    if (sum(h2o_ice(i,:)*subslope_dist(i,:)/cos(pi*def_slope_mean(:)/180._dp)) > threshold_h2oice_cap) then
652        ! There is enough ice to be considered as an infinite reservoir
653        is_h2o_perice(i) = .true.
654    else
655        ! Too little ice to be considered as an infinite reservoir so ice is transferred to the frost
656        is_h2o_perice(i) = .false.
657    end if
658end do
659if (icetable_equilibrium) then
660    call print_msg("> Evolution of ice table (equilibrium method)",LVL_NFO)
661    icetable_thickness_old(:,:) = icetable_thickness(:,:)
662    call computeice_table_equilibrium(is_h2o_perice,h2o_surfdensity_avg,h2o_soildensity_avg,TI(:,1,:),icetable_depth,icetable_thickness)
663else if (icetable_dynamic) then
664    call print_msg("> Evolution of ice table (dynamic method)",LVL_NFO)
665    ice_porefilling_old(:,:,:) = ice_porefilling(:,:,:)
666    icetable_depth_old(:,:) = icetable_depth(:,:)
667    do i = 1,ngrid
668        do islope = 1,nslope
669            call dyn_ss_ice_m(icetable_depth(i,islope),tsurf_avg(i,islope), tsoil_avg(i,:,islope),nsoil,TI(i,1,islope),ps_avg, &
670                              (/sum(q_h2o_ts(i,:))/real(nday,dp)/),ice_porefilling(i,:,islope),ice_porefilling(i,:,islope),icetable_depth(i,islope))
671        end do
672    end do
673end if
674
675! Mass of H2O exchange between the ssi and the atmosphere
676call compute_deltam_icetable(icetable_thickness_old,ice_porefilling_old,tsurf_avg,tsoil_avg,delta_icetable)
677
678END SUBROUTINE evolve_ice_table
679!=======================================================================
680
681END MODULE ice_table
Note: See TracBrowser for help on using the repository browser.