New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

icethd_pnd.F90 in NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icethd_pnd.F90 @ 13879

Last change on this file since 13879 was 13879, checked in by vancop, 4 years ago
More decent computing cost
Property svn:keywords set to `Id`
File size: 59.4 KB

Line
1	MODULE icethd_pnd
2	!!======================================================================
3	!! --- closest to CICE version
4	!! * MODULE icethd_pnd *
5	!! sea-ice: Melt ponds on top of sea ice
6	!!======================================================================
7	!! history : ! 2012 (O. Lecomte) Adaptation from Flocco and Turner
8	!! ! 2017 (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation
9	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
10	!!----------------------------------------------------------------------
11	#if defined key_si3
12	!!----------------------------------------------------------------------
13	!! 'key_si3' : SI3 sea-ice model
14	!!----------------------------------------------------------------------
15	!! ice_thd_pnd_init : some initialization and namelist read
16	!! ice_thd_pnd : main calling routine
17	!!----------------------------------------------------------------------
18	USE phycst ! physical constants
19	USE dom_oce ! ocean space and time domain
20	USE ice ! sea-ice: variables
21	USE ice1D ! sea-ice: thermodynamics variables
22	USE icetab ! sea-ice: 1D <==> 2D transformation
23	USE sbc_ice ! surface energy budget
24	!
25	USE in_out_manager ! I/O manager
26	USE lib_mpp ! MPP library
27	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
28	USE timing ! Timing
29
30	IMPLICIT NONE
31	PRIVATE
32
33	PUBLIC ice_thd_pnd_init ! routine called by icestp.F90
34	PUBLIC ice_thd_pnd ! routine called by icestp.F90
35
36	INTEGER :: nice_pnd ! choice of the type of pond scheme
37	! ! associated indices:
38	INTEGER, PARAMETER :: np_pndNO = 0 ! No pond scheme
39	INTEGER, PARAMETER :: np_pndCST = 1 ! Constant ice pond scheme
40	INTEGER, PARAMETER :: np_pndLEV = 2 ! Level ice pond scheme
41	INTEGER, PARAMETER :: np_pndTOPO = 3 ! Level ice pond scheme
42
43	REAL(wp), PARAMETER :: & ! shared parameters for topographic melt ponds
44	zhi_min = 0.1_wp , & ! minimum ice thickness with ponds (m)
45	zTd = 0.15_wp , & ! temperature difference for freeze-up (C)
46	zvp_min = 1.e-4_wp ! minimum pond volume (m)
47
48	!! * Substitutions
49	# include "vectopt_loop_substitute.h90"
50	!!----------------------------------------------------------------------
51	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
52	!! $Id$
53	!! Software governed by the CeCILL license (see ./LICENSE)
54	!!----------------------------------------------------------------------
55	CONTAINS
56
57	SUBROUTINE ice_thd_pnd
58	!!-------------------------------------------------------------------
59	!! * ROUTINE ice_thd_pnd *
60	!!
61	!! ** Purpose : change melt pond fraction and thickness
62	!!
63	!! Note: Melt ponds affect only radiative transfer for now
64	!!
65	!! No heat, no salt.
66	!! The melt water they carry is collected but
67	!! not removed from fw budget or released to the ocean
68	!!
69	!! A wfx_pnd has been coded for diagnostic purposes
70	!! It is not fully consistent yet.
71	!!
72	!! The current diagnostic lacks a contribution from drainage
73	!!
74	!! Each time wfx_pnd is updated, wfx_sum / wfx_snw_sum must be updated
75	!!
76	!!-------------------------------------------------------------------
77	!
78	SELECT CASE ( nice_pnd )
79	!
80	CASE (np_pndCST) ; CALL pnd_CST !== Constant melt ponds ==!
81	!
82	CASE (np_pndLEV) ; CALL pnd_LEV !== Level ice melt ponds ==!
83	!
84	CASE (np_pndTOPO) ; CALL pnd_TOPO !& !== Topographic melt ponds ==!
85	!
86	END SELECT
87	!
88	END SUBROUTINE ice_thd_pnd
89
90
91	SUBROUTINE pnd_CST
92	!!-------------------------------------------------------------------
93	!! * ROUTINE pnd_CST *
94	!!
95	!! ** Purpose : Compute melt pond evolution
96	!!
97	!! ** Method : Melt pond fraction and thickness are prescribed
98	!! to non-zero values when t_su = 0C
99	!!
100	!! ** Tunable parameters : Pond fraction (rn_apnd) & depth (rn_hpnd)
101	!!
102	!! ** Note : Coupling with such melt ponds is only radiative
103	!! Advection, ridging, rafting... are bypassed
104	!!
105	!! ** References : Bush, G.W., and Trump, D.J. (2017)
106	!!-------------------------------------------------------------------
107	INTEGER :: ji ! loop indices
108	!!-------------------------------------------------------------------
109	DO ji = 1, npti
110	!
111	IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN
112	h_ip_1d(ji) = rn_hpnd
113	a_ip_1d(ji) = rn_apnd * a_i_1d(ji)
114	h_il_1d(ji) = 0._wp ! no pond lids whatsoever
115	ELSE
116	h_ip_1d(ji) = 0._wp
117	a_ip_1d(ji) = 0._wp
118	h_il_1d(ji) = 0._wp
119	ENDIF
120	!
121	END DO
122	!
123	END SUBROUTINE pnd_CST
124
125
126	SUBROUTINE pnd_LEV
127	!!-------------------------------------------------------------------
128	!! * ROUTINE pnd_LEV *
129	!!
130	!! ** Purpose : Compute melt pond evolution
131	!!
132	!! ** Method : A fraction of meltwater is accumulated in ponds and sent to ocean when surface is freezing
133	!! We work with volumes and then redistribute changes into thickness and concentration
134	!! assuming linear relationship between the two.
135	!!
136	!! ** Action : - pond growth: Vp = Vp + dVmelt --- from Holland et al 2012 ---
137	!! dVmelt = (1-r)/rhow * ( rhoidh_i + rhosdh_s ) * a_i
138	!! dh_i = meltwater from ice surface melt
139	!! dh_s = meltwater from snow melt
140	!! (1-r) = fraction of melt water that is not flushed
141	!!
142	!! - limtations: a_ip must not exceed (1-r)*a_i
143	!! h_ip must not exceed 0.5*h_i
144	!!
145	!! - pond shrinking:
146	!! if lids: Vp = Vp -dH * a_ip
147	!! dH = lid thickness change. Retrieved from this eq.: --- from Flocco et al 2010 ---
148	!!
149	!! rhoi * Lf * dH/dt = ki * MAX(Tp-Tsu,0) / H
150	!! H = lid thickness
151	!! Lf = latent heat of fusion
152	!! Tp = -2C
153	!!
154	!! And solved implicitely as:
155	!! H(t+dt)*2 -H(t) H(t+dt) -ki * (Tp-Tsu) * dt / (rhoi*Lf) = 0
156	!!
157	!! if no lids: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp) --- from Holland et al 2012 ---
158	!!
159	!! - Flushing: w = -perm/visc * rho_oce * grav * Hp / Hi --- from Flocco et al 2007 ---
160	!! perm = permability of sea-ice
161	!! visc = water viscosity
162	!! Hp = height of top of the pond above sea-level
163	!! Hi = ice thickness thru which there is flushing
164	!!
165	!! - Corrections: remove melt ponds when lid thickness is 10 times the pond thickness
166	!!
167	!! - pond thickness and area is retrieved from pond volume assuming a linear relationship between h_ip and a_ip:
168	!! a_ip/a_i = a_ip_frac = h_ip / zaspect
169	!!
170	!! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min
171	!!
172	!! ** Note : mostly stolen from CICE but not only
173	!!
174	!! ** References : Holland et al (J. Clim, 2012)
175	!!
176	!!-------------------------------------------------------------------
177
178	REAL(wp), DIMENSION(nlay_i) :: ztmp ! temporary array
179	!!
180	REAL(wp), PARAMETER :: zaspect = 0.8_wp ! pond aspect ratio
181	REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature
182	REAL(wp), PARAMETER :: zvisc = 1.79e-3_wp ! water viscosity
183	!!
184	REAL(wp) :: zfr_mlt, zdv_mlt ! fraction and volume of available meltwater retained for melt ponding
185	REAL(wp) :: zdv_frz, zdv_flush ! Amount of melt pond that freezes, flushes
186	REAL(wp) :: zhp ! heigh of top of pond lid wrt ssh
187	REAL(wp) :: zv_ip_max ! max pond volume allowed
188	REAL(wp) :: zdT ! zTp-t_su
189	REAL(wp) :: zsbr ! Brine salinity
190	REAL(wp) :: zperm ! permeability of sea ice
191	REAL(wp) :: zfac, zdum ! temporary arrays
192	REAL(wp) :: z1_rhow, z1_aspect, z1_Tp ! inverse
193	!!
194	INTEGER :: ji, jj, jk, jl ! loop indices
195
196	!!-------------------------------------------------------------------
197
198	z1_rhow = 1._wp / rhow
199	z1_aspect = 1._wp / zaspect
200	z1_Tp = 1._wp / zTp
201
202	!-----------------------------------------------------------------------------------------------
203	! Identify grid cells with ice
204	!-----------------------------------------------------------------------------------------------
205	at_i(:,:) = SUM( a_i, dim=3 )
206	!
207	npti = 0 ; nptidx(:) = 0
208	DO jj = 1, jpj
209	DO ji = 1, jpi
210	IF ( at_i(ji,jj) > epsi10 ) THEN
211	npti = npti + 1
212	nptidx( npti ) = (jj - 1) * jpi + ji
213	ENDIF
214	END DO
215	END DO
216
217	!-----------------------------------------------------------------------------------------------
218	! Prepare 1D arrays
219	!-----------------------------------------------------------------------------------------------
220
221	IF( npti > 0 ) THEN
222
223	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum )
224	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum )
225	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd )
226
227	DO jl = 1, jpl
228
229	CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,jl) )
230	CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,jl) )
231	CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d (1:npti), t_su (:,:,jl) )
232	CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,jl) )
233	CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,jl) )
234	CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,jl) )
235
236	CALL tab_2d_1d( npti, nptidx(1:npti), dh_i_sum (1:npti), dh_i_sum_2d (:,:,jl) )
237	CALL tab_2d_1d( npti, nptidx(1:npti), dh_s_mlt (1:npti), dh_s_mlt_2d (:,:,jl) )
238
239	DO jk = 1, nlay_i
240	CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl) )
241	END DO
242
243	!-----------------------------------------------------------------------------------------------
244	! Go for ponds
245	!-----------------------------------------------------------------------------------------------
246
247
248	DO ji = 1, npti
249	! !----------------------------------------------------!
250	IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN ! Case ice thickness < rn_himin or tiny ice fraction !
251	! !----------------------------------------------------!
252	!--- Remove ponds on thin ice or tiny ice fractions
253	a_ip_1d(ji) = 0._wp
254	h_ip_1d(ji) = 0._wp
255	h_il_1d(ji) = 0._wp
256	! !--------------------------------!
257	ELSE ! Case ice thickness >= rn_himin !
258	! !--------------------------------!
259	v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! retrieve volume from thickness
260	v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji)
261	!
262	!------------------!
263	! case ice melting !
264	!------------------!
265	!
266	!--- available meltwater for melt ponding ---!
267	zdum = -( dh_i_sum(ji)rhoi + dh_s_mlt(ji)rhos ) * z1_rhow * a_i_1d(ji)
268	zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) ! = ( 1 - r ) = fraction of melt water that is not flushed
269	zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors?
270	!
271	!--- overflow ---!
272	! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume
273	! a_ip_max = zfr_mlt * a_i
274	! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:
275	zv_ip_max = zfr_mlt*2 a_i_1d(ji) * zaspect
276	zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
277
278	! If pond depth exceeds half the ice thickness then reduce the pond volume
279	! h_ip_max = 0.5 * h_i
280	! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:
281	zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) ! MV dimensions are wrong here or comment is unclear
282	zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
283
284	!--- Pond growing ---!
285	v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt
286	!
287	!--- Lid melting ---!
288	IF( ln_pnd_lids ) v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0
289	!
290	!--- mass flux ---!
291	! MV I would recommend to remove that
292	! Since melt ponds carry no freshwater there is no point in modifying water fluxes
293
294	IF( zdv_mlt > 0._wp ) THEN
295	zfac = zdv_mlt * rhow * r1_rdtice ! melt pond mass flux < 0 [kg.m-2.s-1]
296	wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
297	!
298	zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) ! adjust ice/snow melting flux > 0 to balance melt pond flux
299	wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum)
300	wfx_sum_1d(ji) = wfx_sum_1d(ji) * (1._wp + zdum)
301	ENDIF
302
303	!-------------------!
304	! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0)
305	!-------------------!
306	!
307	zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp )
308	!
309	!--- Pond contraction (due to refreezing) ---!
310	IF( ln_pnd_lids ) THEN
311	!
312	!--- Lid growing and subsequent pond shrinking ---!
313	zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0
314	& SQRT( v_il_1d(ji)2 + a_ip_1d(ji)2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors
315
316	! Lid growing
317	v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz )
318
319	! Pond shrinking
320	v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz )
321
322	ELSE
323	! Pond shrinking
324	v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6)
325	ENDIF
326	!
327	!--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
328	! v_ip = h_ip * a_ip
329	! a_ip/a_i = a_ip_frac = h_ip / zaspect (cf Holland 2012, fitting SHEBA so that knowing v_ip we can distribute it to a_ip and h_ip)
330	a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i
331	h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji)
332
333	!---------------!
334	! Pond flushing !
335	!---------------!
336	! height of top of the pond above sea-level
337	zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0
338
339	! Calculate the permeability of the ice (Assur 1958, see Flocco 2010)
340	DO jk = 1, nlay_i
341	! MV Assur is inconsistent with SI3
342	zsbr = - 1.2_wp &
343	& - 21.8_wp * ( t_i_1d(ji,jk) - rt0 ) &
344	& - 0.919_wp * ( t_i_1d(ji,jk) - rt0 )**2 &
345	& - 0.0178_wp * ( t_i_1d(ji,jk) - rt0 )**3
346	! MV linear expression more consistent & simpler zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt
347	ztmp(jk) = sz_i_1d(ji,jk) / zsbr
348	END DO
349	zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 )
350
351	! Do the drainage using Darcy's law
352	zdv_flush = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji)
353	zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) )
354	zdv_flush = 0._wp ! MV remove pond drainage for now
355	v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush
356
357	! MV --- why pond drainage does not give back water into freshwater flux ?
358
359	!--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
360	a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i
361	h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji)
362
363	!--- Corrections and lid thickness ---!
364	IF( ln_pnd_lids ) THEN
365	!--- retrieve lid thickness from volume ---!
366	IF( a_ip_1d(ji) > epsi10 ) THEN ; h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji)
367	ELSE ; h_il_1d(ji) = 0._wp
368	ENDIF
369	!--- remove ponds if lids are much larger than ponds ---!
370	IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN
371	a_ip_1d(ji) = 0._wp
372	h_ip_1d(ji) = 0._wp
373	h_il_1d(ji) = 0._wp
374	ENDIF
375	ENDIF
376	!
377	ENDIF
378
379	END DO ! ji
380
381	!-----------------------------------------------------------------------------------------------
382	! Retrieve 2D arrays
383	!-----------------------------------------------------------------------------------------------
384
385	v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti)
386	v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti)
387	CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,jl) )
388	CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,jl) )
389	CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,jl) )
390	CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d (1:npti), v_ip (:,:,jl) )
391	CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d (1:npti), v_il (:,:,jl) )
392	DO jk = 1, nlay_i
393	CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl) )
394	END DO
395
396	END DO ! ji
397
398	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum )
399	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum )
400	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd )
401
402	!
403	ENDIF
404
405	END SUBROUTINE pnd_LEV
406
407	SUBROUTINE pnd_TOPO
408
409	!!-------------------------------------------------------------------
410	!! * ROUTINE pnd_TOPO *
411	!!
412	!! ** Purpose : Compute melt pond evolution
413	!!
414	!! ** Purpose : Compute melt pond evolution based on the ice
415	!! topography as inferred from the ice thickness
416	!! distribution.
417	!!
418	!! ** Method : This code is initially based on Flocco and Feltham
419	!! (2007) and Flocco et al. (2010). More to come...
420	!!
421	!! ** Tunable parameters :
422	!!
423	!! ** Note :
424	!!
425	!! ** References
426	!!
427	!! Flocco, D. and D. L. Feltham, 2007. A continuum model of melt pond
428	!! evolution on Arctic sea ice. J. Geophys. Res. 112, C08016, doi:
429	!! 10.1029/2006JC003836.
430	!!
431	!! Flocco, D., D. L. Feltham and A. K. Turner, 2010. Incorporation of
432	!! a physically based melt pond scheme into the sea ice component of a
433	!! climate model. J. Geophys. Res. 115, C08012,
434	!! doi: 10.1029/2009JC005568.
435	!!
436	!!-------------------------------------------------------------------
437
438	! local variables
439	REAL(wp) :: &
440	zdHui, & ! change in thickness of ice lid (m)
441	zomega, & ! conduction
442	zdTice, & ! temperature difference across ice lid (C)
443	zdvice, & ! change in ice volume (m)
444	zTavg, & ! mean surface temperature across categories (C)
445	zfsurf, & ! net heat flux, excluding conduction and transmitted radiation (W/m2)
446	zTp, & ! pond freezing temperature (C)
447	zrhoi_L, & ! volumetric latent heat of sea ice (J/m^3)
448	zfr_mlt, & ! fraction and volume of available meltwater retained for melt ponding
449	z1_rhow, & ! inverse water density
450	zpond , & ! dummy variable
451	zdum ! dummy variable
452
453	REAL(wp), DIMENSION(jpi,jpj) :: zvolp, & !! total melt pond water available before redistribution and drainage
454	zvolp_res
455
456	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i
457
458	INTEGER :: ji, jj, jk, jl ! loop indices
459
460	INTEGER :: i_test
461
462	! Note
463	! equivalent for CICE translation
464	! a_ip -> apond
465	! a_ip_frac -> apnd
466
467	!---------------------------------------------------------------
468	! Initialise
469	!---------------------------------------------------------------
470
471	! Parameters & constants (move to parameters)
472	zrhoi_L = rhoi * rLfus ! volumetric latent heat (J/m^3)
473	zTp = rt0 - 0.15_wp ! pond freezing point, slightly below 0C (ponds are bid saline)
474	z1_rhow = 1._wp / rhow
475
476	! Set required ice variables (hard-coded here for now)
477	! zfpond(:,:) = 0._wp ! contributing freshwater flux (?)
478
479	at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) ! ice fraction
480	vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) ! volume per grid area
481	vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) ! pond volume per grid area
482	vt_il(:,:) = SUM( v_il(:,:,:), dim=3 ) ! lid volume per grid area
483
484	! thickness
485	WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
486	ELSEWHERE ; z1_a_i(:,:,:) = 0._wp
487	END WHERE
488	h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
489
490	!---------------------------------------------------------------
491	! Change 2D to 1D
492	!---------------------------------------------------------------
493
494	!--------------------------------------------------------------
495	! Collect total available pond water
496	!--------------------------------------------------------------
497
498	zvolp(:,:) = 0._wp
499
500	DO jl = 1, jpl
501	DO jj = 1, jpj
502	DO ji = 1, jpi
503
504	IF ( a_i(ji,jj,jl) > epsi10 ) THEN
505
506	!--- Available meltwater for melt ponding ---!
507	zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * a_i(ji,jj,jl) ! = ( 1 - r ) = fraction of melt water that is not flushed
508	zdum = -( dh_i_sum_2d(ji,jj,jl)rhoi + dh_s_mlt_2d(ji,jj,jl)rhos ) * z1_rhow * a_i(ji,jj,jl)
509	zpond = zfr_mlt * zdum ! check
510
511	!--- Create possible new ponds
512	! if pond does not exist, create new pond over full ice area
513	IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN
514	h_ip(ji,jj,jl) = 0._wp
515	a_ip_frac(ji,jj,jl) = 1.0_wp ! pond fraction of sea ice (apnd for CICE)
516	ENDIF
517
518	!--- Deepen existing ponds before redistribution and drainage(later on)
519	h_ip(ji,jj,jl) = ( zpond + h_ip(ji,jj,jl) * a_ip_frac(ji,jj,jl) ) / a_ip_frac(ji,jj,jl)
520	zvolp(ji,jj) = zvolp(ji,jj) + h_ip(ji,jj,jl) * a_ip_frac(ji,jj,jl) * a_i(ji,jj,jl)
521	! zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) !useless
522
523	ENDIF ! a_i
524
525	END DO! jl
526	END DO ! jj
527	END DO ! ji
528
529	h_ip(:,:,:) = 0._wp ! pond thickness
530	a_ip(:,:,:) = 0._wp ! pond fraction per grid area
531
532	!--------------------------------------------------------------
533	! Redistribute and drain water from ponds
534	!--------------------------------------------------------------
535	CALL ice_thd_pnd_area( zvolp, zvolp_res )
536
537	!--------------------------------------------------------------
538	! Freeze and melt lid
539	!--------------------------------------------------------------
540	DO jj = 1, jpj
541	DO ji = 1, jpi
542
543	IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. vt_ip(ji,jj) > zvp_min *at_i(ji,jj) ) THEN
544
545	!--------------------------
546	! Pond lid growth and melt
547	!--------------------------
548	! Mean surface temperature
549	zTavg = 0._wp
550	DO jl = 1, jpl
551	zTavg = zTavg + t_su(ji,jj,jl)*a_i(ji,jj,jl)
552	END DO
553	zTavg = zTavg / a_i(ji,jj,jl) !!! could get a division by zero here
554
555	DO jl = 1, jpl-1
556
557	IF ( v_il(ji,jj,jl) > epsi10 ) THEN
558
559	!----------------------------------------------------------------
560	! Lid melting: floating upper ice layer melts in whole or part
561	!----------------------------------------------------------------
562	! Use Tsfc for each category
563	IF ( t_su(ji,jj,jl) > zTp ) THEN
564
565	zdvice = MIN( dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) )
566	IF ( zdvice > epsi10 ) then
567	v_il (ji,jj,jl) = v_il (ji,jj,jl) - zdvice
568	v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zdvice
569	! zvolp(ji,jj) = zvolp(ji,jj) + zdvice ! pointless to calculate total volume (done in icevar)
570	! zfpond(ji,jj) = fpond(ji,jj) + zdvice ! pointless to follow fw budget (ponds have no fw)
571
572	IF ( v_il(ji,jj,jl) < epsi10 .AND. v_ip(ji,jj,jl) > epsi10) THEN
573	! ice lid melted and category is pond covered
574	v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + v_il(ji,jj,jl)
575	! zfpond(ji,jj) = zfpond (ji,jj) + v_il(ji,jj,jl)
576	v_il(ji,jj,jl) = 0._wp
577	ENDIF
578	h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) !!! could get a division by zero here
579	ENDIF
580
581	!----------------------------------------------------------------
582	! Freeze pre-existing lid
583	!----------------------------------------------------------------
584
585	ELSE IF ( v_ip(ji,jj,jl) > epsi10 ) THEN ! Tsfcn(i,j,n) <= Tp
586
587	! differential growth of base of surface floating ice layer
588	zdTice = MAX( - t_su(ji,jj,jl) - zTd , 0._wp ) ! > 0
589	zomega = rcnd_i * zdTice / zrhoi_L
590	zdHui = SQRT( 2._wp * zomega * rdt_ice + ( v_il(ji,jj,jl) / a_i(ji,jj,jl) )**2 ) &
591	- v_il(ji,jj,jl) / a_i(ji,jj,jl)
592	zdvice = min( zdHui*a_ip(ji,jj,jl) , v_ip(ji,jj,jl) )
593
594	IF ( zdvice > epsi10 ) THEN
595	v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice
596	v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice
597	! zvolp(ji,jj) = zvolp(ji,jj) - zdvice
598	! zfpond(ji,jj) = zfpond(ji,jj) - zdvice
599	h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
600	ENDIF
601
602	ENDIF ! Tsfcn(i,j,n)
603
604	!----------------------------------------------------------------
605	! Freeze new lids
606	!----------------------------------------------------------------
607	! upper ice layer begins to form
608	! note: albedo does not change
609
610	ELSE ! v_il < epsi10
611
612	! thickness of newly formed ice
613	! the surface temperature of a meltpond is the same as that
614	! of the ice underneath (0C), and the thermodynamic surface
615	! flux is the same
616
617	!!! we need net surface energy flux, excluding conduction
618	!!! fsurf is summed over categories in CICE
619	!!! we have the category-dependent flux, let us use it ?
620	zfsurf = qns_ice(ji,jj,jl) + qsr_ice(ji,jj,jl)
621	zdHui = MAX ( -zfsurf * rdt_ice/zrhoi_L , 0._wp )
622	zdvice = MIN ( zdHui * a_ip(ji,jj,jl) , v_ip(ji,jj,jl) )
623	IF ( zdvice > epsi10 ) THEN
624	v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice
625	v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice
626	! zvolp(ji,jj) = zvolp(ji,jj) - zdvice
627	! zfpond(ji,jj) = zfpond(ji,jj) - zdvice
628	h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) ! in principle, this is useless as h_ip is computed in icevar
629	ENDIF
630
631	ENDIF ! v_il
632
633	END DO ! jl
634
635	ELSE ! remove ponds on thin ice
636
637	v_ip(ji,jj,:) = 0._wp
638	v_il(ji,jj,:) = 0._wp
639	! zfpond(ji,jj) = zfpond(ji,jj)- zvolp(ji,jj)
640	! zvolp(ji,jj) = 0._wp
641
642	ENDIF
643
644	END DO ! jj
645	END DO ! ji
646
647	!---------------------------------------------------------------
648	! Clean-up variables (probably duplicates what icevar would do)
649	!---------------------------------------------------------------
650	! MV comment
651	! In the ideal world, the lines above should update only v_ip, a_ip, v_il
652	! icevar should recompute all other variables (if needed at all)
653
654	DO jl = 1, jpl
655	DO jj = 1, jpj
656	DO ji = 1, jpi
657
658	IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 &
659	.AND. v_il (ji,jj,jl) > epsi10) THEN
660	v_il(ji,jj,jl) = 0._wp
661	ENDIF
662
663	! reload tracers
664	IF ( a_ip(ji,jj,jl) > epsi10) THEN
665	h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar
666	ELSE
667	v_il(ji,jj,jl) = 0._wp
668	h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as a_ip_frac computed in icevar
669	ENDIF
670
671	IF ( a_ip(ji,jj,jl) > epsi10 ) THEN
672	a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar
673	!h_ip(ji,jj,jl) = zhpondn(ji,jj,jl)
674	ELSE
675	a_ip_frac(ji,jj,jl) = 0._wp
676	h_ip(ji,jj,jl) = 0._wp ! MV in principle, useless as computed in icevar
677	h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as omputed in icevar
678	ENDIF
679
680	END DO ! ji
681	END DO ! jj
682	END DO ! jl
683
684	END SUBROUTINE pnd_TOPO
685
686
687	SUBROUTINE ice_thd_pnd_area( zvolp , zdvolp )
688
689	!!-------------------------------------------------------------------
690	!! * ROUTINE ice_thd_pnd_area *
691	!!
692	!! ** Purpose : Given the total volume of available pond water,
693	!! redistribute and drain water
694	!!
695	!! **
696	!!
697	!!------------------------------------------------------------------
698
699	REAL (wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: &
700	zvolp, & ! total available pond water
701	zdvolp ! remaining meltwater after redistribution
702
703	INTEGER :: &
704	ns, &
705	m_index, &
706	permflag
707
708	REAL (wp), DIMENSION(jpl) :: &
709	hicen, &
710	hsnon, &
711	asnon, &
712	alfan, &
713	betan, &
714	cum_max_vol, &
715	reduced_aicen
716
717	REAL (wp), DIMENSION(0:jpl) :: &
718	cum_max_vol_tmp
719
720	REAL (wp) :: &
721	hpond, &
722	drain, &
723	floe_weight, &
724	pressure_head, &
725	hsl_rel, &
726	deltah, &
727	perm, &
728	msno
729
730	REAL (wp), parameter :: &
731	viscosity = 1.79e-3_wp ! kinematic water viscosity in kg/m/s
732
733	INTEGER :: ji, jj, jk, jl ! loop indices
734
735	!-----------\|
736	! \|
737	! \|-----------\|
738	!___________\|___________\|______________________________________sea-level
739	! \| \|
740	! \| \|---^--------\|
741	! \| \| \| \|
742	! \| \| \| \|-----------\| \|-------
743	! \| \| \|alfan(jl)\| \| \|
744	! \| \| \| \| \|--------------\|
745	! \| \| \| \| \| \|
746	!---------------------------v-------------------------------------------
747	! \| \| ^ \| \| \|
748	! \| \| \| \| \|--------------\|
749	! \| \| \|betan(jl)\| \| \|
750	! \| \| \| \|-----------\| \|-------
751	! \| \| \| \|
752	! \| \|---v------- \|
753	! \| \|
754	! \|-----------\|
755	! \|
756	!-----------\|
757
758	a_ip(:,:,:) = 0._wp
759	h_ip(:,:,:) = 0._wp
760
761	DO jj = 1, jpj
762	DO ji = 1, jpi
763
764	IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. zvolp(ji,jj) > zvp_min * at_i(ji,jj) ) THEN
765
766	!-------------------------------------------------------------------
767	! initialize
768	!-------------------------------------------------------------------
769
770	DO jl = 1, jpl
771
772	a_ip(ji,jj,jl) = 0._wp
773	h_ip(ji,jj,jl) = 0._wp
774
775	!----------------------------------------
776	! compute the effective snow fraction
777	!----------------------------------------
778
779	IF (a_i(ji,jj,jl) < epsi10) THEN
780	hicen(jl) = 0._wp
781	hsnon(jl) = 0._wp
782	reduced_aicen(jl) = 0._wp
783	asnon(jl) = 0._wp !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables
784	ELSE
785	hicen(jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
786	hsnon(jl) = v_s(ji,jj,jl) / a_i(ji,jj,jl)
787	reduced_aicen(jl) = 1._wp ! n=jpl
788
789	!js: initial code in NEMO_DEV
790	!IF (n < jpl) reduced_aicen(jl) = aicen(jl) &
791	! * (-0.024_wp*hicen(jl) + 0.832_wp)
792
793	!js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large
794	IF (jl < jpl) reduced_aicen(jl) = a_i(ji,jj,jl) &
795	* max(0.2_wp,(-0.024_wp*hicen(jl) + 0.832_wp))
796
797	asnon(jl) = reduced_aicen(jl) ! effective snow fraction (empirical)
798	! MV should check whether this makes sense to have the same effective snow fraction in here
799	! OLI: it probably doesn't
800	END IF
801
802	! This choice for alfa and beta ignores hydrostatic equilibium of categories.
803	! Hydrostatic equilibium of the entire ITD is accounted for below, assuming
804	! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all
805	! categories. alfa and beta partition the ITD - they are areas not thicknesses!
806	! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area.
807	! Here, alfa = 60% of the ice area (and since hice is constant in a category,
808	! alfan = 60% of the ice volume) in each category lies above the reference line,
809	! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required.
810
811	! MV:
812	! Note that this choice is not in the original FF07 paper and has been adopted in CICE
813	! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe
814
815	! Where does that choice come from ? => OLI : Coz' Chuck Norris said so...
816
817	alfan(jl) = 0.6 * hicen(jl)
818	betan(jl) = 0.4 * hicen(jl)
819
820	cum_max_vol(jl) = 0._wp
821	cum_max_vol_tmp(jl) = 0._wp
822
823	END DO ! jpl
824
825	cum_max_vol_tmp(0) = 0._wp
826	drain = 0._wp
827	zdvolp(ji,jj) = 0._wp
828
829	!----------------------------------------------------------
830	! Drain overflow water, update pond fraction and volume
831	!----------------------------------------------------------
832
833	!--------------------------------------------------------------------------
834	! the maximum amount of water that can be contained up to each ice category
835	!--------------------------------------------------------------------------
836	! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW)
837	! Then the excess volume cum_max_vol(jl) drains out of the system
838	! It should be added to wfx_pnd_out
839
840	DO jl = 1, jpl-1 ! last category can not hold any volume
841
842	IF (alfan(jl+1) >= alfan(jl) .AND. alfan(jl+1) > 0._wp ) THEN
843
844	! total volume in level including snow
845	cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) + &
846	(alfan(jl+1) - alfan(jl)) * sum(reduced_aicen(1:jl))
847
848	! subtract snow solid volumes from lower categories in current level
849	DO ns = 1, jl
850	cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) &
851	- rhos/rhow * & ! free air fraction that can be filled by water
852	asnon(ns) * & ! effective areal fraction of snow in that category
853	max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp)
854	END DO
855
856	ELSE ! assume higher categories unoccupied
857	cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1)
858	END IF
859	!IF (cum_max_vol_tmp(jl) < z0) THEN
860	! CALL abort_ice('negative melt pond volume')
861	!END IF
862	END DO
863	cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1) ! last category holds no volume
864	cum_max_vol (1:jpl) = cum_max_vol_tmp(1:jpl)
865
866	!----------------------------------------------------------------
867	! is there more meltwater than can be held in the floe?
868	!----------------------------------------------------------------
869	IF (zvolp(ji,jj) >= cum_max_vol(jpl)) THEN
870	drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10
871	zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available
872	zdvolp(ji,jj) = drain ! this is the drained water
873	IF (zvolp(ji,jj) < epsi10) THEN
874	zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj)
875	zvolp(ji,jj) = 0._wp
876	END IF
877	END IF
878
879	! height and area corresponding to the remaining volume
880	CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index)
881
882	DO jl = 1, m_index
883	!h_ip(jl) = hpond - alfan(jl) + alfan(1) ! here oui choulde update
884	! ! volume instead, no ?
885	h_ip(ji,jj,jl) = max((hpond - alfan(jl) + alfan(1)), 0._wp) !js: from CICE 5.1.2
886	a_ip(ji,jj,jl) = reduced_aicen(jl)
887	! in practise, pond fraction depends on the empirical snow fraction
888	! so in turn on ice thickness
889	END DO
890	!zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d
891
892	!------------------------------------------------------------------------
893	! Drainage through brine network (permeability)
894	!------------------------------------------------------------------------
895	!!! drainage due to ice permeability - Darcy's law
896
897	! sea water level
898	msno = 0._wp
899	DO jl = 1 , jpl
900	msno = msno + v_s(ji,jj,jl) * rhos
901	END DO
902	floe_weight = ( msno + rhoivt_i(ji,jj) + rau0zvolp(ji,jj) ) / at_i(ji,jj)
903	hsl_rel = floe_weight / rau0 &
904	- ( ( sum(betan(:)*a_i(ji,jj,:)) / at_i(ji,jj) ) + alfan(1) )
905
906	deltah = hpond - hsl_rel
907	pressure_head = grav * rau0 * max(deltah, 0._wp)
908
909	! drain if ice is permeable
910	permflag = 0
911
912	IF (pressure_head > 0._wp) THEN
913	DO jl = 1, jpl-1
914	IF ( hicen(jl) /= 0._wp ) THEN
915
916	!IF (hicen(jl) > 0._wp) THEN !js: from CICE 5.1.2
917
918	perm = 0._wp ! MV ugly dummy patch
919	! CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl), sz_i(ji,jj,:,jl), perm) ! bof
920	IF (perm > 0._wp) permflag = 1
921
922	drain = perma_ip(ji,jj,jl)pressure_head*rdt_ice / &
923	(viscosity*hicen(jl))
924	zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj))
925	zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp)
926	IF (zvolp(ji,jj) < epsi10) THEN
927	zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj)
928	zvolp(ji,jj) = 0._wp
929	END IF
930	END IF
931	END DO
932
933	! adjust melt pond dimensions
934	IF (permflag > 0) THEN
935	! recompute pond depth
936	CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index)
937	DO jl = 1, m_index
938	h_ip(ji,jj,jl) = hpond - alfan(jl) + alfan(1)
939	a_ip(ji,jj,jl) = reduced_aicen(jl)
940	END DO
941	!zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d
942	END IF
943	END IF ! pressure_head
944
945	!-------------------------------
946	! remove water from the snow
947	!-------------------------------
948	!------------------------------------------------------------------------
949	! total melt pond volume in category does not include snow volume
950	! snow in melt ponds is not melted
951	!------------------------------------------------------------------------
952
953	! Calculate pond volume for lower categories
954	DO jl = 1,m_index-1
955	v_ip(ji,jj,jl) = a_ip(ji,jj,jl) * h_ip(ji,jj,jl) & ! what is not in the snow
956	- (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl))
957	END DO
958
959	! Calculate pond volume for highest category = remaining pond volume
960
961	! The following is completely unclear to Martin at least
962	! Could we redefine properly and recode in a more readable way ?
963
964	! m_index = last category with melt pond
965
966	IF (m_index == 1) v_ip(ji,jj,m_index) = zvolp(ji,jj) ! volume of mw in 1st category is the total volume of melt water
967
968	IF (m_index > 1) THEN
969	IF (zvolp(ji,jj) > sum( v_ip(ji,jj,1:m_index-1))) THEN
970	v_ip(ji,jj,m_index) = zvolp(ji,jj) - sum(v_ip(ji,jj,1:m_index-1))
971	ELSE
972	v_ip(ji,jj,m_index) = 0._wp
973	h_ip(ji,jj,m_index) = 0._wp
974	a_ip(ji,jj,m_index) = 0._wp
975	! If remaining pond volume is negative reduce pond volume of
976	! lower category
977	IF ( zvolp(ji,jj) + epsi10 < SUM(v_ip(ji,jj,1:m_index-1))) &
978	v_ip(ji,jj,m_index-1) = v_ip(ji,jj,m_index-1) - sum(v_ip(ji,jj,1:m_index-1)) + zvolp(ji,jj)
979	END IF
980	END IF
981
982	DO jl = 1,m_index
983	IF (a_ip(ji,jj,jl) > epsi10) THEN
984	h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
985	ELSE
986	zdvolp(ji,jj) = zdvolp(ji,jj) + v_ip(ji,jj,jl)
987	h_ip(ji,jj,jl) = 0._wp
988	v_ip(ji,jj,jl) = 0._wp
989	a_ip(ji,jj,jl) = 0._wp
990	END IF
991	END DO
992	DO jl = m_index+1, jpl
993	h_ip(ji,jj,jl) = 0._wp
994	a_ip(ji,jj,jl) = 0._wp
995	v_ip(ji,jj,jl) = 0._wp
996	END DO
997
998	ENDIF
999	END DO ! ji
1000	END DO ! jj
1001
1002	END SUBROUTINE ice_thd_pnd_area
1003
1004
1005	SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index)
1006	!!-------------------------------------------------------------------
1007	!! * ROUTINE ice_thd_pnd_depth *
1008	!!
1009	!! ** Purpose : Compute melt pond depth
1010	!!-------------------------------------------------------------------
1011
1012	REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
1013	aicen, &
1014	asnon, &
1015	hsnon, &
1016	alfan, &
1017	cum_max_vol
1018
1019	REAL (wp), INTENT(IN) :: &
1020	zvolp
1021
1022	REAL (wp), INTENT(OUT) :: &
1023	hpond
1024
1025	INTEGER, INTENT(OUT) :: &
1026	m_index
1027
1028	INTEGER :: n, ns
1029
1030	REAL (wp), DIMENSION(0:jpl+1) :: &
1031	hitl, &
1032	aicetl
1033
1034	REAL (wp) :: &
1035	rem_vol, &
1036	area, &
1037	vol, &
1038	tmp, &
1039	z0 = 0.0_wp
1040
1041	!----------------------------------------------------------------
1042	! hpond is zero if zvolp is zero - have we fully drained?
1043	!----------------------------------------------------------------
1044
1045	IF (zvolp < epsi10) THEN
1046	hpond = z0
1047	m_index = 0
1048	ELSE
1049
1050	!----------------------------------------------------------------
1051	! Calculate the category where water fills up to
1052	!----------------------------------------------------------------
1053
1054	!----------\|
1055	! \|
1056	! \|
1057	! \|----------\| -- --
1058	!__________\|__________\|_________________________________________ ^
1059	! \| \| rem_vol ^ \| Semi-filled
1060	! \| \|----------\|-- -- -- - ---\|-- ---- -- -- --v layer
1061	! \| \| \| \|
1062	! \| \| \| \|hpond
1063	! \| \| \|----------\| \| \|-------
1064	! \| \| \| \| \| \|
1065	! \| \| \| \|---v-----\|
1066	! \| \| m_index \| \| \|
1067	!-------------------------------------------------------------
1068
1069	m_index = 0 ! 1:m_index categories have water in them
1070	DO n = 1, jpl
1071	IF (zvolp <= cum_max_vol(n)) THEN
1072	m_index = n
1073	IF (n == 1) THEN
1074	rem_vol = zvolp
1075	ELSE
1076	rem_vol = zvolp - cum_max_vol(n-1)
1077	END IF
1078	exit ! to break out of the loop
1079	END IF
1080	END DO
1081	m_index = min(jpl-1, m_index)
1082
1083	!----------------------------------------------------------------
1084	! semi-filled layer may have m_index different snow in it
1085	!----------------------------------------------------------------
1086
1087	!----------------------------------------------------------- ^
1088	! \| alfan(m_index+1)
1089	! \|
1090	!hitl(3)--> \|----------\| \|
1091	!hitl(2)--> \|------------\| * * * * *\| \|
1092	!hitl(1)--> \|----------\|* * * * * * \|* * * * * \| \|
1093	!hitl(0)-->------------------------------------------------- \| ^
1094	! various snow from lower categories \| \|alfa(m_index)
1095
1096	! hitl - heights of the snow layers from thinner and current categories
1097	! aicetl - area of each snow depth in this layer
1098
1099	hitl(:) = z0
1100	aicetl(:) = z0
1101	DO n = 1, m_index
1102	hitl(n) = max(min(hsnon(n) + alfan(n) - alfan(m_index), &
1103	alfan(m_index+1) - alfan(m_index)), z0)
1104	aicetl(n) = asnon(n)
1105
1106	aicetl(0) = aicetl(0) + (aicen(n) - asnon(n))
1107	END DO
1108
1109	hitl(m_index+1) = alfan(m_index+1) - alfan(m_index)
1110	aicetl(m_index+1) = z0
1111
1112	!----------------------------------------------------------------
1113	! reorder array according to hitl
1114	! snow heights not necessarily in height order
1115	!----------------------------------------------------------------
1116
1117	DO ns = 1, m_index+1
1118	DO n = 0, m_index - ns + 1
1119	IF (hitl(n) > hitl(n+1)) THEN ! swap order
1120	tmp = hitl(n)
1121	hitl(n) = hitl(n+1)
1122	hitl(n+1) = tmp
1123	tmp = aicetl(n)
1124	aicetl(n) = aicetl(n+1)
1125	aicetl(n+1) = tmp
1126	END IF
1127	END DO
1128	END DO
1129
1130	!----------------------------------------------------------------
1131	! divide semi-filled layer into set of sublayers each vertically homogenous
1132	!----------------------------------------------------------------
1133
1134	!hitl(3)----------------------------------------------------------------
1135	! \| * * * * * * * *
1136	! \|* * * * * * * * *
1137	!hitl(2)----------------------------------------------------------------
1138	! \| * * * * * * * * \| * * * * * * * *
1139	! \|* * * * * * * * * \|* * * * * * * * *
1140	!hitl(1)----------------------------------------------------------------
1141	! \| * * * * * * * * \| * * * * * * * * \| * * * * * * * *
1142	! \|* * * * * * * * * \|* * * * * * * * * \|* * * * * * * * *
1143	!hitl(0)----------------------------------------------------------------
1144	! aicetl(0) aicetl(1) aicetl(2) aicetl(3)
1145
1146	! move up over layers incrementing volume
1147	DO n = 1, m_index+1
1148
1149	area = sum(aicetl(:)) - & ! total area of sub-layer
1150	(rhos/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow
1151
1152	vol = (hitl(n) - hitl(n-1)) * area ! thickness of sub-layer times area
1153
1154	IF (vol >= rem_vol) THEN ! have reached the sub-layer with the depth within
1155	hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1)
1156
1157	exit
1158	ELSE ! still in sub-layer below the sub-layer with the depth
1159	rem_vol = rem_vol - vol
1160	END IF
1161
1162	END DO
1163
1164	END IF
1165
1166	END SUBROUTINE ice_thd_pnd_depth
1167
1168
1169	SUBROUTINE ice_thd_pnd_perm(ticen, salin, perm)
1170	!!-------------------------------------------------------------------
1171	!! * ROUTINE ice_thd_pnd_perm *
1172	!!
1173	!! ** Purpose : Determine the liquid fraction of brine in the ice
1174	!! and its permeability
1175	!!-------------------------------------------------------------------
1176
1177	REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: &
1178	ticen, & ! internal ice temperature (K)
1179	salin ! salinity (ppt) !js: ppt according to cice
1180
1181	REAL (wp), INTENT(OUT) :: &
1182	perm ! permeability
1183
1184	REAL (wp) :: &
1185	Sbr ! brine salinity
1186
1187	REAL (wp), DIMENSION(nlay_i) :: &
1188	Tin, & ! ice temperature
1189	phi ! liquid fraction
1190
1191	INTEGER :: k
1192
1193	!-----------------------------------------------------------------
1194	! Compute ice temperatures from enthalpies using quadratic formula
1195	!-----------------------------------------------------------------
1196
1197	DO k = 1,nlay_i
1198	Tin(k) = ticen(k) - rt0 !js: from K to degC
1199	END DO
1200
1201	!-----------------------------------------------------------------
1202	! brine salinity and liquid fraction
1203	!-----------------------------------------------------------------
1204
1205	DO k = 1, nlay_i
1206
1207	Sbr = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus)
1208	! Best expression to date is that one
1209	Sbr = - 18.7 * Tin(k) - 0.519 * Tin(k)*2 - 0.00535 Tin(k) **3
1210	phi(k) = salin(k) / Sbr
1211
1212	END DO
1213
1214	!-----------------------------------------------------------------
1215	! permeability
1216	!-----------------------------------------------------------------
1217
1218	perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007)
1219
1220	END SUBROUTINE ice_thd_pnd_perm
1221
1222
1223	!----------------------------------------------------------------------------------------------------------------------
1224
1225	SUBROUTINE ice_thd_pnd_init
1226	!!-------------------------------------------------------------------
1227	!! * ROUTINE ice_thd_pnd_init *
1228	!!
1229	!! ** Purpose : Physical constants and parameters linked to melt ponds
1230	!! over sea ice
1231	!!
1232	!! ** Method : Read the namthd_pnd namelist and check the melt pond
1233	!! parameter values called at the first timestep (nit000)
1234	!!
1235	!! ** input : Namelist namthd_pnd
1236	!!-------------------------------------------------------------------
1237	INTEGER :: ios, ioptio ! Local integer
1238	!!
1239	NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, &
1240	& ln_pnd_CST , rn_apnd, rn_hpnd, &
1241	& ln_pnd_TOPO , &
1242	& ln_pnd_lids, ln_pnd_alb
1243	!!-------------------------------------------------------------------
1244	!
1245	REWIND( numnam_ice_ref ) ! Namelist namthd_pnd in reference namelist : Melt Ponds
1246	READ ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901)
1247	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_pnd in reference namelist' )
1248	REWIND( numnam_ice_cfg ) ! Namelist namthd_pnd in configuration namelist : Melt Ponds
1249	READ ( numnam_ice_cfg, namthd_pnd, IOSTAT = ios, ERR = 902 )
1250	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namthd_pnd in configuration namelist' )
1251	IF(lwm) WRITE ( numoni, namthd_pnd )
1252	!
1253	IF(lwp) THEN ! control print
1254	WRITE(numout,*)
1255	WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
1256	WRITE(numout,*) '~~~~~~~~~~~~~~~~'
1257	WRITE(numout,*) ' Namelist namicethd_pnd:'
1258	WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd
1259	WRITE(numout,*) ' Topographic melt pond scheme ln_pnd_TOPO = ', ln_pnd_TOPO
1260	WRITE(numout,*) ' Level ice melt pond scheme ln_pnd_LEV = ', ln_pnd_LEV
1261	WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_apnd_min = ', rn_apnd_min
1262	WRITE(numout,*) ' Maximum ice fraction that contributes to melt ponds rn_apnd_max = ', rn_apnd_max
1263	WRITE(numout,*) ' Constant ice melt pond scheme ln_pnd_CST = ', ln_pnd_CST
1264	WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd
1265	WRITE(numout,*) ' Prescribed pond depth rn_hpnd = ', rn_hpnd
1266	WRITE(numout,*) ' Frozen lids on top of melt ponds ln_pnd_lids = ', ln_pnd_lids
1267	WRITE(numout,*) ' Melt ponds affect albedo or not ln_pnd_alb = ', ln_pnd_alb
1268	ENDIF
1269	!
1270	! !== set the choice of ice pond scheme ==!
1271	ioptio = 0
1272	IF( .NOT.ln_pnd ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndNO ; ENDIF
1273	IF( ln_pnd_CST ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndCST ; ENDIF
1274	IF( ln_pnd_LEV ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndLEV ; ENDIF
1275	IF( ln_pnd_TOPO ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndTOPO ; ENDIF
1276	IF( ioptio /= 1 ) &
1277	& CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_LEV or ln_pnd_CST)' )
1278	!
1279	SELECT CASE( nice_pnd )
1280	CASE( np_pndNO )
1281	IF( ln_pnd_alb ) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF
1282	IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when no ponds' ) ; ENDIF
1283	CASE( np_pndCST )
1284	IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when constant ponds' ) ; ENDIF
1285	END SELECT
1286	!
1287	END SUBROUTINE ice_thd_pnd_init
1288
1289	#else
1290	!!----------------------------------------------------------------------
1291	!! Default option Empty module NO SI3 sea-ice model
1292	!!----------------------------------------------------------------------
1293	#endif
1294
1295	!!======================================================================
1296	END MODULE icethd_pnd

Note: See TracBrowser for help on using the repository browser.

Download in other formats: