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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

icethd_pnd.F90 in NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/ICE – NEMO

Context Navigation

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

Last change on this file since 14451 was 14451, checked in by edblockley, 3 years ago
Removing duplicate definition of pond flushing logical that crept in as a bug with rev=14302
File size: 67.3 KB

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

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

Download in other formats: