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 @ 13891

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

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

Download in other formats: