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

Last change on this file since 4066 was 4065, checked in by jbclement, 2 weeks ago

PEM:
Major refactor following the previous ones (r3989 and r3991) completing the large structural reorganization and cleanup of the PEM codebase. This revision introduces newly designed modules, standardizes interfaces with explicit ini/end APIs and adds native NetCDF I/O together with explicit PCM/PEM adapters. In detail:

  • Some PEM models were corrected or improved:
    • Frost/perennial ice semantics are clarified via renaming;
    • Soil temperature remapping clarified, notably by removing the rescaling of temperature deviation;
    • Geothermal flux for the PCM is computed based on the PEM state;
  • New explicit PEM/PCM adapters ("set_*"/"build4PCM_*") to decouple PEM internal representation from PCM file layouts and reconstruct consistent fields returned to the PCM;
  • New explicit build/teardown routines that centralize allocation and initialization ordering, reducing accidental use of uninitialized data and making the model lifecycle explicit;
  • Add native read/write helpers for NetCDF that centralize all low-level NetCDF interactions with major improvements (and more simplicity) compared to legacy PEM/PCM I/O (see the modules "io_netcdf" and "output"). They support reading, creation and writing of "diagevol.nc" (renamed from "diagpem.nc") and start/restart files;
  • Provide well-focused modules ("numerics"/"maths"/"utility"/"display") to host commonly-used primitives:
    • "numerics" defines numerical types and constants for reproducibility, portability across compilers and future transitions (e.g. quadruple precision experiments);
    • "display" provides a single controlled interface for runtime messages, status output and diagnostics, avoiding direct 'print'/'write' to enable silent mode, log redirection, and MPI-safe output in the future.
    • "utility" (new module) hosts generic helpers used throughout the code (e.g. "int2str" or "real2str");
  • Add modules "clim_state_init"/"clim_state_rec" which provide robust read/write logic for "start/startfi/startpem", including 1D fallbacks, mesh conversions and dimension checks. NetCDF file creation is centralized and explicit. Restart files are now self-consistent and future-proof, requiring changes only to affected variables;
  • Add module "atmosphere" which computes pressure fields, reconstructs potential temperature and air mass. It also holds the whole logic to define sigma or hybrid coordinates for altitudes;
  • Add module "geometry" to centrilize dimensions logic and grid conversions routines (including 2 new ones "dyngrd2vect"/"vect2dyngrd");
  • Add module "slopes" to isolate slopes handling;
  • Add module "surface" to isolate surface management. Notably, albedo and emissivity are now fully reconstructed following the PCM settings;
  • Add module "allocation" to check the dimension initialization and centrilize allocation/deallocation;
  • Finalize module decomposition and renaming to consolidate domain-based modules, purpose-based routines and physics/process-based variables;
  • The main program now drives a clearer sequence of conceptual steps (initialization / reading / evolution / update / build / writing) and fails explicitly instead of silently defaulting;
  • Ice table logic is made restart-safe;
  • 'Open'/'read' intrinsic logic is made safe and automatic;
  • Improve discoverability and standardize the data handling (private vs protected vs public);
  • Apply consistent documentation/header style already introduced;
  • Update deftank scripts to reflect new names and launch wrappers;

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