Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

icethd_pnd.F90

Last change on this file was 15123, checked in by edblockley, 3 years ago
Bug fix to TOPO melt-pond scheme to fix issue of non-melting lids (identified by David Livings@CPOM)
File size: 67.6 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	zvolp_res(:,:) = 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 and contributing meltwater for melt ponding ---!
593	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
594	& * z1_rhow * a_i(ji,jj,jl)
595	! MV -> could move this directly in ice_thd_dh and get an array (ji,jj,jl) for surface melt water volume per grid area
596	zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i(ji,jj) ! fraction of surface meltwater going to ponds
597	zv_pnd = zfr_mlt * zv_mlt ! contributing meltwater volume for category jl
598
599	diag_dvpn_mlt(ji,jj) = diag_dvpn_mlt(ji,jj) + zv_mlt * r1_rdtice ! diags
600	diag_dvpn_rnf(ji,jj) = diag_dvpn_rnf(ji,jj) + ( 1. - zfr_mlt ) * zv_mlt * r1_rdtice
601
602	!--- Create possible new ponds
603	! if pond does not exist, create new pond over full ice area
604	!IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN
605	IF ( a_ip(ji,jj,jl) < epsi10 ) THEN
606	a_ip(ji,jj,jl) = a_i(ji,jj,jl)
607	a_ip_frac(ji,jj,jl) = 1.0_wp ! pond fraction of sea ice (apnd for CICE)
608	ENDIF
609
610	!--- Deepen existing ponds with no change in pond fraction, before redistribution and drainage
611	v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zv_pnd ! use pond water to increase thickness
612	h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
613
614	!--- Total available pond water volume (pre-existing + newly produced)j
615	zvolp(ji,jj) = zvolp(ji,jj) + v_ip(ji,jj,jl)
616	! zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) ! useless for now
617
618	ENDIF ! a_i
619
620	END DO! jl
621	END DO ! jj
622	END DO ! ji
623
624	zvolp_ini(:,:) = zvolp(:,:)
625
626	!--------------------------------------------------------------
627	! Redistribute and drain water from ponds
628	!--------------------------------------------------------------
629	CALL ice_thd_pnd_area( zvolp, zvolp_res )
630
631	!--------------------------------------------------------------
632	! Melt pond lid growth and melt
633	!--------------------------------------------------------------
634
635	IF( ln_pnd_lids ) THEN
636
637	DO jj = 1, jpj
638	DO ji = 1, jpi
639
640	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
641
642	!--------------------------
643	! Pond lid growth and melt
644	!--------------------------
645	! Mean surface temperature
646	zTavg = 0._wp
647	DO jl = 1, jpl
648	zTavg = zTavg + t_su(ji,jj,jl)*a_i(ji,jj,jl)
649	END DO
650	zTavg = zTavg / at_i(ji,jj) !!! could get a division by zero here
651
652	DO jl = 1, jpl-1
653
654	IF ( v_il(ji,jj,jl) > epsi10 ) THEN
655
656	!----------------------------------------------------------------
657	! Lid melting: floating upper ice layer melts in whole or part
658	!----------------------------------------------------------------
659	! Use Tsfc for each category
660	IF ( t_su(ji,jj,jl) > zTp ) THEN
661
662	zdvice = MIN( -dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) )
663
664	IF ( zdvice > epsi10 ) THEN
665
666	v_il (ji,jj,jl) = v_il (ji,jj,jl) - zdvice
667	v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zdvice ! MV: not sure i understand dh_i_sum seems counted twice -
668	! as it is already counted in surface melt
669	! zvolp(ji,jj) = zvolp(ji,jj) + zdvice ! pointless to calculate total volume (done in icevar)
670	! zfpond(ji,jj) = fpond(ji,jj) + zdvice ! pointless to follow fw budget (ponds have no fw)
671
672	IF ( v_il(ji,jj,jl) < epsi10 .AND. v_ip(ji,jj,jl) > epsi10) THEN
673	! ice lid melted and category is pond covered
674	v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + v_il(ji,jj,jl)
675	! zfpond(ji,jj) = zfpond (ji,jj) + v_il(ji,jj,jl)
676	v_il(ji,jj,jl) = 0._wp
677	ENDIF
678	h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) !!! could get a division by zero here
679
680	diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice ! diag
681
682	ENDIF
683
684	!----------------------------------------------------------------
685	! Freeze pre-existing lid
686	!----------------------------------------------------------------
687
688	ELSE IF ( v_ip(ji,jj,jl) > epsi10 ) THEN ! Tsfcn(i,j,n) <= Tp
689
690	! differential growth of base of surface floating ice layer
691	! zdTice = MAX( - t_su(ji,jj,jl) - zTd , 0._wp ) ! > 0 ! line valid for temp in C
692	zdTice = MAX( - ( t_su(ji,jj,jl) - zTd ) , 0._wp ) ! > 0
693	zomega = rcnd_i * zdTice / zrhoi_L
694	zdHui = SQRT( 2._wp * zomega * rdt_ice + ( v_il(ji,jj,jl) / a_i(ji,jj,jl) )**2 ) &
695	- v_il(ji,jj,jl) / a_i(ji,jj,jl)
696	zdvice = min( zdHui*a_ip(ji,jj,jl) , v_ip(ji,jj,jl) )
697
698	IF ( zdvice > epsi10 ) THEN
699	v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice
700	v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice
701	! zvolp(ji,jj) = zvolp(ji,jj) - zdvice
702	! zfpond(ji,jj) = zfpond(ji,jj) - zdvice
703	h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
704
705	diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag
706
707	ENDIF
708
709	ENDIF ! Tsfcn(i,j,n)
710
711	!----------------------------------------------------------------
712	! Freeze new lids
713	!----------------------------------------------------------------
714	! upper ice layer begins to form
715	! note: albedo does not change
716
717	ELSE ! v_il < epsi10
718
719	! thickness of newly formed ice
720	! the surface temperature of a meltpond is the same as that
721	! of the ice underneath (0C), and the thermodynamic surface
722	! flux is the same
723
724	!!! we need net surface energy flux, excluding conduction
725	!!! fsurf is summed over categories in CICE
726	!!! we have the category-dependent flux, let us use it ?
727	zfsurf = qns_ice(ji,jj,jl) + qsr_ice(ji,jj,jl)
728	zdHui = MAX ( -zfsurf * rdt_ice / zrhoi_L , 0._wp )
729	zdvice = MIN ( zdHui * a_ip(ji,jj,jl) , v_ip(ji,jj,jl) )
730
731	IF ( zdvice > epsi10 ) THEN
732
733	v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice
734	v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice
735	diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag
736
737	! zvolp(ji,jj) = zvolp(ji,jj) - zdvice
738	! zfpond(ji,jj) = zfpond(ji,jj) - zdvice
739
740	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
741	ENDIF
742
743	ENDIF ! v_il
744
745	!----------------------------------------------------------------
746	! Remove ice lid if there is no liquid pond
747	!----------------------------------------------------------------
748
749	IF (v_ip(ji,jj,jl) < epsi10) THEN
750	v_il(ji,jj,jl) = 0._wp
751	ENDIF
752
753	END DO ! jl
754
755	ELSE ! remove ponds on thin ice
756
757	v_ip(ji,jj,:) = 0._wp
758	v_il(ji,jj,:) = 0._wp
759	! zfpond(ji,jj) = zfpond(ji,jj)- zvolp(ji,jj)
760	! zvolp(ji,jj) = 0._wp
761
762	ENDIF
763
764	END DO ! jj
765	END DO ! ji
766
767	ENDIF ! ln_pnd_lids
768
769	!---------------------------------------------------------------
770	! Clean-up variables (probably duplicates what icevar would do)
771	!---------------------------------------------------------------
772	! MV comment
773	! In the ideal world, the lines above should update only v_ip, a_ip, v_il
774	! icevar should recompute all other variables (if needed at all)
775
776	DO jl = 1, jpl
777
778	DO jj = 1, jpj
779	DO ji = 1, jpi
780
781	! ! zap lids on small ponds
782	! IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 &
783	! .AND. v_il(ji,jj,jl) > epsi10) THEN
784	! v_il(ji,jj,jl) = 0._wp ! probably uselesss now since we get zap_small
785	! ENDIF
786
787	! recalculate equivalent pond variables
788	IF ( a_ip(ji,jj,jl) > epsi10) THEN
789	h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
790	a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar
791	h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar
792	ENDIF
793	! h_ip(ji,jj,jl) = 0._wp ! MV in principle, useless as computed in icevar
794	! h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as omputed in icevar
795	! ENDIF
796
797	END DO ! ji
798	END DO ! jj
799
800	END DO ! jl
801
802	END SUBROUTINE pnd_TOPO
803
804
805	SUBROUTINE ice_thd_pnd_area( zvolp , zdvolp )
806
807	!!-------------------------------------------------------------------
808	!! * ROUTINE ice_thd_pnd_area *
809	!!
810	!! ** Purpose : Given the total volume of available pond water,
811	!! redistribute and drain water
812	!!
813	!! ** Method
814	!!
815	!-----------\|
816	! \|
817	! \|-----------\|
818	!___________\|___________\|______________________________________sea-level
819	! \| \|
820	! \| \|---^--------\|
821	! \| \| \| \|
822	! \| \| \| \|-----------\| \|-------
823	! \| \| \| alfan \| \| \|
824	! \| \| \| \| \|--------------\|
825	! \| \| \| \| \| \|
826	!---------------------------v-------------------------------------------
827	! \| \| ^ \| \| \|
828	! \| \| \| \| \|--------------\|
829	! \| \| \| betan \| \| \|
830	! \| \| \| \|-----------\| \|-------
831	! \| \| \| \|
832	! \| \|---v------- \|
833	! \| \|
834	! \|-----------\|
835	! \|
836	!-----------\|
837	!
838	!!
839	!!------------------------------------------------------------------
840
841	REAL (wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: &
842	zvolp, & ! total available pond water
843	zdvolp ! remaining meltwater after redistribution
844
845	INTEGER :: &
846	ns, &
847	m_index, &
848	permflag
849
850	REAL (wp), DIMENSION(jpl) :: &
851	hicen, &
852	hsnon, &
853	asnon, &
854	alfan, &
855	betan, &
856	cum_max_vol, &
857	reduced_aicen
858
859	REAL (wp), DIMENSION(0:jpl) :: &
860	cum_max_vol_tmp
861
862	REAL (wp) :: &
863	hpond, &
864	drain, &
865	floe_weight, &
866	pressure_head, &
867	hsl_rel, &
868	deltah, &
869	perm, &
870	msno
871
872	REAL (wp), parameter :: &
873	viscosity = 1.79e-3_wp ! kinematic water viscosity in kg/m/s
874
875	INTEGER :: ji, jj, jk, jl ! loop indices
876
877	a_ip(:,:,:) = 0._wp
878	h_ip(:,:,:) = 0._wp
879
880	DO jj = 1, jpj
881	DO ji = 1, jpi
882
883	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
884
885	!-------------------------------------------------------------------
886	! initialize
887	!-------------------------------------------------------------------
888
889	DO jl = 1, jpl
890
891	!----------------------------------------
892	! compute the effective snow fraction
893	!----------------------------------------
894
895	IF (a_i(ji,jj,jl) < epsi10) THEN
896	hicen(jl) = 0._wp
897	hsnon(jl) = 0._wp
898	reduced_aicen(jl) = 0._wp
899	asnon(jl) = 0._wp !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables
900	ELSE
901	hicen(jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
902	hsnon(jl) = v_s(ji,jj,jl) / a_i(ji,jj,jl)
903	reduced_aicen(jl) = 1._wp ! n=jpl
904
905	!js: initial code in NEMO_DEV
906	!IF (n < jpl) reduced_aicen(jl) = aicen(jl) &
907	! * (-0.024_wp*hicen(jl) + 0.832_wp)
908
909	!js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large
910	IF (jl < jpl) reduced_aicen(jl) = a_i(ji,jj,jl) &
911	* max(0.2_wp,(-0.024_wp*hicen(jl) + 0.832_wp))
912
913	asnon(jl) = reduced_aicen(jl) ! effective snow fraction (empirical)
914	! MV should check whether this makes sense to have the same effective snow fraction in here
915	! OLI: it probably doesn't
916	END IF
917
918	! This choice for alfa and beta ignores hydrostatic equilibium of categories.
919	! Hydrostatic equilibium of the entire ITD is accounted for below, assuming
920	! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all
921	! categories. alfa and beta partition the ITD - they are areas not thicknesses!
922	! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area.
923	! Here, alfa = 60% of the ice area (and since hice is constant in a category,
924	! alfan = 60% of the ice volume) in each category lies above the reference line,
925	! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required.
926
927	! MV:
928	! Note that this choice is not in the original FF07 paper and has been adopted in CICE
929	! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe
930
931	! Where does that choice come from ? => OLI : Coz' Chuck Norris said so...
932
933	alfan(jl) = 0.6 * hicen(jl)
934	betan(jl) = 0.4 * hicen(jl)
935
936	cum_max_vol(jl) = 0._wp
937	cum_max_vol_tmp(jl) = 0._wp
938
939	END DO ! jpl
940
941	cum_max_vol_tmp(0) = 0._wp
942	drain = 0._wp
943	zdvolp(ji,jj) = 0._wp
944
945	!----------------------------------------------------------
946	! Drain overflow water, update pond fraction and volume
947	!----------------------------------------------------------
948
949	!--------------------------------------------------------------------------
950	! the maximum amount of water that can be contained up to each ice category
951	!--------------------------------------------------------------------------
952	! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW)
953	! Then the excess volume cum_max_vol(jl) drains out of the system
954	! It should be added to wfx_pnd_out
955
956	DO jl = 1, jpl-1 ! last category can not hold any volume
957
958	IF (alfan(jl+1) >= alfan(jl) .AND. alfan(jl+1) > 0._wp ) THEN
959
960	! total volume in level including snow
961	cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) + &
962	(alfan(jl+1) - alfan(jl)) * sum(reduced_aicen(1:jl))
963
964	! subtract snow solid volumes from lower categories in current level
965	DO ns = 1, jl
966	cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) &
967	- rhos/rhow * & ! free air fraction that can be filled by water
968	asnon(ns) * & ! effective areal fraction of snow in that category
969	max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp)
970	END DO
971
972	ELSE ! assume higher categories unoccupied
973	cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1)
974	END IF
975	!IF (cum_max_vol_tmp(jl) < z0) THEN
976	! CALL abort_ice('negative melt pond volume')
977	!END IF
978	END DO
979	cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1) ! last category holds no volume
980	cum_max_vol (1:jpl) = cum_max_vol_tmp(1:jpl)
981
982	!----------------------------------------------------------------
983	! is there more meltwater than can be held in the floe?
984	!----------------------------------------------------------------
985	IF (zvolp(ji,jj) >= cum_max_vol(jpl)) THEN
986	drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10
987	zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available
988
989	diag_dvpn_rnf(ji,jj) = - drain ! diag - overflow counted in the runoff part (arbitrary choice)
990
991	zdvolp(ji,jj) = drain ! this is the drained water
992	IF (zvolp(ji,jj) < epsi10) THEN
993	zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj)
994	zvolp(ji,jj) = 0._wp
995	END IF
996	END IF
997
998	! height and area corresponding to the remaining volume
999	! routine leaves zvolp unchanged
1000	CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index)
1001
1002	DO jl = 1, m_index
1003	!h_ip(jl) = hpond - alfan(jl) + alfan(1) ! here oui choulde update
1004	! ! volume instead, no ?
1005	h_ip(ji,jj,jl) = max((hpond - alfan(jl) + alfan(1)), 0._wp) !js: from CICE 5.1.2
1006	a_ip(ji,jj,jl) = reduced_aicen(jl)
1007	! in practise, pond fraction depends on the empirical snow fraction
1008	! so in turn on ice thickness
1009	END DO
1010	!zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d
1011
1012	!------------------------------------------------------------------------
1013	! Drainage through brine network (permeability)
1014	!------------------------------------------------------------------------
1015	!!! drainage due to ice permeability - Darcy's law
1016
1017	! sea water level
1018	msno = 0._wp
1019	DO jl = 1 , jpl
1020	msno = msno + v_s(ji,jj,jl) * rhos
1021	END DO
1022	floe_weight = ( msno + rhoivt_i(ji,jj) + rau0zvolp(ji,jj) ) / at_i(ji,jj)
1023	hsl_rel = floe_weight / rau0 &
1024	- ( ( sum(betan(:)*a_i(ji,jj,:)) / at_i(ji,jj) ) + alfan(1) )
1025
1026	deltah = hpond - hsl_rel
1027	pressure_head = grav * rau0 * max(deltah, 0._wp)
1028
1029	! drain if ice is permeable
1030	permflag = 0
1031
1032	IF (pressure_head > 0._wp) THEN
1033	DO jl = 1, jpl-1
1034	IF ( hicen(jl) /= 0._wp ) THEN
1035
1036	!IF (hicen(jl) > 0._wp) THEN !js: from CICE 5.1.2
1037
1038	perm = 0._wp ! MV ugly dummy patch
1039	CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl), sz_i(ji,jj,:,jl), perm) ! bof
1040	IF (perm > 0._wp) permflag = 1
1041
1042	drain = perma_ip(ji,jj,jl)pressure_head*rdt_ice / &
1043	(viscosity*hicen(jl))
1044	zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj))
1045	zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp)
1046
1047	diag_dvpn_drn(ji,jj) = - drain ! diag (could be better coded)
1048
1049	IF (zvolp(ji,jj) < epsi10) THEN
1050	zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj)
1051	zvolp(ji,jj) = 0._wp
1052	END IF
1053	END IF
1054	END DO
1055
1056	! adjust melt pond dimensions
1057	IF (permflag > 0) THEN
1058	! recompute pond depth
1059	CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index)
1060	DO jl = 1, m_index
1061	h_ip(ji,jj,jl) = hpond - alfan(jl) + alfan(1)
1062	a_ip(ji,jj,jl) = reduced_aicen(jl)
1063	END DO
1064	!zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d
1065	END IF
1066	END IF ! pressure_head
1067
1068	!-------------------------------
1069	! remove water from the snow
1070	!-------------------------------
1071	!------------------------------------------------------------------------
1072	! total melt pond volume in category does not include snow volume
1073	! snow in melt ponds is not melted
1074	!------------------------------------------------------------------------
1075
1076	! MV here, it seems that we remove some meltwater from the ponds, but I can't really tell
1077	! how much, so I did not diagnose it
1078	! so if there is a problem here, nobody is going to see it...
1079
1080
1081	! Calculate pond volume for lower categories
1082	DO jl = 1,m_index-1
1083	v_ip(ji,jj,jl) = a_ip(ji,jj,jl) * h_ip(ji,jj,jl) & ! what is not in the snow
1084	- (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl))
1085	END DO
1086
1087	! Calculate pond volume for highest category = remaining pond volume
1088
1089	! The following is completely unclear to Martin at least
1090	! Could we redefine properly and recode in a more readable way ?
1091
1092	! m_index = last category with melt pond
1093
1094	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
1095
1096	IF (m_index > 1) THEN
1097	IF (zvolp(ji,jj) > sum( v_ip(ji,jj,1:m_index-1))) THEN
1098	v_ip(ji,jj,m_index) = zvolp(ji,jj) - sum(v_ip(ji,jj,1:m_index-1))
1099	ELSE
1100	v_ip(ji,jj,m_index) = 0._wp
1101	h_ip(ji,jj,m_index) = 0._wp
1102	a_ip(ji,jj,m_index) = 0._wp
1103	! If remaining pond volume is negative reduce pond volume of
1104	! lower category
1105	IF ( zvolp(ji,jj) + epsi10 < SUM(v_ip(ji,jj,1:m_index-1))) &
1106	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)
1107	END IF
1108	END IF
1109
1110	DO jl = 1,m_index
1111	IF (a_ip(ji,jj,jl) > epsi10) THEN
1112	h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
1113	ELSE
1114	zdvolp(ji,jj) = zdvolp(ji,jj) + v_ip(ji,jj,jl)
1115	h_ip(ji,jj,jl) = 0._wp
1116	v_ip(ji,jj,jl) = 0._wp
1117	a_ip(ji,jj,jl) = 0._wp
1118	END IF
1119	END DO
1120	DO jl = m_index+1, jpl
1121	h_ip(ji,jj,jl) = 0._wp
1122	a_ip(ji,jj,jl) = 0._wp
1123	v_ip(ji,jj,jl) = 0._wp
1124	END DO
1125
1126	ENDIF
1127	END DO ! ji
1128	END DO ! jj
1129
1130	END SUBROUTINE ice_thd_pnd_area
1131
1132
1133	SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index)
1134	!!-------------------------------------------------------------------
1135	!! * ROUTINE ice_thd_pnd_depth *
1136	!!
1137	!! ** Purpose : Compute melt pond depth
1138	!!-------------------------------------------------------------------
1139
1140	REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
1141	aicen, &
1142	asnon, &
1143	hsnon, &
1144	alfan, &
1145	cum_max_vol
1146
1147	REAL (wp), INTENT(IN) :: &
1148	zvolp
1149
1150	REAL (wp), INTENT(OUT) :: &
1151	hpond
1152
1153	INTEGER, INTENT(OUT) :: &
1154	m_index
1155
1156	INTEGER :: n, ns
1157
1158	REAL (wp), DIMENSION(0:jpl+1) :: &
1159	hitl, &
1160	aicetl
1161
1162	REAL (wp) :: &
1163	rem_vol, &
1164	area, &
1165	vol, &
1166	tmp, &
1167	z0 = 0.0_wp
1168
1169	!----------------------------------------------------------------
1170	! hpond is zero if zvolp is zero - have we fully drained?
1171	!----------------------------------------------------------------
1172
1173	IF (zvolp < epsi10) THEN
1174	hpond = z0
1175	m_index = 0
1176	ELSE
1177
1178	!----------------------------------------------------------------
1179	! Calculate the category where water fills up to
1180	!----------------------------------------------------------------
1181
1182	!----------\|
1183	! \|
1184	! \|
1185	! \|----------\| -- --
1186	!__________\|__________\|_________________________________________ ^
1187	! \| \| rem_vol ^ \| Semi-filled
1188	! \| \|----------\|-- -- -- - ---\|-- ---- -- -- --v layer
1189	! \| \| \| \|
1190	! \| \| \| \|hpond
1191	! \| \| \|----------\| \| \|-------
1192	! \| \| \| \| \| \|
1193	! \| \| \| \|---v-----\|
1194	! \| \| m_index \| \| \|
1195	!-------------------------------------------------------------
1196
1197	m_index = 0 ! 1:m_index categories have water in them
1198	DO n = 1, jpl
1199	IF (zvolp <= cum_max_vol(n)) THEN
1200	m_index = n
1201	IF (n == 1) THEN
1202	rem_vol = zvolp
1203	ELSE
1204	rem_vol = zvolp - cum_max_vol(n-1)
1205	END IF
1206	exit ! to break out of the loop
1207	END IF
1208	END DO
1209	m_index = min(jpl-1, m_index)
1210
1211	!----------------------------------------------------------------
1212	! semi-filled layer may have m_index different snow in it
1213	!----------------------------------------------------------------
1214
1215	!----------------------------------------------------------- ^
1216	! \| alfan(m_index+1)
1217	! \|
1218	!hitl(3)--> \|----------\| \|
1219	!hitl(2)--> \|------------\| * * * * *\| \|
1220	!hitl(1)--> \|----------\|* * * * * * \|* * * * * \| \|
1221	!hitl(0)-->------------------------------------------------- \| ^
1222	! various snow from lower categories \| \|alfa(m_index)
1223
1224	! hitl - heights of the snow layers from thinner and current categories
1225	! aicetl - area of each snow depth in this layer
1226
1227	hitl(:) = z0
1228	aicetl(:) = z0
1229	DO n = 1, m_index
1230	hitl(n) = max(min(hsnon(n) + alfan(n) - alfan(m_index), &
1231	alfan(m_index+1) - alfan(m_index)), z0)
1232	aicetl(n) = asnon(n)
1233
1234	aicetl(0) = aicetl(0) + (aicen(n) - asnon(n))
1235	END DO
1236
1237	hitl(m_index+1) = alfan(m_index+1) - alfan(m_index)
1238	aicetl(m_index+1) = z0
1239
1240	!----------------------------------------------------------------
1241	! reorder array according to hitl
1242	! snow heights not necessarily in height order
1243	!----------------------------------------------------------------
1244
1245	DO ns = 1, m_index+1
1246	DO n = 0, m_index - ns + 1
1247	IF (hitl(n) > hitl(n+1)) THEN ! swap order
1248	tmp = hitl(n)
1249	hitl(n) = hitl(n+1)
1250	hitl(n+1) = tmp
1251	tmp = aicetl(n)
1252	aicetl(n) = aicetl(n+1)
1253	aicetl(n+1) = tmp
1254	END IF
1255	END DO
1256	END DO
1257
1258	!----------------------------------------------------------------
1259	! divide semi-filled layer into set of sublayers each vertically homogenous
1260	!----------------------------------------------------------------
1261
1262	!hitl(3)----------------------------------------------------------------
1263	! \| * * * * * * * *
1264	! \|* * * * * * * * *
1265	!hitl(2)----------------------------------------------------------------
1266	! \| * * * * * * * * \| * * * * * * * *
1267	! \|* * * * * * * * * \|* * * * * * * * *
1268	!hitl(1)----------------------------------------------------------------
1269	! \| * * * * * * * * \| * * * * * * * * \| * * * * * * * *
1270	! \|* * * * * * * * * \|* * * * * * * * * \|* * * * * * * * *
1271	!hitl(0)----------------------------------------------------------------
1272	! aicetl(0) aicetl(1) aicetl(2) aicetl(3)
1273
1274	! move up over layers incrementing volume
1275	DO n = 1, m_index+1
1276
1277	area = sum(aicetl(:)) - & ! total area of sub-layer
1278	(rhos/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow
1279
1280	vol = (hitl(n) - hitl(n-1)) * area ! thickness of sub-layer times area
1281
1282	IF (vol >= rem_vol) THEN ! have reached the sub-layer with the depth within
1283	hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1)
1284
1285	exit
1286	ELSE ! still in sub-layer below the sub-layer with the depth
1287	rem_vol = rem_vol - vol
1288	END IF
1289
1290	END DO
1291
1292	END IF
1293
1294	END SUBROUTINE ice_thd_pnd_depth
1295
1296
1297	SUBROUTINE ice_thd_pnd_perm(ticen, salin, perm)
1298	!!-------------------------------------------------------------------
1299	!! * ROUTINE ice_thd_pnd_perm *
1300	!!
1301	!! ** Purpose : Determine the liquid fraction of brine in the ice
1302	!! and its permeability
1303	!!-------------------------------------------------------------------
1304
1305	REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: &
1306	ticen, & ! internal ice temperature (K)
1307	salin ! salinity (ppt) !js: ppt according to cice
1308
1309	REAL (wp), INTENT(OUT) :: &
1310	perm ! permeability
1311
1312	REAL (wp) :: &
1313	Sbr ! brine salinity
1314
1315	REAL (wp), DIMENSION(nlay_i) :: &
1316	Tin, & ! ice temperature
1317	phi ! liquid fraction
1318
1319	INTEGER :: k
1320
1321	!-----------------------------------------------------------------
1322	! Compute ice temperatures from enthalpies using quadratic formula
1323	!-----------------------------------------------------------------
1324
1325	DO k = 1,nlay_i
1326	Tin(k) = ticen(k) - rt0 !js: from K to degC
1327	END DO
1328
1329	!-----------------------------------------------------------------
1330	! brine salinity and liquid fraction
1331	!-----------------------------------------------------------------
1332
1333	DO k = 1, nlay_i
1334
1335	! Sbr = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus)
1336	! Best expression to date is that one
1337	Sbr = - 18.7 * Tin(k) - 0.519 * Tin(k)*2 - 0.00535 Tin(k) **3
1338	phi(k) = salin(k) / Sbr
1339
1340	END DO
1341
1342	!-----------------------------------------------------------------
1343	! permeability
1344	!-----------------------------------------------------------------
1345
1346	perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007)
1347
1348	END SUBROUTINE ice_thd_pnd_perm
1349
1350
1351	!----------------------------------------------------------------------------------------------------------------------
1352
1353	SUBROUTINE ice_thd_pnd_init
1354	!!-------------------------------------------------------------------
1355	!! * ROUTINE ice_thd_pnd_init *
1356	!!
1357	!! ** Purpose : Physical constants and parameters linked to melt ponds
1358	!! over sea ice
1359	!!
1360	!! ** Method : Read the namthd_pnd namelist and check the melt pond
1361	!! parameter values called at the first timestep (nit000)
1362	!!
1363	!! ** input : Namelist namthd_pnd
1364	!!-------------------------------------------------------------------
1365	INTEGER :: ios, ioptio ! Local integer
1366	!!
1367	NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, rn_pnd_flush, &
1368	& ln_pnd_CST , rn_apnd, rn_hpnd, &
1369	& ln_pnd_TOPO , &
1370	& ln_pnd_lids, ln_pnd_alb
1371	!!-------------------------------------------------------------------
1372	!
1373	REWIND( numnam_ice_ref ) ! Namelist namthd_pnd in reference namelist : Melt Ponds
1374	READ ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901)
1375	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_pnd in reference namelist' )
1376	REWIND( numnam_ice_cfg ) ! Namelist namthd_pnd in configuration namelist : Melt Ponds
1377	READ ( numnam_ice_cfg, namthd_pnd, IOSTAT = ios, ERR = 902 )
1378	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namthd_pnd in configuration namelist' )
1379	IF(lwm) WRITE ( numoni, namthd_pnd )
1380	!
1381	IF(lwp) THEN ! control print
1382	WRITE(numout,*)
1383	WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
1384	WRITE(numout,*) '~~~~~~~~~~~~~~~~'
1385	WRITE(numout,*) ' Namelist namicethd_pnd:'
1386	WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd
1387	WRITE(numout,*) ' Topographic melt pond scheme ln_pnd_TOPO = ', ln_pnd_TOPO
1388	WRITE(numout,*) ' Level ice melt pond scheme ln_pnd_LEV = ', ln_pnd_LEV
1389	WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_apnd_min = ', rn_apnd_min
1390	WRITE(numout,*) ' Maximum ice fraction that contributes to melt ponds rn_apnd_max = ', rn_apnd_max
1391	WRITE(numout,*) ' Pond flushing efficiency rn_pnd_flush = ', rn_pnd_flush
1392	WRITE(numout,*) ' Constant ice melt pond scheme ln_pnd_CST = ', ln_pnd_CST
1393	WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd
1394	WRITE(numout,*) ' Prescribed pond depth rn_hpnd = ', rn_hpnd
1395	WRITE(numout,*) ' Frozen lids on top of melt ponds ln_pnd_lids = ', ln_pnd_lids
1396	WRITE(numout,*) ' Melt ponds affect albedo or not ln_pnd_alb = ', ln_pnd_alb
1397	ENDIF
1398	!
1399	! !== set the choice of ice pond scheme ==!
1400	ioptio = 0
1401	IF( .NOT.ln_pnd ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndNO ; ENDIF
1402	IF( ln_pnd_CST ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndCST ; ENDIF
1403	IF( ln_pnd_LEV ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndLEV ; ENDIF
1404	IF( ln_pnd_TOPO ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndTOPO ; ENDIF
1405	IF( ioptio /= 1 ) &
1406	& 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)' )
1407	!
1408	SELECT CASE( nice_pnd )
1409	CASE( np_pndNO )
1410	IF( ln_pnd_alb ) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF
1411	IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when no ponds' ) ; ENDIF
1412	CASE( np_pndCST )
1413	IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when constant ponds' ) ; ENDIF
1414	END SELECT
1415	!
1416	END SUBROUTINE ice_thd_pnd_init
1417
1418	#else
1419	!!----------------------------------------------------------------------
1420	!! Default option Empty module NO SI3 sea-ice model
1421	!!----------------------------------------------------------------------
1422	#endif
1423
1424	!!======================================================================
1425	END MODULE icethd_pnd

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

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

Context Navigation

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

Download in other formats: