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

Last change on this file since 3512 was 3512, checked in by jbclement, 9 days ago

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