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/UKMO/NEMO_4.0.4_GO8_package/src/ICE – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/ICE/icethd_pnd.F90 @ 14335

Last change on this file since 14335 was 14302, checked in by edblockley, 3 years ago
Porting topographic melt-ponds from Martin's branch branches/2020/SI3-05_MeltPonds_topo@14301
File size: 67.2 KB

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

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

Download in other formats: