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_martin_ponds/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/SI3_martin_ponds/src/ICE/icethd_pnd.F90 @ 13912

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

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

Download in other formats: