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

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

PEM:

  • Correction of the variable name for the ice table depth in "pemetat0.F90". So it is now got as intended from the "startpem.nc" file;
  • Renaming of the tendencies in the PEM with the prefix 'd_' instead of 'tend_';
  • Modification of the PEM time step type from integer to real. As a consequence, all time variables are now of real type. This change adds the possibility to consider fractions of year as time step.

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