source: trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90 @ 3527

Last change on this file since 3527 was 3527, checked in by jbclement, 26 hours ago

PEM:

  • Removing redundant Norbert Schorghofer's subroutines/parameters (follow-up of r3526);
  • Making all Norbert Schorghofer's subroutines with modern explicit interface via modules;
  • Cleaning of "glaciers_mod.F90";
  • Optimization for the computation of ice density according to temperature by using a function.

JBC
File size: 20.8 KB
Line 
1MODULE ice_table_mod
2
3implicit none
4
5!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
6!!!
7!!! Purpose: Ice table (pore-filling) variables and methods to compute it (dynamic and static)
8!!! Author:  LL, 02/2023
9!!!
10!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
11
12logical                             :: icetable_equilibrium ! Boolean to say if the PEM needs to recompute the icetable depth when at equilibrium
13logical                             :: icetable_dynamic     ! Boolean to say if the PEM needs to recompute the icetable depth with the dynamic method
14real, allocatable, dimension(:,:)   :: icetable_depth       ! ngrid x nslope: Depth of the ice table [m]
15real, allocatable, dimension(:,:)   :: icetable_thickness   ! ngrid x nslope: Thickness of the ice table [m]
16real, allocatable, dimension(:,:,:) :: ice_porefilling      ! the amout of porefilling in each layer in each grid [m^3/m^3]
17
18!-----------------------------------------------------------------------
19contains
20
21!-----------------------------------------------------------------------
22SUBROUTINE ini_ice_table(ngrid,nslope,nsoil)
23
24implicit none
25
26integer, intent(in) :: ngrid  ! number of atmospheric columns
27integer, intent(in) :: nslope ! number of slope within a mesh
28integer, intent(in) :: nsoil  ! number of soil layers
29
30allocate(icetable_depth(ngrid,nslope))
31allocate(icetable_thickness(ngrid,nslope))
32allocate(ice_porefilling(ngrid,nsoil,nslope))
33
34END SUBROUTINE ini_ice_table
35
36!-----------------------------------------------------------------------
37SUBROUTINE end_ice_table
38
39implicit none
40
41if (allocated(icetable_depth)) deallocate(icetable_depth)
42if (allocated(icetable_thickness)) deallocate(icetable_thickness)
43if (allocated(ice_porefilling)) deallocate(ice_porefilling)
44
45END SUBROUTINE end_ice_table
46
47!------------------------------------------------------------------------------------------------------
48SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness)
49!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
50!!!
51!!! Purpose: Compute the ice table depth (pore-filling) knowing the yearly average water
52!!! density at the surface and at depth.
53!!! Computations are made following the methods in Schorgofer et al., 2005
54!!! This subroutine only gives the ice table at equilibrium and does not consider exchange with the atmosphere
55!!!
56!!! Author: LL
57!!!
58!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
59use math_mod,                   only: findroot
60use comsoil_h_PEM,              only: mlayer_PEM, layer_PEM ! Depth of the vertical grid
61use soil_thermalproperties_mod, only: ice_thermal_properties
62
63implicit none
64
65! Inputs
66! ------
67integer,                             intent(in) :: ngrid, nslope, nsoil ! Size of the physical grid, number of subslope, number of soil layer in the PEM
68logical, dimension(ngrid),           intent(in) :: watercaptag          ! Boolean to check the presence of a perennial glacier
69real, dimension(ngrid,nslope),       intent(in) :: rhowatersurf_ave     ! Water density at the surface, yearly averaged [kg/m^3]
70real, 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]
71real, dimension(ngrid,nslope),       intent(in) :: regolith_inertia     ! Thermal inertia of the regolith layer [SI]
72! Ouputs
73! ------
74real, dimension(ngrid,nslope), intent(out) :: ice_table_beg       ! ice table depth [m]
75real, dimension(ngrid,nslope), intent(out) :: ice_table_thickness ! ice table thickness [m]
76! Locals
77! ------
78integer                       :: ig, islope, isoil, isoilend ! loop variables
79real, dimension(nsoil)        :: diff_rho                    ! difference of water vapor density between the surface and at depth [kg/m^3]
80real                          :: ice_table_end               ! depth of the end of the ice table  [m]
81real, dimension(ngrid,nslope) :: previous_icetable_depth     ! Ice table computed at previous ice depth [m]
82real                          :: stretch                     ! stretch factor to improve the convergence of the ice table
83real                          :: wice_inertia                ! Water Ice thermal Inertia [USI]
84! Code
85! ----
86previous_icetable_depth = ice_table_beg
87do ig = 1,ngrid
88    if (watercaptag(ig)) then
89        ice_table_beg(ig,:) = 0.
90        ice_table_thickness(ig,:) = layer_PEM(nsoil) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
91    else
92        do islope = 1,nslope
93            ice_table_beg(ig,islope) = -1.
94            ice_table_thickness(ig,islope) = 0.
95            do isoil = 1,nsoil
96                diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope)
97            enddo
98            if (diff_rho(1) > 0) then ! ice is at the surface
99                ice_table_beg(ig,islope) = 0.
100                do isoilend = 2,nsoil - 1
101                    if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend + 1) < 0.) then
102                        call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend),mlayer_PEM(isoilend + 1),ice_table_end)
103                        ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope)
104                        exit
105                    endif
106                enddo
107            else
108                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
109                    if (diff_rho(isoil) < 0 .and. diff_rho(isoil + 1) > 0.) then
110                        call findroot(diff_rho(isoil),diff_rho(isoil + 1),mlayer_PEM(isoil),mlayer_PEM(isoil + 1),ice_table_beg(ig,islope))
111                        ! Now let's find the end of the ice table
112                        ice_table_thickness(ig,islope) = layer_PEM(nsoil) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
113                        do isoilend = isoil + 1,nsoil - 1
114                            if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend + 1) < 0.) then
115                                call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend),mlayer_PEM(isoilend + 1),ice_table_end)
116                                ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope)
117                                exit
118                            endif
119                        enddo
120                    endif ! diff_rho(z) <0 & diff_rho(z+1) > 0
121                enddo
122            endif ! diff_rho(1) > 0
123        enddo
124    endif ! watercaptag
125enddo
126
127! Small trick to speed up the convergence, Oded's idea.
128do islope = 1,nslope
129    do ig = 1,ngrid
130        if (ice_table_beg(ig,islope) > previous_icetable_depth(ig,islope) .and. previous_icetable_depth(ig,islope) >= 0) then
131            call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),wice_inertia)
132            stretch = (regolith_inertia(ig,islope)/wice_inertia)**2
133            ice_table_thickness(ig,islope) = ice_table_thickness(ig,islope) + (ice_table_beg(ig,islope) - &
134                                             previous_icetable_depth(ig,islope) + (ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch)
135            ice_table_beg(ig,islope) = previous_icetable_depth(ig,islope) + (ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch
136        endif
137    enddo
138enddo
139
140END SUBROUTINE computeice_table_equilibrium
141
142!-----------------------------------------------------------------------
143SUBROUTINE compute_massh2o_exchange_ssi(ngrid,nslope,nsoil,icetable_thickness_old,ice_porefilling_old,tsurf,tsoil,delta_m_h2o)
144!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
145!!!
146!!! Purpose: Compute the mass of H2O that has sublimated from the ice table / condensed
147!!!
148!!! Author: LL
149!!!
150!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
151use comsoil_h_PEM,         only: mlayer_PEM
152use comslope_mod,          only: subslope_dist, def_slope_mean
153use constants_marspem_mod, only: porosity
154#ifndef CPP_STD
155    use comcstfi_h,   only: pi
156#else
157    use comcstfi_mod, only: pi
158#endif
159
160implicit none
161
162! Inputs
163! ------
164integer,                             intent(in) :: ngrid, nslope, nsoil   ! Size of the physical grid, number of subslope, number of soil layer in the PEM
165real, dimension(ngrid,nslope),       intent(in) :: icetable_thickness_old ! Ice table thickness at the previous iteration [m]
166real, dimension(ngrid,nsoil,nslope), intent(in) :: ice_porefilling_old    ! Ice pore filling at the previous iteration [m]
167real, dimension(ngrid,nslope),       intent(in) :: tsurf                  ! Surface temperature [K]
168real, dimension(ngrid,nsoil,nslope), intent(in) :: tsoil                  ! Soil temperature [K]
169! Outputs
170! -------
171real, 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]
172! Locals
173!-------
174integer :: ig, islope, isoil ! loop index
175real    :: rho               ! density of water ice [kg/m^3]
176real    :: Tice              ! ice temperature [k]
177! Code
178! ----
179
180! Let's compute the amount of ice that has sublimated in each subslope
181if (icetable_equilibrium) then
182    delta_m_h2o = 0.
183    do ig = 1,ngrid
184        do islope = 1,nslope
185            call compute_Tice_pem(nsoil,tsoil(ig,:,islope),tsurf(ig,islope),icetable_depth(ig,islope),Tice)
186            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
187                              *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
188        enddo
189    enddo
190else if (icetable_dynamic) then
191    delta_m_h2o = 0.
192    do ig = 1,ngrid
193        do islope = 1,nslope
194            do isoil = 1,nsoil
195                call compute_Tice_pem(nsoil,tsoil(ig,:,islope),tsurf(ig,islope),mlayer_PEM(isoil - 1),Tice)
196                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
197                                  *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
198            enddo
199        enddo
200    enddo
201endif
202
203END SUBROUTINE compute_massh2o_exchange_ssi
204
205!-----------------------------------------------------------------------
206SUBROUTINE 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)
207!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
208!!!
209!!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM
210!!!
211!!! Author: LL
212!!!
213!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
214use comsoil_h_PEM, only: nsoilmx_PEM, mlayer_PEM
215use math_mod,      only: deriv1, deriv1_onesided, colint, findroot, deriv2_simple
216
217implicit none
218
219! Inputs
220! ------
221real, dimension(nsoilmx_PEM), intent(in) :: porefill    ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice
222real, dimension(nsoilmx_PEM), intent(in) :: psat_soil   ! Soil water pressure at saturation, yearly averaged [Pa]
223real,                         intent(in) :: psat_surf   ! surface water pressure at saturation, yearly averaged [Pa]
224real,                         intent(in) :: pwat_surf   ! Water vapor pressure  at the surface, not necesseraly at saturation, yearly averaged [Pa]
225real,                         intent(in) :: psat_bottom ! Boundary conditions for soil vapor pressure [Pa]
226real,                         intent(in) :: B           ! constant (Eq. 8 from  Schorgofer, Icarus (2010).)
227integer,                      intent(in) :: index_IS    ! index of the soil layer where the ice sheet begins [1]
228real, intent(inout) :: depth_filling ! depth where pore filling begins [m]
229! Outputs
230! -------
231integer,                      intent(out) :: index_filling    ! index where the pore filling begins [1]
232integer,                      intent(out) :: index_geothermal ! index where the ice table stops [1]
233real,                         intent(out) :: depth_geothermal ! depth where the ice table stops [m]
234real, dimension(nsoilmx_PEM), intent(out) :: dz_etadz_rho     ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase
235! Locals
236!-------
237real, dimension(nsoilmx_PEM) :: eta                          ! constriction
238real, dimension(nsoilmx_PEM) :: dz_psat                      ! first derivative of the vapor pressure at saturation
239real, dimension(nsoilmx_PEM) :: dz_eta                       ! \partial z \eta
240real, dimension(nsoilmx_PEM) :: dzz_psat                     ! \partial \partial psat
241integer                      :: ilay, index_tmp              ! index for loop
242real                         :: old_depth_filling            ! depth_filling saved
243real                         :: Jdry                         ! flux trought the dry layer
244real                         :: Jsat                         ! flux trought the ice layer
245real                         :: Jdry_prevlay, Jsat_prevlay   ! same but for the previous ice layer
246integer                      :: index_firstice               ! first index where ice appears (i.e., f > 0)
247real                         :: massfillabove, massfillafter ! h2O mass above and after index_geothermal
248! Constant
249!---------
250real, parameter :: pvap2rho = 18.e-3/8.314
251! Code
252!-----
253! 0. Compute constriction over the layer
254! Within the ice sheet, constriction is set to 0. Elsewhere, constriction =  (1-porefilling)**2
255if (index_IS < 0) then
256    index_tmp = nsoilmx_PEM
257    do ilay = 1,nsoilmx_PEM
258        eta(ilay) = constriction(porefill(ilay))
259    enddo
260else
261    index_tmp = index_IS
262    do ilay = 1,index_IS - 1
263        eta(ilay) = constriction(porefill(ilay))
264    enddo
265    do ilay = index_IS,nsoilmx_PEM
266        eta(ilay) = 0.
267    enddo
268endif
269
270! 1. Depth at which pore filling occurs. We solve Eq. 9 from  Schorgofer, Icarus  (2010)
271old_depth_filling = depth_filling
272
273call deriv1(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dz_psat)
274
275do ilay = 1,index_tmp
276    Jdry = (psat_soil(ilay) - pwat_surf)/mlayer_PEM(ilay) ! left member of Eq. 9 from Schorgofer, Icarus  (2010)
277    Jsat =  eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010)
278    if (Jdry - Jsat <= 0) then
279        index_filling = ilay
280        exit
281    endif
282enddo
283
284if (index_filling == 1) then
285    depth_filling = mlayer_PEM(1)
286else if (index_filling > 1) then
287    Jdry_prevlay = (psat_soil(index_filling - 1) - pwat_surf)/mlayer_PEM(index_filling - 1)
288    Jsat_prevlay = eta(index_filling - 1)*dz_psat(index_filling - 1)
289    call findroot(Jdry - Jsat,Jdry_prevlay - Jsat_prevlay,mlayer_PEM(index_filling),mlayer_PEM(index_filling - 1),depth_filling)
290endif
291
292! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010))
293! 2.0 preliminary: depth to shallowest ice (discontinuity at interface)
294index_firstice = -1
295do ilay = 1,nsoilmx_PEM
296    if (porefill(ilay) <= 0.) then
297        index_firstice = ilay  ! first point with ice
298        exit
299    endif
300enddo
301
302! 2.1: now we can compute
303call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM - 1),dz_eta)
304
305if (index_firstice > 0 .and. index_firstice < nsoilmx_PEM - 2) call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice))
306
307call deriv2_simple(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dzz_psat)
308dz_etadz_rho = pvap2rho*(dz_eta*dz_psat + eta*dzz_psat)
309
310! 3. Ice table boundary due to geothermal heating
311if (index_IS > 0) index_geothermal = -1
312if (index_geothermal < 0) depth_geothermal = -1.
313if ((index_geothermal > 0).and.(index_IS < 0)) then ! Eq. 21 from Schorfoger, Icarus (2010)
314    index_geothermal = -1
315    do ilay = 2,nsoilmx_PEM
316        if (dz_psat(ilay) > 0.) then  ! first point with reversed flux
317            index_geothermal = ilay
318            call findroot(dz_psat(ilay - 1),dz_psat(ilay),mlayer_PEM(ilay - 1),mlayer_PEM(ilay),depth_geothermal)
319            exit
320        endif
321    enddo
322else
323    index_geothermal = -1
324endif
325if (index_geothermal > 0 .and. index_IS < 0) then ! Eq. 24 from Schorgofer, Icarus (2010)
326    call  colint(porefill/eta,mlayer_PEM,nsoilmx_PEM,index_geothermal - 1,nsoilmx_PEM,massfillabove)
327    index_tmp = -1
328    do ilay = index_geothermal,nsoilmx_PEM
329        if (minval(eta(ilay:nsoilmx_PEM)).le.0.) cycle ! eta=0 means completely full
330        call colint(porefill/eta,mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter)
331        if (massfillafter < dz_psat(ilay)*pvap2rho*B) then ! usually executes on i=typeG
332            if (ilay > index_geothermal) then
333!                write(34,*) '# adjustment to geotherm depth by',ilay-index_geothermal
334                call findroot(dz_psat(ilay - 1)*pvap2rho*B - massfillabove, &
335                              dz_psat(ilay)*pvap2rho*B - massfillafter,mlayer_PEM(ilay - 1),mlayer_PEM(ilay),depth_geothermal)
336!                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)
337              index_tmp = ilay
338           endif
339           ! otherwise leave depth_geothermal unchanged
340           exit
341        endif
342        massfillabove = massfillafter
343    enddo
344    if (index_tmp > 0) index_geothermal = index_tmp
345end if
346
347END SUBROUTINE find_layering_icetable
348
349!-----------------------------------------------------------------------
350FUNCTION constriction(porefill) RESULT(eta)
351
352!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
353!!!
354!!! Purpose: Compute the constriction of vapor flux by pore ice
355!!!
356!!! Author: LL
357!!!
358!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
359
360implicit none
361
362real, intent(in) :: porefill ! pore filling fraction
363real :: eta ! constriction
364
365if (porefill <= 0.) then
366    eta = 1.
367else if (0 < porefill .and. porefill < 1.) then
368    eta = (1-porefill)**2 ! Hudson et al., JGR, 2009
369else
370    eta = 0.
371endif
372
373END FUNCTION constriction
374
375!-----------------------------------------------------------------------
376SUBROUTINE compute_Tice_pem(nsoil,ptsoil,ptsurf,ice_depth,Tice)
377!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
378!!!
379!!! Purpose: Compute subsurface ice temperature by interpolating the temperatures between the two adjacent cells.
380!!!
381!!! Author: LL
382!!!
383!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
384
385use comsoil_h_PEM, only: layer_PEM, mlayer_PEM
386use abort_pem_mod, only: abort_pem
387
388implicit none
389
390! Inputs
391! ------
392integer,                intent(in) :: nsoil     ! Number of soil layers
393real, dimension(nsoil), intent(in) :: ptsoil    ! Soil temperature (K)
394real,                   intent(in) :: ptsurf    ! Surface temperature (K)
395real,                   intent(in) :: ice_depth ! Ice depth (m)
396
397! Outputs
398! -------
399real, intent(out) :: Tice ! Ice temperatures (K)
400
401! Locals
402! ------
403integer :: ik       ! Loop variables
404integer :: indexice ! Index of the ice
405
406! Code
407!-----
408indexice = -1
409if (ice_depth >= mlayer_PEM(nsoil - 1)) then
410    Tice = ptsoil(nsoil)
411else
412    if(ice_depth < mlayer_PEM(0)) then
413        indexice = 0.
414    else     
415        do ik = 0,nsoil - 2 ! go through all the layers to find the ice locations
416            if (mlayer_PEM(ik) <= ice_depth .and. mlayer_PEM(ik + 1) > ice_depth) then
417                indexice = ik + 1
418                exit
419            endif
420        enddo
421    endif
422    if (indexice < 0) then
423        call abort_pem("compute_Tice_pem","Subsurface ice is below the last soil layer!",1)
424    else
425        if(indexice >= 1) then ! Linear inteprolation between soil temperature
426            Tice = (ptsoil(indexice) - ptsoil(indexice + 1))/(mlayer_PEM(indexice - 1) - mlayer_PEM(indexice))*(ice_depth - mlayer_PEM(indexice)) + ptsoil(indexice + 1)
427        else ! Linear inteprolation between the 1st soil temperature and the surface temperature
428            Tice = (ptsoil(1) - ptsurf)/mlayer_PEM(0)*(ice_depth - mlayer_PEM(0)) + ptsoil(1)
429        endif ! index ice >= 1
430    endif !indexice < 0
431endif ! icedepth > mlayer_PEM(nsoil - 1)
432
433END SUBROUTINE compute_Tice_pem
434
435!-----------------------------------------------------------------------
436FUNCTION rho_ice(T,ice) RESULT(rho)
437!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
438!!!
439!!! Purpose: Compute subsurface ice density
440!!!
441!!! Author: JBC
442!!!
443!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
444
445use abort_pem_mod, only: abort_pem
446
447implicit none
448
449real,         intent(in) :: T
450character(3), intent(in) :: ice
451real :: rho
452
453select case (trim(adjustl(ice)))
454    case('h2o')
455        rho = -3.5353e-4*T**2 + 0.0351*T + 933.5030 ! Rottgers 2012
456    case('co2')
457        rho = (1.72391 - 2.53e-4*T-2.87*1e-7*T**2)*1.e3 ! Mangan et al. 2017
458    case default
459        call abort_pem("rho_ice","Type of ice not known!",1)
460end select
461
462END FUNCTION rho_ice
463
464END MODULE ice_table_mod
Note: See TracBrowser for help on using the repository browser.