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

Last change on this file since 4092 was 4074, checked in by jbclement, 2 weeks ago

PEM:

  • Correct management of H2O ice tendency in 1D when there is not enough ice anymore.
  • Clean initialization of allocatable module arrays (especially needed when no slope)
  • One more renaming for consistency + few small updates thoughout the code.

JBC
File size: 26.1 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
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.')
141    icetable_equilibrium = .false.
142end if
143call print_msg('icetable_equilibrium = '//bool2str(icetable_equilibrium_in))
144call print_msg('icetable_dynamic     = '//bool2str(icetable_dynamic_in))
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
279use slopes,   only: subslope_dist, def_slope_mean
280use planet,   only: porosity
281use maths,    only: pi
282
283! DECLARATION
284! -----------
285implicit none
286
287! ARGUMENTS
288! ---------
289real(dp),    dimension(:,:),   intent(in)  :: icetable_thickness_old   ! Ice table thickness at the previous iteration [m]
290real(dp),    dimension(:,:,:), intent(in)  :: ice_porefilling_old      ! Ice pore filling at the previous iteration [m]
291real(dp),    dimension(:,:),   intent(in)  :: tsurf                    ! Surface temperature [K]
292real(dp),    dimension(:,:,:), intent(in)  :: tsoil                    ! Soil temperature [K]
293real(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]
294
295! LOCAL VARIABLES
296! ---------------
297integer(di) :: ig, islope, isoil
298real(dp)    :: Tice ! Ice temperature [k]
299
300! CODE
301! ----
302! Let's compute the amount of ice that has sublimated in each subslope
303if (icetable_equilibrium) then
304    delta_icetable = 0._dp
305    do ig = 1,ngrid
306        do islope = 1,nslope
307            call compute_Tice(tsoil(ig,:,islope),tsurf(ig,islope),icetable_depth(ig,islope),Tice)
308            delta_icetable(ig) = delta_icetable(ig) + porosity*rho_ice(Tice,'H2O')*(icetable_thickness(ig,islope) - icetable_thickness_old(ig,islope)) & ! convention > 0. <=> it condenses
309                              *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180._dp)
310        end do
311    end do
312else if (icetable_dynamic) then
313    delta_icetable = 0._dp
314    do ig = 1,ngrid
315        do islope = 1,nslope
316            do isoil = 1,nsoil
317                call compute_Tice(tsoil(ig,:,islope),tsurf(ig,islope),mlayer(isoil - 1),Tice)
318                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
319                                  *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180._dp)
320            end do
321        end do
322    end do
323end if
324
325END SUBROUTINE compute_deltam_icetable
326!=======================================================================
327
328!=======================================================================
329SUBROUTINE 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)
330!-----------------------------------------------------------------------
331! NAME
332!     find_layering_icetable
333!
334! DESCRIPTION
335!     Compute layering between dry soil, pore filling ice and ice
336!     sheet based on Schorghofer, Icarus (2010). Adapted from NS MSIM.
337!
338! AUTHORS & DATE
339!     L. Lange
340!     JB Clement, 2023-2025
341!
342! NOTES
343!
344!-----------------------------------------------------------------------
345
346! DEPENDENCIES
347! ------------
348use soil,           only: mlayer
349use geometry,       only: nsoil
350use maths,          only: deriv1, deriv1_onesided, colint, findroot, deriv2_simple
351use subsurface_ice, only: constriction
352
353! DECLARATION
354! -----------
355implicit none
356
357! ARGUMENTS
358! ---------
359real(dp), dimension(:), intent(in)    :: porefill         ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice
360real(dp), dimension(:), intent(in)    :: psat_soil        ! Soil water pressure at saturation, yearly averaged [Pa]
361real(dp),               intent(in)    :: psat_surf        ! Surface water pressure at saturation, yearly averaged [Pa]
362real(dp),               intent(in)    :: pwat_surf        ! Water vapor pressure  at the surface, not necesseraly at saturation, yearly averaged [Pa]
363real(dp),               intent(in)    :: psat_bottom      ! Boundary conditions for soil vapor pressure [Pa]
364real(dp),               intent(in)    :: B                ! Constant (Eq. 8 from  Schorgofer, Icarus (2010).)
365integer(di),            intent(in)    :: index_IS         ! Index of the soil layer where the ice sheet begins [1]
366real(dp),               intent(inout) :: depth_filling    ! Depth where pore filling begins [m]
367integer(di),            intent(out)   :: index_filling    ! Index where the pore filling begins [1]
368integer(di),            intent(out)   :: index_geothermal ! Index where the ice table stops [1]
369real(dp),               intent(out)   :: depth_geothermal ! Depth where the ice table stops [m]
370real(dp), dimension(:), intent(out)   :: dz_etadz_rho     ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase
371
372! LOCAL VARIABLES
373! ---------------
374real(dp), dimension(nsoil) :: eta                          ! Constriction
375real(dp), dimension(nsoil) :: dz_psat                      ! First derivative of the vapor pressure at saturation
376real(dp), dimension(nsoil) :: dz_eta                       ! \partial z \eta
377real(dp), dimension(nsoil) :: dzz_psat                     ! \partial \partial psat
378integer(di)                :: ilay, index_tmp              ! Index for loop
379real(dp)                   :: old_depth_filling            ! Depth_filling saved
380real(dp)                   :: Jdry                         ! Flux trought the dry layer
381real(dp)                   :: Jsat                         ! Flux trought the ice layer
382real(dp)                   :: Jdry_prevlay, Jsat_prevlay   ! Same but for the previous ice layer
383integer(di)                :: index_firstice               ! First index where ice appears (i.e., f > 0)
384real(dp)                   :: massfillabove, massfillafter ! H2O mass above and after index_geothermal
385real(dp), parameter        :: pvap2rho = 18.e-3/8.314
386
387! CODE
388! ----
389! 0. Compute constriction over the layer
390! Within the ice sheet, constriction is set to 0. Elsewhere, constriction =  (1-porefilling)**2
391if (index_IS < 0) then
392    index_tmp = nsoil
393    do ilay = 1,nsoil
394        eta(ilay) = constriction(porefill(ilay))
395    end do
396else
397    index_tmp = index_IS
398    do ilay = 1,index_IS - 1
399        eta(ilay) = constriction(porefill(ilay))
400    end do
401    do ilay = index_IS,nsoil
402        eta(ilay) = 0._dp
403    end do
404end if
405
406! 1. Depth at which pore filling occurs. We solve Eq. 9 from  Schorgofer, Icarus  (2010)
407old_depth_filling = depth_filling
408
409call deriv1(mlayer,nsoil,psat_soil,psat_surf,psat_bottom,dz_psat)
410
411do ilay = 1,index_tmp
412    Jdry = (psat_soil(ilay) - pwat_surf)/mlayer(ilay) ! left member of Eq. 9 from Schorgofer, Icarus  (2010)
413    Jsat =  eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010)
414    if (Jdry - Jsat <= 0) then
415        index_filling = ilay
416        exit
417    end if
418end do
419
420if (index_filling == 1) then
421    depth_filling = mlayer(1)
422else if (index_filling > 1) then
423    Jdry_prevlay = (psat_soil(index_filling - 1) - pwat_surf)/mlayer(index_filling - 1)
424    Jsat_prevlay = eta(index_filling - 1)*dz_psat(index_filling - 1)
425    call findroot(Jdry - Jsat,Jdry_prevlay - Jsat_prevlay,mlayer(index_filling),mlayer(index_filling - 1),depth_filling)
426end if
427
428! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010))
429! 2.0 preliminary: depth to shallowest ice (discontinuity at interface)
430index_firstice = -1
431do ilay = 1,nsoil
432    if (porefill(ilay) <= 0._dp) then
433        index_firstice = ilay  ! first point with ice
434        exit
435    end if
436end do
437
438! 2.1: now we can compute
439call deriv1(mlayer,nsoil,eta,1.,eta(nsoil - 1),dz_eta)
440
441if (index_firstice > 0 .and. index_firstice < nsoil - 2) call deriv1_onesided(index_firstice,mlayer,nsoil,eta,dz_eta(index_firstice))
442
443call deriv2_simple(mlayer,nsoil,psat_soil,psat_surf,psat_bottom,dzz_psat)
444dz_etadz_rho = pvap2rho*(dz_eta*dz_psat + eta*dzz_psat)
445
446! 3. Ice table boundary due to geothermal heating
447if (index_IS > 0) index_geothermal = -1
448if (index_geothermal < 0) depth_geothermal = -1._dp
449if ((index_geothermal > 0).and.(index_IS < 0)) then ! Eq. 21 from Schorfoger, Icarus (2010)
450    index_geothermal = -1
451    do ilay = 2,nsoil
452        if (dz_psat(ilay) > 0._dp) then  ! first point with reversed flux
453            index_geothermal = ilay
454            call findroot(dz_psat(ilay - 1),dz_psat(ilay),mlayer(ilay - 1),mlayer(ilay),depth_geothermal)
455            exit
456        end if
457    end do
458else
459    index_geothermal = -1
460end if
461if (index_geothermal > 0 .and. index_IS < 0) then ! Eq. 24 from Schorgofer, Icarus (2010)
462    call  colint(porefill/eta,mlayer,nsoil,index_geothermal - 1,nsoil,massfillabove)
463    index_tmp = -1
464    do ilay = index_geothermal,nsoil
465        if (minval(eta(ilay:nsoil)) <= 0._dp) cycle ! eta=0 means completely full
466        call colint(porefill/eta,mlayer,nsoil,ilay,nsoil,massfillafter)
467        if (massfillafter < dz_psat(ilay)*pvap2rho*B) then ! usually executes on i=typeG
468            if (ilay > index_geothermal) then
469                call findroot(dz_psat(ilay - 1)*pvap2rho*B - massfillabove, &
470                              dz_psat(ilay)*pvap2rho*B - massfillafter,mlayer(ilay - 1),mlayer(ilay),depth_geothermal)
471!                if (depth_geothermal > mlayer(ilay) .or. depth_geothermal < mlayer(ilay - 1)) write(*,*) '# WARNING: zdepthG interpolation failed',ilay,mlayer(ilay - 1),depth_geothermal,mlayer(ilay)
472              index_tmp = ilay
473           end if
474           ! otherwise leave depth_geothermal unchanged
475           exit
476        end if
477        massfillabove = massfillafter
478    end do
479    if (index_tmp > 0) index_geothermal = index_tmp
480end if
481
482END SUBROUTINE find_layering_icetable
483!=======================================================================
484
485!=======================================================================
486SUBROUTINE compute_Tice(ptsoil,ptsurf,ice_depth,Tice)
487!-----------------------------------------------------------------------
488! NAME
489!     compute_Tice
490!
491! DESCRIPTION
492!     Compute subsurface ice temperature by interpolating the
493!     temperatures between the two adjacent cells.
494!
495! AUTHORS & DATE
496!     L. Lange
497!     JB Clement, 12/12/2025
498!
499! NOTES
500!
501!-----------------------------------------------------------------------
502
503! DEPENDENCIES
504! ------------
505use geometry, only: nsoil
506use soil,     only: mlayer
507use stoppage, only: stop_clean
508
509! DECLARATION
510! -----------
511implicit none
512
513! ARGUMENTS
514! ---------
515real(dp), dimension(:), intent(in)  :: ptsoil    ! Soil temperature (K)
516real(dp),               intent(in)  :: ptsurf    ! Surface temperature (K)
517real(dp),               intent(in)  :: ice_depth ! Ice depth (m)
518real(dp),               intent(out) :: Tice      ! Ice temperatures (K)
519
520! LOCAL VARIABLES
521! ---------------
522integer(di) :: ik       ! Loop variables
523integer(di) :: indexice ! Index of the ice
524
525! CODE
526! ----
527indexice = -1
528if (ice_depth >= mlayer(nsoil - 1)) then
529    Tice = ptsoil(nsoil)
530else
531    if(ice_depth < mlayer(0)) then
532        indexice = 0
533    else
534        do ik = 0,nsoil - 2 ! go through all the layers to find the ice locations
535            if (mlayer(ik) <= ice_depth .and. mlayer(ik + 1) > ice_depth) then
536                indexice = ik + 1
537                exit
538            end if
539        end do
540    end if
541    if (indexice < 0) then
542        call stop_clean(__FILE__,__LINE__,"subsurface ice is below the last soil layer!",1)
543    else
544        if(indexice >= 1) then ! Linear inteprolation between soil temperature
545            Tice = (ptsoil(indexice) - ptsoil(indexice + 1))/(mlayer(indexice - 1) - mlayer(indexice))*(ice_depth - mlayer(indexice)) + ptsoil(indexice + 1)
546        else ! Linear inteprolation between the 1st soil temperature and the surface temperature
547            Tice = (ptsoil(1) - ptsurf)/mlayer(0)*(ice_depth - mlayer(0)) + ptsoil(1)
548        end if ! index ice >= 1
549    end if !indexice < 0
550end if ! icedepth > mlayer(nsoil - 1)
551
552END SUBROUTINE compute_Tice
553!=======================================================================
554
555!=======================================================================
556FUNCTION rho_ice(T,ice) RESULT(rho)
557!-----------------------------------------------------------------------
558! NAME
559!     rho_ice
560!
561! DESCRIPTION
562!     Compute subsurface ice density.
563!
564! AUTHORS & DATE
565!     JB Clement, 12/12/2025
566!
567! NOTES
568!
569!-----------------------------------------------------------------------
570
571! DEPENDENCIES
572! ------------
573use stoppage, only: stop_clean
574
575! DECLARATION
576! -----------
577implicit none
578
579! ARGUMENTS
580! ---------
581real(dp),     intent(in) :: T
582character(3), intent(in) :: ice
583
584! LOCAL VARIABLES
585! ---------------
586real(dp) :: rho
587
588! CODE
589! ----
590select case (trim(adjustl(ice)))
591    case('H2O')
592        rho = -3.5353e-4_dp*T**2 + 0.0351_dp*T + 933.5030_dp ! Rottgers 2012
593    case('CO2')
594        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
595    case default
596        call stop_clean(__FILE__,__LINE__,"type of ice unknown!",1)
597end select
598
599END FUNCTION rho_ice
600!=======================================================================
601
602!=======================================================================
603SUBROUTINE evolve_ice_table(h2o_ice,h2o_surfdensity_avg,h2o_soildensity_avg,tsoil_avg,tsurf_avg,delta_icetable, &
604                            q_h2o_ts,ps_avg,icetable_depth,icetable_thickness,ice_porefilling,icetable_depth_old)
605!-----------------------------------------------------------------------
606! NAME
607!     evolve_ice_table
608!
609! DESCRIPTION
610!     Update the subsurface ice table using either the equilibrium or
611!     dynamic method and calculate the resulting H2O mass exchange.
612!
613! AUTHORS & DATE
614!     JB Clement, 12/2025
615!     C. Metz, 02/2026
616!-----------------------------------------------------------------------
617
618! DEPENDENCIES
619! ------------
620use geometry,       only: ngrid, nsoil, nslope, nday
621use slopes,         only: subslope_dist, def_slope_mean
622use surf_ice,       only: threshold_h2oice_cap
623use maths,          only: pi
624use soil,           only: TI
625use subsurface_ice, only: dyn_ss_ice_m
626use display,        only: print_msg
627
628! DECLARATION
629! -----------
630implicit none
631
632! ARGUMENTS
633! ---------
634real(dp), dimension(:),     intent(in)    :: ps_avg
635real(dp), dimension(:,:),   intent(in)    :: h2o_ice, h2o_surfdensity_avg, tsurf_avg, q_h2o_ts
636real(dp), dimension(:,:,:), intent(in)    :: h2o_soildensity_avg, tsoil_avg
637real(dp), dimension(:,:),   intent(out)   :: icetable_depth_old
638real(dp), dimension(:),     intent(out)   :: delta_icetable
639real(dp), dimension(:,:),   intent(inout) :: icetable_depth, icetable_thickness
640real(dp), dimension(:,:,:), intent(inout) :: ice_porefilling
641
642! LOCAL VARIABLES
643! ---------------
644integer(di)                                :: i, islope
645logical(k4), dimension(ngrid)              :: is_h2o_perice
646real(dp),    dimension(ngrid,nslope)       :: icetable_thickness_old
647real(dp),    dimension(ngrid,nsoil,nslope) :: ice_porefilling_old
648
649! CODE
650! ----
651do i = 1,ngrid
652    if (sum(h2o_ice(i,:)*subslope_dist(i,:)/cos(pi*def_slope_mean(:)/180._dp)) > threshold_h2oice_cap) then
653        ! There is enough ice to be considered as an infinite reservoir
654        is_h2o_perice(i) = .true.
655    else
656        ! Too little ice to be considered as an infinite reservoir so ice is transferred to the frost
657        is_h2o_perice(i) = .false.
658    end if
659end do
660if (icetable_equilibrium) then
661    call print_msg("> Evolution of ice table (equilibrium method)")
662    icetable_thickness_old(:,:) = icetable_thickness(:,:)
663    call computeice_table_equilibrium(is_h2o_perice,h2o_surfdensity_avg,h2o_soildensity_avg,TI(:,1,:),icetable_depth,icetable_thickness)
664else if (icetable_dynamic) then
665    call print_msg("> Evolution of ice table (dynamic method)")
666    ice_porefilling_old(:,:,:) = ice_porefilling(:,:,:)
667    icetable_depth_old(:,:) = icetable_depth(:,:)
668    do i = 1,ngrid
669        do islope = 1,nslope
670            call dyn_ss_ice_m(icetable_depth(i,islope),tsurf_avg(i,islope), tsoil_avg(i,:,islope),nsoil,TI(i,1,islope),ps_avg, &
671                              (/sum(q_h2o_ts(i,:))/real(nday,dp)/),ice_porefilling(i,:,islope),ice_porefilling(i,:,islope),icetable_depth(i,islope))
672        end do
673    end do
674end if
675
676! Mass of H2O exchange between the ssi and the atmosphere
677call compute_deltam_icetable(icetable_thickness_old,ice_porefilling_old,tsurf_avg,tsoil_avg,delta_icetable)
678
679END SUBROUTINE evolve_ice_table
680!=======================================================================
681
682END MODULE ice_table
Note: See TracBrowser for help on using the repository browser.