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

Last change on this file since 4138 was 4138, checked in by jbclement, 20 hours ago

PEM:

  • Move ice table variables from "ice_table" to the main program.
  • Merge "job_id_mod" and "job_timelimit_mod" into "job" which is relocated to the PEM folder.
  • Rename local variables in procedures to avoid masking variables in parent scope.
  • Few cleanings to delete remaining PEM-external "include" and "use".

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