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

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

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

Download in other formats: