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

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

icethd_pnd.F90 in NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icethd_pnd.F90 @ 14345

Last change on this file since 14345 was 14251, checked in by vancop, 4 years ago
Fix pond lids in topo ponds for christmas
Property svn:keywords set to `Id`
File size: 67.2 KB

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

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

Download in other formats: