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

Last change on this file was 3991, checked in by jbclement, 4 weeks ago

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