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

Last change on this file since 4175 was 4170, checked in by jbclement, 6 days ago

PEM:

  • Deletion of 'flux_ssice' ('zqdsdif_tot') from the "startfi.nc" as it is not needed.
  • Using the yearly average flux for the sublimating subsurface ice instead of the value at last timestep.
  • Keeping a clear separation between the subsurface ice flux and the surface ice tendency.
  • Making sure that subsurface ice depth is well given to the PCM at the end of the PEM (ice table VS layering).

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