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

Context Navigation

source: NEMO/branches/2020/SI3-05_MP/src/ICE/icethd_pnd.F90 @ 13809

Last change on this file since 13809 was 13809, checked in by vancop, 3 years ago
Original melt pond code from LLN
Property svn:keywords set to `Id`
File size: 56.9 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	!
23	USE in_out_manager ! I/O manager
24	USE lib_mpp ! MPP library
25	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
26	USE timing ! Timing
27
28	IMPLICIT NONE
29	PRIVATE
30
31	PUBLIC ice_thd_pnd_init ! routine called by icestp.F90
32	PUBLIC ice_thd_pnd ! routine called by icestp.F90
33
34	INTEGER :: nice_pnd ! choice of the type of pond scheme
35	! ! associated indices:
36	INTEGER, PARAMETER :: np_pndNO = 0 ! No pond scheme
37	INTEGER, PARAMETER :: np_pndCST = 1 ! Constant ice pond scheme
38	INTEGER, PARAMETER :: np_pndLEV = 2 ! Level ice pond scheme
39	INTEGER, PARAMETER :: np_pndTOPO = 3 ! Level ice pond scheme
40
41	!! * Substitutions
42	# include "vectopt_loop_substitute.h90"
43	!!----------------------------------------------------------------------
44	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
45	!! $Id$
46	!! Software governed by the CeCILL license (see ./LICENSE)
47	!!----------------------------------------------------------------------
48	CONTAINS
49
50	SUBROUTINE ice_thd_pnd
51	!!-------------------------------------------------------------------
52	!! * ROUTINE ice_thd_pnd *
53	!!
54	!! ** Purpose : change melt pond fraction and thickness
55	!!
56	!!-------------------------------------------------------------------
57	!
58	SELECT CASE ( nice_pnd )
59	!
60	CASE (np_pndCST) ; CALL pnd_CST !== Constant melt ponds ==!
61	!
62	CASE (np_pndLEV) ; CALL pnd_LEV !== Level ice melt ponds ==!
63	!
64	CASE (np_pndTOPO) ; CALL pnd_TOPO !== Topographic melt ponds ==!
65	!
66	END SELECT
67	!
68	END SUBROUTINE ice_thd_pnd
69
70
71	SUBROUTINE pnd_CST
72	!!-------------------------------------------------------------------
73	!! * ROUTINE pnd_CST *
74	!!
75	!! ** Purpose : Compute melt pond evolution
76	!!
77	!! ** Method : Melt pond fraction and thickness are prescribed
78	!! to non-zero values when t_su = 0C
79	!!
80	!! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd)
81	!!
82	!! ** Note : Coupling with such melt ponds is only radiative
83	!! Advection, ridging, rafting... are bypassed
84	!!
85	!! ** References : Bush, G.W., and Trump, D.J. (2017)
86	!!-------------------------------------------------------------------
87	INTEGER :: ji ! loop indices
88	!!-------------------------------------------------------------------
89	DO ji = 1, npti
90	!
91	IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN
92	h_ip_1d(ji) = rn_hpnd
93	a_ip_1d(ji) = rn_apnd * a_i_1d(ji)
94	h_il_1d(ji) = 0._wp ! no pond lids whatsoever
95	ELSE
96	h_ip_1d(ji) = 0._wp
97	a_ip_1d(ji) = 0._wp
98	h_il_1d(ji) = 0._wp
99	ENDIF
100	!
101	END DO
102	!
103	END SUBROUTINE pnd_CST
104
105
106	SUBROUTINE pnd_LEV
107	!!-------------------------------------------------------------------
108	!! * ROUTINE pnd_LEV *
109	!!
110	!! ** Purpose : Compute melt pond evolution
111	!!
112	!! ** Method : A fraction of meltwater is accumulated in ponds and sent to ocean when surface is freezing
113	!! We work with volumes and then redistribute changes into thickness and concentration
114	!! assuming linear relationship between the two.
115	!!
116	!! ** Action : - pond growth: Vp = Vp + dVmelt --- from Holland et al 2012 ---
117	!! dVmelt = (1-r)/rhow * ( rhoidh_i + rhosdh_s ) * a_i
118	!! dh_i = meltwater from ice surface melt
119	!! dh_s = meltwater from snow melt
120	!! (1-r) = fraction of melt water that is not flushed
121	!!
122	!! - limtations: a_ip must not exceed (1-r)*a_i
123	!! h_ip must not exceed 0.5*h_i
124	!!
125	!! - pond shrinking:
126	!! if lids: Vp = Vp -dH * a_ip
127	!! dH = lid thickness change. Retrieved from this eq.: --- from Flocco et al 2010 ---
128	!!
129	!! rhoi * Lf * dH/dt = ki * MAX(Tp-Tsu,0) / H
130	!! H = lid thickness
131	!! Lf = latent heat of fusion
132	!! Tp = -2C
133	!!
134	!! And solved implicitely as:
135	!! H(t+dt)*2 -H(t) H(t+dt) -ki * (Tp-Tsu) * dt / (rhoi*Lf) = 0
136	!!
137	!! if no lids: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp) --- from Holland et al 2012 ---
138	!!
139	!! - Flushing: w = -perm/visc * rho_oce * grav * Hp / Hi --- from Flocco et al 2007 ---
140	!! perm = permability of sea-ice
141	!! visc = water viscosity
142	!! Hp = height of top of the pond above sea-level
143	!! Hi = ice thickness thru which there is flushing
144	!!
145	!! - Corrections: remove melt ponds when lid thickness is 10 times the pond thickness
146	!!
147	!! - pond thickness and area is retrieved from pond volume assuming a linear relationship between h_ip and a_ip:
148	!! a_ip/a_i = a_ip_frac = h_ip / zaspect
149	!!
150	!! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min
151	!!
152	!! ** Note : mostly stolen from CICE
153	!!
154	!! ** References : Flocco and Feltham (JGR, 2007)
155	!! Flocco et al (JGR, 2010)
156	!! Holland et al (J. Clim, 2012)
157	!!-------------------------------------------------------------------
158
159	REAL(wp), DIMENSION(nlay_i) :: ztmp ! temporary array
160	!!
161	REAL(wp), PARAMETER :: zaspect = 0.8_wp ! pond aspect ratio
162	REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature
163	REAL(wp), PARAMETER :: zvisc = 1.79e-3_wp ! water viscosity
164	!!
165	REAL(wp) :: zfr_mlt, zdv_mlt ! fraction and volume of available meltwater retained for melt ponding
166	REAL(wp) :: zdv_frz, zdv_flush ! Amount of melt pond that freezes, flushes
167	REAL(wp) :: zhp ! heigh of top of pond lid wrt ssh
168	REAL(wp) :: zv_ip_max ! max pond volume allowed
169	REAL(wp) :: zdT ! zTp-t_su
170	REAL(wp) :: zsbr ! Brine salinity
171	REAL(wp) :: zperm ! permeability of sea ice
172	REAL(wp) :: zfac, zdum ! temporary arrays
173	REAL(wp) :: z1_rhow, z1_aspect, z1_Tp ! inverse
174	!!
175	INTEGER :: ji, jk ! loop indices
176
177	!!-------------------------------------------------------------------
178
179	z1_rhow = 1._wp / rhow
180	z1_aspect = 1._wp / zaspect
181	z1_Tp = 1._wp / zTp
182
183	DO ji = 1, npti
184	! !----------------------------------------------------!
185	IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN ! Case ice thickness < rn_himin or tiny ice fraction !
186	! !----------------------------------------------------!
187	!--- Remove ponds on thin ice or tiny ice fractions
188	a_ip_1d(ji) = 0._wp
189	h_ip_1d(ji) = 0._wp
190	h_il_1d(ji) = 0._wp
191	! !--------------------------------!
192	ELSE ! Case ice thickness >= rn_himin !
193	! !--------------------------------!
194	v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! retrieve volume from thickness
195	v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji)
196	!
197	!------------------!
198	! case ice melting !
199	!------------------!
200	!
201	!--- available meltwater for melt ponding ---!
202	zdum = -( dh_i_sum(ji)rhoi + dh_s_mlt(ji)rhos ) * z1_rhow * a_i_1d(ji)
203	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
204	zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! >0, max for roundoff errors?
205	!
206	!--- overflow ---!
207	! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume
208	! a_ip_max = zfr_mlt * a_i
209	! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:
210	zv_ip_max = zfr_mlt*2 a_i_1d(ji) * zaspect
211	zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
212
213	! If pond depth exceeds half the ice thickness then reduce the pond volume
214	! h_ip_max = 0.5 * h_i
215	! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:
216	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
217	zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
218
219	!--- Pond growing ---!
220	v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt
221	!
222	!--- Lid melting ---!
223	IF( ln_pnd_lids ) v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0
224
225	!
226	!--- mass flux ---!
227	IF( zdv_mlt > 0._wp ) THEN
228	! MV add comment on what that mass flux means
229	! water removed from fw flux due to melt pond growth
230	zfac = zdv_mlt * rhow * r1_rdtice ! melt pond mass flux < 0 [kg.m-2.s-1]
231	wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
232	!
233	! MV
234	! why surface melt and snow fluxes must be adjusted is not clear
235	! sounds like things are counted twice
236	!
237	zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) ! adjust ice/snow melting flux > 0 to balance melt pond flux
238	wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum)
239	wfx_sum_1d(ji) = wfx_sum_1d(ji) * (1._wp + zdum)
240	ENDIF
241
242	!-------------------!
243	! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0)
244	!-------------------!
245	!
246	zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp )
247	!
248	!--- Pond contraction (due to refreezing) ---!
249	IF( ln_pnd_lids ) THEN
250	!
251	!--- Lid growing and subsequent pond shrinking ---!
252	zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0
253	& SQRT( v_il_1d(ji)2 + a_ip_1d(ji)2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors
254
255	! Lid growing
256	v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz )
257
258	! Pond shrinking
259	v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz )
260
261	ELSE
262	! Pond shrinking
263	v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6)
264	ENDIF
265	!
266	!--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
267	! v_ip = h_ip * a_ip
268	! 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)
269	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
270	h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji)
271
272	!---------------!
273	! Pond flushing !
274	!---------------!
275	! height of top of the pond above sea-level
276	zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0
277
278	! Calculate the permeability of the ice (Assur 1958, see Flocco 2010)
279	! MV the expression for zsbr has wrong sign!!!!
280	! MV note - the sign-corrected expression is inconsistent
281	! with the rest of the SI3 code which assumes linear liquidus
282	! best expression (most consistent for SI3 should be)
283	! zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt
284	DO jk = 1, nlay_i
285	zsbr = - 1.2_wp &
286	& - 21.8_wp * ( t_i_1d(ji,jk) - rt0 ) &
287	& - 0.919_wp * ( t_i_1d(ji,jk) - rt0 )**2 &
288	& - 0.0178_wp * ( t_i_1d(ji,jk) - rt0 )**3
289
290	ztmp(jk) = sz_i_1d(ji,jk) / zsbr ! MV -> this is brine fraction and as it reads is always negative
291	END DO
292	zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) ! MV at present since ztmp is negative, this is always zero
293
294	! Do the drainage using Darcy's law
295	zdv_flush = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji)
296	zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) )
297	v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush
298
299	!--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
300	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
301	h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji)
302
303	!--- Corrections and lid thickness ---!
304	IF( ln_pnd_lids ) THEN
305	!--- retrieve lid thickness from volume ---!
306	IF( a_ip_1d(ji) > epsi10 ) THEN ; h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji)
307	ELSE ; h_il_1d(ji) = 0._wp
308	ENDIF
309	!--- remove ponds if lids are much larger than ponds ---!
310	IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN
311	a_ip_1d(ji) = 0._wp
312	h_ip_1d(ji) = 0._wp
313	h_il_1d(ji) = 0._wp
314	ENDIF
315	ENDIF
316	!
317	ENDIF
318
319	END DO
320	!
321	END SUBROUTINE pnd_LEV
322
323	SUBROUTINE pnd_TOPO (aice, aicen, &
324	vice, vicen, &
325	vsnon, &
326	ticen, salin, &
327	a_ip_frac, h_ip, &
328	Tsfc )
329
330	!!-------------------------------------------------------------------
331	!! * ROUTINE pnd_TOPO *
332	!!
333	!! ** Purpose : Compute melt pond evolution
334	!!
335	!! ** Purpose : Compute melt pond evolution based on the ice
336	!! topography as inferred from the ice thickness
337	!! distribution.
338	!!
339	!! ** Method : This code is initially based on Flocco and Feltham
340	!! (2007) and Flocco et al. (2010). More to come...
341	!!
342	!! ** Tunable parameters :
343	!!
344	!! ** Note :
345	!!
346	!! ** References
347	!! Flocco, D. and D. L. Feltham, 2007. A continuum model of melt pond
348	!! evolution on Arctic sea ice. J. Geophys. Res. 112, C08016, doi:
349	!! 10.1029/2006JC003836.
350	!! Flocco, D., D. L. Feltham and A. K. Turner, 2010. Incorporation of
351	!! a physically based melt pond scheme into the sea ice component of a
352	!! climate model. J. Geophys. Res. 115, C08012,
353	!! doi: 10.1029/2009JC005568.
354	!!
355	!!-------------------------------------------------------------------
356
357	!js 190423: the lid on melt ponds appears only in the analog subroutine of CICE 5.1.2
358
359	REAL (wp), DIMENSION (jpi,jpj), &
360	INTENT(IN) :: &
361	aice, & ! total ice area fraction
362	vice ! total ice volume (m)
363
364	REAL (wp), DIMENSION (jpi,jpj,jpl), &
365	INTENT(IN) :: &
366	aicen, & ! ice area fraction, per category
367	vsnon, & ! snow volume, per category (m)
368	vicen ! ice volume, per category (m)
369
370	REAL (wp), DIMENSION (jpi,jpj,nlay_i,jpl), &
371	INTENT(IN) :: &
372	ticen, & ! ice temperature per category (K)
373	salin
374
375	REAL (wp), DIMENSION (jpi,jpj,jpl), &
376	INTENT(INOUT) :: &
377	a_ip_frac , & ! pond area fraction of ice, per ice category
378	h_ip ! pond depth, per ice category (m)
379
380	REAL (wp), DIMENSION (jpi,jpj,jpl), &
381	INTENT(IN) :: &
382	Tsfc ! snow/sea ice surface temperature
383
384	! local variables
385	REAL (wp), DIMENSION (jpi,jpj,jpl) :: &
386	zTsfcn, & ! ice/snow surface temperature (C)
387	zvolpn ! pond volume per unit area, per category (m)
388
389	REAL (wp), DIMENSION (jpi,jpj,jpl) :: &
390	zrfrac, & ! fraction of available meltwater retained for melt ponding
391	zapondn,& ! pond area fraction, per category
392	zhpondn ! pond depth, per category (m)
393
394	REAL (wp), DIMENSION (jpi,jpj) :: &
395	zvolp, & ! total volume of pond, per unit area of pond (m)
396	zwfx_tmp ! temporary array for melt water
397
398	REAL (wp) :: &
399	zhi, & ! ice thickness (m)
400	zTavg, & ! mean surface temperature across categories (C)
401	z1_rhow, & ! inverse freshwater density
402	zdTs, & ! temperature difference for freeze-up (C)
403	zvpold, & ! dummy pond volume
404	zdvn ! change in melt pond volume for fresh water budget
405
406	INTEGER, DIMENSION (jpi*jpj) :: &
407	indxi, indxj ! compressed indices for cells with ice melting
408
409	INTEGER :: ij,icells,ji,jj,jl ! loop indices
410
411	REAL (wp), PARAMETER :: &
412	zr1_rlfus = 1._wp / 0.334e+6 / 917._wp , & ! (J/m^3)
413	zTp = -0.15_wp, & ! pond freezing temperature (C)
414	zmin_volp = 1.e-4_wp, & ! minimum pond volume (m)
415	zrexp = 0.01_wp, & ! constant melt pond freeze-up rate
416	z01 = 0.01_wp, &
417	z25 = 0.25_wp, &
418	z5 = 0.5_wp
419
420	z1_rhow = 1. / rhow
421
422	!---------------------------------------------------------------
423	! Initialization
424	!---------------------------------------------------------------
425	zhpondn (:,:,:) = 0._wp
426	zapondn (:,:,:) = 0._wp
427	zvolpn (:,:,:) = 0._wp
428
429	zTsfcn(:,:,:) = Tsfc(:,:,:) - rt0 ! Convert in Celsius
430
431	IF ( ln_pnd_fw ) THEN
432	v_ip_b(:,:,:) = v_ip(:,:,:)
433	ELSE
434	v_ip_b(:,:,:) = 0._wp
435	ENDIF
436
437	!------------------------------------------------------------------
438	! Available melt water for melt ponding and corresponding fraction
439	!------------------------------------------------------------------
440	!js 03/05/19 unset restriction on sign of wfx_pnd_in; mask values close to zero for future division
441	!wfx_pnd_in(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:), 0._wp ) ! available meltwater for melt ponding
442	!wfx_pnd_in(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:) , epsi10 ) * MAX( 0._wp, SIGN( 1._wp, wfx_sum(:,:) + wfx_snw_sum(:,:) - epsi10 ) )
443	!wfx_pnd_in(:,:) = (wfx_sum(:,:) + wfx_snw_sum(:,:)) * MAX( 0._wp, SIGN( 1._wp, ABS(wfx_sum(:,:) + wfx_snw_sum(:,:)) - epsi10 ) )
444
445	! MV
446	! NB: wfx_pnd_in can be slightly negative for very small values (why?)
447	! This can in some occasions give negative
448	! v_ip in the first category, which then gives crazy pond
449	! fractions and crashes the code as soon as the melt-pond
450	! radiative coupling is activated
451	! if we understand and remove why wfx_sum or wfx_snw could be
452	! negative, then, we can remove the MAX
453	! NB: I now changed to wfx_snw_sum, this may fix the problem.
454	! We should check
455
456
457	! OLI 07/2017: when we (MV & OL) first started the inclusion of melt
458	! ponds in the model, we removed the Holland et al. (2012, see
459	! CESM scheme above) parameterization. I put it back here,
460	! because I think it is needed. In summary, the sinks of FW for
461	! ponds are: 1/ Runoff through cracks/leads => depends on the
462	! total ice area only
463	! 2/ overflow, including Lüthje et al. (2006) limitation
464	! (max a_ip fraction function of h_i). This is in
465	! fact an other form of runoff that depends on the
466	! ITD
467	! 3/ Flushing - losses by permeability
468	! 4/ Refreezing
469	! 5/ Removal of ponds on thin ice
470	! I think 1 is needed because it is different from 2. However,
471	! test runs should/could be done, to check the sensitivity and
472	! the real usefulness of that stuff.
473	! Note : the Holland et al. param was wrongly wired in NEMO3.1 (using
474	! a_i instead of at_i), which might well explain why I had a too
475	! weak melt pond cover in my simulations (compared to MODIS, in
476	! situ obs. and CICE simulations.
477
478	!js 23/04/19: rewired back to a fraction with a_i
479	zrfrac(:,:,:) = rn_pnd_fracmin + ( rn_pnd_fracmax - rn_pnd_fracmin ) * aicen(:,:,:)
480	zwfx_tmp(:,:) = 0._wp
481
482	! MV ---> use expression from level-ice ponds, clarify what is needed
483
484	!--- Add retained melt water to melt ponds
485	! v_ip should never be negative, otherwise code crashes
486	! Rq MV: as far as I saw, UM5 can create very small negative v_ip values
487	! hence I added the max, which was not required with Prather (1 yr run)
488	! OLI: Here I use vt_ip, so I don't know if the max is
489	! required...
490	!zvolp(:,:) = MAX(vt_ip(:,:),0._wp) + zrfrac(:,:) * wfx_pnd_in(:,:) * z1_rhow * rdt_ice ! Total available melt water, to be distributed as melt ponds
491	! OLI: 07/2017 Bugfix above, removed " * aice(:,:)"
492	!js 19/04/18: change zrfrac to use aicen
493
494	zvolp(:,:) = vt_ip(:,:)
495
496	DO jl = 1, jpl
497	! Melt water, to be distributed as melt ponds
498	zvolp(:,:) = zvolp(:,:) - zrfrac(:,:,jl) &
499	* ( dh_i_pnd(:,:,jl)rhoic + dh_s_pnd(:,:,jl)rhosn ) &
500	* z1_rhow * a_i(:,:,jl)
501	END DO
502
503	! MV ---> use expression from level ice melt ponds (dv_mlt)
504
505	!js 03/05/19: we truncate negative values after calculating zvolp, in a
506	! similar manner to the subroutine lim_mp_cesm. Variation dh_i_pnd and
507	! dh_s_pnd are negative, indicating a loss of ice or snow. But we can expect them
508	! to be negative for some reasons. We keep this behaviour as it is, for
509	! fluxes conservation reasons. If some dh are positive, then we remove water
510	! indirectly from the ponds.
511	zvolp(:,:) = MAX( zvolp(:,:) , 0._wp )
512
513
514	! Fresh water flux going into the ponds
515	wfx_pnd_in(:,:) = wfx_pnd_in(:,:) + rhow * ( zvolp(:,:) - vt_ip(:,:) ) * r1_rdtice
516
517	!--- Remove retained meltwater from surface fluxes
518	IF ( ln_pnd_fw ) THEN
519	!wfx_snw_sum(:,:) = wfx_snw_sum(:,:) * ( 1. - zrfrac(:,:) )
520	!wfx_sum(:,:) = wfx_sum(:,:) * ( 1. - zrfrac(:,:) )
521
522	!js 190419: we change the code to use a_i in zrfrac. To be tested, but
523	!it should be conservative. zwfx_tmp is the flux accumulated in the
524	!ponds. wfx_pnd_in is the total surface melt fluxes.
525	zwfx_tmp(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:) , epsi10 ) &
526	& * MAX( 0._wp, SIGN( 1._wp, ABS(wfx_sum(:,:) + wfx_snw_sum(:,:)) - epsi10 ) )
527	WHERE ( ABS(zwfx_tmp(:,:)) > epsi10 )
528	zwfx_tmp(:,:) = wfx_pnd_in(:,:) / zwfx_tmp(:,:)
529	ELSEWHERE
530	zwfx_tmp(:,:) = 0._wp
531	ENDWHERE
532
533	wfx_sum(:,:) = ( 1._wp - zwfx_tmp(:,:) ) * wfx_sum(:,:)
534	wfx_snw_sum(:,:) = ( 1._wp - zwfx_tmp(:,:) ) * wfx_snw_sum(:,:)
535	ENDIF
536
537	!-----------------------------------------------------------------
538	! Identify grid cells with ponds
539	!-----------------------------------------------------------------
540
541	icells = 0
542	DO jj = 1, jpj
543	DO ji = 1, jpi
544	IF ( aice(ji,jj) > epsi10 ) THEN
545	zhi = vice(ji,jj) / aice(ji,jj)
546	ELSE
547	zhi = 0._wp
548	END IF
549
550	IF ( aice(ji,jj) > z01 .and. zhi > rn_himin .and. &
551	zvolp(ji,jj) > zmin_volp*aice(ji,jj)) THEN
552	icells = icells + 1
553	indxi(icells) = ji
554	indxj(icells) = jj
555	ELSE ! remove ponds on thin ice, or too small ponds
556	zvolpn(ji,jj,:) = 0._wp
557	zvolp (ji,jj) = 0._wp
558
559	a_ip(ji,jj,:) = 0._wp
560	v_ip(ji,jj,:) = 0._wp
561	a_ip_frac(ji,jj,:) = 0._wp
562	h_ip(ji,jj,:) = 0._wp
563
564	vt_ip(ji,jj) = 0._wp
565	at_ip(ji,jj) = 0._wp
566
567	! IF ( ln_pnd_fw ) & !--- Give freshwater to the ocean
568	! wfx_pnd_out(ji,jj) = wfx_pnd_out(ji,jj) + zvolp(ji,jj) * rhow * r1_rdtice
569	END IF
570	END DO ! ji
571	END DO ! jj
572
573
574	DO ij = 1, icells
575
576	ji = indxi(ij)
577	jj = indxj(ij)
578
579	!--------------------------------------------------------------
580	! Shrink pond due to refreezing
581	!--------------------------------------------------------------
582	! OLI 07/2017: Done like for empirical melt pond scheme (CESM).
583	! Therefore, I chose to put this part of the code before the main
584	! routines lim_mp_area/depth (contrary to the original code),
585	! seeing the freeze-up as a global sink of
586	! freshwater for melt ponds in the whole grid cell. If this was done
587	! after, I would need to make an additional assumption on the shape of
588	! melt ponds, which I don't want to do (for the CESM scheme, this
589	! assumption was on the aspect ratio). So I remove some water due to
590	! refreezing first (using zTavg instead of zTsfcn in each category) and
591	! then let the FF07 routines do their job for the fractional areas and
592	! depths of melt ponds.
593	! The whole ice lid related stuff from FF07 was thus removed and replaced
594	! by this. As mentionned below, this should be improved, but is much
595	! easier to conserve heat and freshwater this way.
596
597	! Average surface temperature is needed to compute freeze-up at the cell
598	! scale
599	zTavg = 0._wp
600	DO jl = 1, jpl
601	zTavg = zTavg + zTsfcn(ji,jj,jl)*aicen(ji,jj,jl)
602	END DO
603	zTavg = zTavg / aice(ji,jj)
604
605	! The freezing temperature for meltponds is assumed slightly below 0C,
606	! as if meltponds had a little salt in them (hence the use of zTp).
607	! The salt budget is not
608	! altered for meltponds, but if it were then an actual pond freezing
609	! temperature could be computed.
610
611	zdTs = MAX ( zTp - zTavg, 0. )
612
613	zvpold = zvolp(ji,jj)
614
615	zvolp(ji,jj) = zvolp(ji,jj) * EXP( zrexp * zdTs / zTp )
616
617	!--- Dump meltwater due to refreezing ( of course this is wrong
618	!--- but this parameterization is too simple )
619	! IF ( ln_pnd_fw ) &
620	! wfx_pnd_out(ji,jj) = wfx_pnd_out(ji,jj) + rhow * ( zvpold - zvolp(ji,jj) ) * r1_rdtice
621	! ! OLI 07/2017 : Bugfix above, zvpold - zvolp instead of the
622	! ! opposite, otherwise negative contribution
623
624	!--------------------------------------------------------------
625	! calculate pond area and depth
626	!--------------------------------------------------------------
627	zdvn = 0._wp
628	CALL lim_mp_area(aice(ji,jj),vice(ji,jj), &
629	aicen(ji,jj,:), vicen(ji,jj,:), vsnon(ji,jj,:), &
630	ticen(ji,jj,:,:), salin(ji,jj,:,:), &
631	zvolpn(ji,jj,:), zvolp(ji,jj), &
632	zapondn(ji,jj,:),zhpondn(ji,jj,:), zdvn)
633	! outputs are
634	! - zdvn
635	! - zvolpn
636	! - zvolp
637	! - zapondn
638	! - zhpondn
639
640	! IF ( ln_pnd_fw ) &
641	! wfx_pnd_out(ji,jj) = wfx_pnd_out(ji,jj) + zdvn * rhow * r1_rdtice ! update flux from ponds to ocean
642
643	!---------------------------------------------------------------
644	! Update pond volume and fraction
645	!---------------------------------------------------------------
646	DO jl = 1, jpl
647	a_ip(ji,jj,jl) = zapondn(ji,jj,jl)
648	v_ip(ji,jj,jl) = zvolpn(ji,jj,jl)
649	a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / MAX(aicen(ji,jj,jl), epsi10) &
650	* MAX(0._wp, SIGN(1._wp, aicen(ji,jj,jl) - epsi10))
651	h_ip (ji,jj,jl) = zhpondn(ji,jj,jl)
652	END DO
653	END DO ! ij
654
655	IF ( ln_pnd_fw ) THEN
656	!js 15/05/19: water going out of the ponds give a positive freshwater
657	! flux.
658	wfx_pnd_out(:,:) = SUM(MAX(0._wp, v_ip_b(:,:,:) - v_ip(:,:,:)), DIM=3) * rhow * r1_rdtice
659	ELSE
660	wfx_pnd_out(:,:) = 0._wp
661	ENDIF
662
663	END SUBROUTINE pnd_TOPO
664
665	SUBROUTINE lim_mp_area(aice, vice, &
666	aicen, vicen, vsnon, ticen, &
667	salin, zvolpn, zvolp, &
668	zapondn,zhpondn,dvolp)
669
670	!!-------------------------------------------------------------------
671	!! * ROUTINE lim_mp_area *
672	!!
673	!! ** Purpose : Given the total volume of meltwater, update
674	!! pond fraction (a_ip) and depth (should be volume)
675	!!
676	!! **
677	!!
678	!!------------------------------------------------------------------
679
680	REAL (wp), INTENT(IN) :: &
681	aice, vice
682
683	REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
684	aicen, vicen, vsnon
685
686	REAL (wp), DIMENSION(nlay_i,jpl), INTENT(IN) :: &
687	ticen, salin
688
689	REAL (wp), DIMENSION(jpl), INTENT(INOUT) :: &
690	zvolpn
691
692	REAL (wp), INTENT(INOUT) :: &
693	zvolp, dvolp
694
695	REAL (wp), DIMENSION(jpl), INTENT(OUT) :: &
696	zapondn, zhpondn
697
698	INTEGER :: &
699	n, ns, &
700	m_index, &
701	permflag
702
703	REAL (wp), DIMENSION(jpl) :: &
704	hicen, &
705	hsnon, &
706	asnon, &
707	rhos, & ! OLI 07/2017 : for now this is useless, but will be useful with new snow scheme
708	alfan, &
709	betan, &
710	cum_max_vol, &
711	reduced_aicen
712
713	REAL (wp), DIMENSION(0:jpl) :: &
714	cum_max_vol_tmp
715
716	REAL (wp) :: &
717	hpond, &
718	drain, &
719	floe_weight, &
720	pressure_head, &
721	hsl_rel, &
722	deltah, &
723	perm, &
724	msno
725
726	REAL (wp), parameter :: &
727	viscosity = 1.79e-3_wp, & ! kinematic water viscosity in kg/m/s
728	z0 = 0.0_wp , &
729	c1 = 1.0_wp , &
730	p4 = 0.4_wp , &
731	p6 = 0.6_wp
732
733	!-----------\|
734	! \|
735	! \|-----------\|
736	!___________\|___________\|______________________________________sea-level
737	! \| \|
738	! \| \|---^--------\|
739	! \| \| \| \|
740	! \| \| \| \|-----------\| \|-------
741	! \| \| \|alfan(n)\| \| \|
742	! \| \| \| \| \|--------------\|
743	! \| \| \| \| \| \|
744	!---------------------------v-------------------------------------------
745	! \| \| ^ \| \| \|
746	! \| \| \| \| \|--------------\|
747	! \| \| \|betan(n)\| \| \|
748	! \| \| \| \|-----------\| \|-------
749	! \| \| \| \|
750	! \| \|---v------- \|
751	! \| \|
752	! \|-----------\|
753	! \|
754	!-----------\|
755
756	!-------------------------------------------------------------------
757	! initialize
758	!-------------------------------------------------------------------
759
760	rhos(:) = rhosn ! OLI 07/2017 : same has above
761
762	DO n = 1, jpl
763
764	zapondn(n) = z0
765	zhpondn(n) = z0
766
767	!----------------------------------------
768	! X) compute the effective snow fraction
769	!----------------------------------------
770	IF (aicen(n) < epsi10) THEN
771	hicen(n) = z0
772	hsnon(n) = z0
773	reduced_aicen(n) = z0
774	asnon(n) = z0 !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables
775	ELSE
776	hicen(n) = vicen(n) / aicen(n)
777	hsnon(n) = vsnon(n) / aicen(n)
778	reduced_aicen(n) = c1 ! n=jpl
779
780	!js: initial code in NEMO_DEV
781	!IF (n < jpl) reduced_aicen(n) = aicen(n) &
782	! * (-0.024_wp*hicen(n) + 0.832_wp)
783
784	!js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large
785	IF (n < jpl) reduced_aicen(n) = aicen(n) &
786	* max(0.2_wp,(-0.024_wp*hicen(n) + 0.832_wp))
787
788	asnon(n) = reduced_aicen(n) ! effective snow fraction (empirical)
789	! MV should check whether this makes sense to have the same effective snow fraction in here
790	! OLI: it probably doesn't
791	END IF
792
793	! This choice for alfa and beta ignores hydrostatic equilibium of categories.
794	! Hydrostatic equilibium of the entire ITD is accounted for below, assuming
795	! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all
796	! categories. alfa and beta partition the ITD - they are areas not thicknesses!
797	! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area.
798	! Here, alfa = 60% of the ice area (and since hice is constant in a category,
799	! alfan = 60% of the ice volume) in each category lies above the reference line,
800	! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required.
801
802	! MV:
803	! Note that this choice is not in the original FF07 paper and has been adopted in CICE
804	! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe
805
806	! Where does that choice come from ? => OLI : Coz' Chuck Norris said so...
807
808	alfan(n) = 0.6 * hicen(n)
809	betan(n) = 0.4 * hicen(n)
810
811	cum_max_vol(n) = z0
812	cum_max_vol_tmp(n) = z0
813
814	END DO ! jpl
815
816	cum_max_vol_tmp(0) = z0
817	drain = z0
818	dvolp = z0
819
820	!----------------------------------------------------------
821	! x) Drain overflow water, update pond fraction and volume
822	!----------------------------------------------------------
823
824	!--------------------------------------------------------------------------
825	! the maximum amount of water that can be contained up to each ice category
826	!--------------------------------------------------------------------------
827
828	! MV
829	! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW)
830	! Then the excess volume cum_max_vol(jl) drains out of the system
831	! It should be added to wfx_pnd_out
832	! END MV
833	!js 18/04/19: XXX do something about this flux thing
834
835	DO n = 1, jpl-1 ! last category can not hold any volume
836
837	IF (alfan(n+1) >= alfan(n) .and. alfan(n+1) > z0) THEN
838
839	! total volume in level including snow
840	cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) + &
841	(alfan(n+1) - alfan(n)) * sum(reduced_aicen(1:n))
842
843	! subtract snow solid volumes from lower categories in current level
844	DO ns = 1, n
845	cum_max_vol_tmp(n) = cum_max_vol_tmp(n) &
846	- rhos(ns)/rhow * & ! free air fraction that can be filled by water
847	asnon(ns) * & ! effective areal fraction of snow in that category
848	max(min(hsnon(ns)+alfan(ns)-alfan(n), alfan(n+1)-alfan(n)), z0)
849	END DO
850
851	ELSE ! assume higher categories unoccupied
852	cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1)
853	END IF
854	!IF (cum_max_vol_tmp(n) < z0) THEN
855	! CALL abort_ice('negative melt pond volume')
856	!END IF
857	END DO
858	cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1) ! last category holds no volume
859	cum_max_vol (1:jpl) = cum_max_vol_tmp(1:jpl)
860
861	!----------------------------------------------------------------
862	! is there more meltwater than can be held in the floe?
863	!----------------------------------------------------------------
864	IF (zvolp >= cum_max_vol(jpl)) THEN
865	drain = zvolp - cum_max_vol(jpl) + epsi10
866	zvolp = zvolp - drain ! update meltwater volume available
867	dvolp = drain ! this is the drained water
868	IF (zvolp < epsi10) THEN
869	dvolp = dvolp + zvolp
870	zvolp = z0
871	END IF
872	END IF
873
874	! height and area corresponding to the remaining volume
875
876	CALL lim_mp_depth(reduced_aicen, asnon, hsnon, rhos, alfan, zvolp, cum_max_vol, hpond, m_index)
877
878	DO n=1, m_index
879	!zhpondn(n) = hpond - alfan(n) + alfan(1) ! here oui choulde update
880	! ! volume instead, no ?
881	zhpondn(n) = max((hpond - alfan(n) + alfan(1)), z0) !js: from CICE 5.1.2
882	zapondn(n) = reduced_aicen(n)
883	! in practise, pond fraction depends on the empirical snow fraction
884	! so in turn on ice thickness
885	END DO
886	!zapond = sum(zapondn(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d
887
888	!------------------------------------------------------------------------
889	! Drainage through brine network (permeability)
890	!------------------------------------------------------------------------
891	!!! drainage due to ice permeability - Darcy's law
892
893	! sea water level
894	msno = z0
895	DO n=1,jpl
896	msno = msno + vsnon(n) * rhos(n)
897	END DO
898	floe_weight = (msno + rhoicvice + rau0zvolp) / aice
899	hsl_rel = floe_weight / rau0 &
900	- ((sum(betan(:)*aicen(:))/aice) + alfan(1))
901
902	deltah = hpond - hsl_rel
903	pressure_head = grav * rau0 * max(deltah, z0)
904
905	! drain if ice is permeable
906	permflag = 0
907	IF (pressure_head > z0) THEN
908	DO n = 1, jpl-1
909	IF (hicen(n) /= z0) THEN
910	!IF (hicen(n) > z0) THEN !js: from CICE 5.1.2
911	perm = 0._wp ! MV ugly dummy patch
912	CALL lim_mp_perm(ticen(:,n), salin(:,n), perm)
913	IF (perm > z0) permflag = 1
914
915	drain = permzapondn(n)pressure_head*rdt_ice / &
916	(viscosity*hicen(n))
917	dvolp = dvolp + min(drain, zvolp)
918	zvolp = max(zvolp - drain, z0)
919	IF (zvolp < epsi10) THEN
920	dvolp = dvolp + zvolp
921	zvolp = z0
922	END IF
923	END IF
924	END DO
925
926	! adjust melt pond dimensions
927	IF (permflag > 0) THEN
928	! recompute pond depth
929	CALL lim_mp_depth(reduced_aicen, asnon, hsnon, rhos, alfan, zvolp, cum_max_vol, hpond, m_index)
930	DO n=1, m_index
931	zhpondn(n) = hpond - alfan(n) + alfan(1)
932	zapondn(n) = reduced_aicen(n)
933	END DO
934	!zapond = sum(zapondn(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d
935	END IF
936	END IF ! pressure_head
937
938	!-------------------------------
939	! X) remove water from the snow
940	!-------------------------------
941	!------------------------------------------------------------------------
942	! total melt pond volume in category does not include snow volume
943	! snow in melt ponds is not melted
944	!------------------------------------------------------------------------
945
946	! Calculate pond volume for lower categories
947	DO n=1,m_index-1
948	zvolpn(n) = zapondn(n) * zhpondn(n) & ! what is not in the snow
949	- (rhos(n)/rhow) * asnon(n) * min(hsnon(n), zhpondn(n))
950	END DO
951
952	! Calculate pond volume for highest category = remaining pond volume
953
954	! The following is completely unclear to Martin at least
955	! Could we redefine properly and recode in a more readable way ?
956
957	! m_index = last category with melt pond
958
959	IF (m_index == 1) zvolpn(m_index) = zvolp ! volume of mw in 1st category is the total volume of melt water
960
961	IF (m_index > 1) THEN
962	IF (zvolp > sum(zvolpn(1:m_index-1))) THEN
963	zvolpn(m_index) = zvolp - sum(zvolpn(1:m_index-1))
964	ELSE
965	zvolpn(m_index) = z0
966	zhpondn(m_index) = z0
967	zapondn(m_index) = z0
968	! If remaining pond volume is negative reduce pond volume of
969	! lower category
970	IF (zvolp+epsi10 < sum(zvolpn(1:m_index-1))) &
971	zvolpn(m_index-1) = zvolpn(m_index-1) - sum(zvolpn(1:m_index-1)) + zvolp
972	END IF
973	END IF
974
975	DO n=1,m_index
976	IF (zapondn(n) > epsi10) THEN
977	zhpondn(n) = zvolpn(n) / zapondn(n)
978	ELSE
979	dvolp = dvolp + zvolpn(n)
980	zhpondn(n) = z0
981	zvolpn(n) = z0
982	zapondn(n) = z0
983	END IF
984	END DO
985	DO n = m_index+1, jpl
986	zhpondn(n) = z0
987	zapondn(n) = z0
988	zvolpn (n) = z0
989	END DO
990
991	END SUBROUTINE lim_mp_area
992
993
994	SUBROUTINE lim_mp_depth(aicen, asnon, hsnon, rhos, alfan, zvolp, cum_max_vol, hpond, m_index)
995	!!-------------------------------------------------------------------
996	!! * ROUTINE lim_mp_depth *
997	!!
998	!! ** Purpose : Compute melt pond depth
999	!!-------------------------------------------------------------------
1000
1001	REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
1002	aicen, &
1003	asnon, &
1004	hsnon, &
1005	rhos, &
1006	alfan, &
1007	cum_max_vol
1008
1009	REAL (wp), INTENT(IN) :: &
1010	zvolp
1011
1012	REAL (wp), INTENT(OUT) :: &
1013	hpond
1014
1015	INTEGER, INTENT(OUT) :: &
1016	m_index
1017
1018	INTEGER :: n, ns
1019
1020	REAL (wp), DIMENSION(0:jpl+1) :: &
1021	hitl, &
1022	aicetl
1023
1024	REAL (wp) :: &
1025	rem_vol, &
1026	area, &
1027	vol, &
1028	tmp, &
1029	z0 = 0.0_wp
1030
1031	!----------------------------------------------------------------
1032	! hpond is zero if zvolp is zero - have we fully drained?
1033	!----------------------------------------------------------------
1034
1035	IF (zvolp < epsi10) THEN
1036	hpond = z0
1037	m_index = 0
1038	ELSE
1039
1040	!----------------------------------------------------------------
1041	! Calculate the category where water fills up to
1042	!----------------------------------------------------------------
1043
1044	!----------\|
1045	! \|
1046	! \|
1047	! \|----------\| -- --
1048	!__________\|__________\|_________________________________________ ^
1049	! \| \| rem_vol ^ \| Semi-filled
1050	! \| \|----------\|-- -- -- - ---\|-- ---- -- -- --v layer
1051	! \| \| \| \|
1052	! \| \| \| \|hpond
1053	! \| \| \|----------\| \| \|-------
1054	! \| \| \| \| \| \|
1055	! \| \| \| \|---v-----\|
1056	! \| \| m_index \| \| \|
1057	!-------------------------------------------------------------
1058
1059	m_index = 0 ! 1:m_index categories have water in them
1060	DO n = 1, jpl
1061	IF (zvolp <= cum_max_vol(n)) THEN
1062	m_index = n
1063	IF (n == 1) THEN
1064	rem_vol = zvolp
1065	ELSE
1066	rem_vol = zvolp - cum_max_vol(n-1)
1067	END IF
1068	exit ! to break out of the loop
1069	END IF
1070	END DO
1071	m_index = min(jpl-1, m_index)
1072
1073	!----------------------------------------------------------------
1074	! semi-filled layer may have m_index different snow in it
1075	!----------------------------------------------------------------
1076
1077	!----------------------------------------------------------- ^
1078	! \| alfan(m_index+1)
1079	! \|
1080	!hitl(3)--> \|----------\| \|
1081	!hitl(2)--> \|------------\| * * * * *\| \|
1082	!hitl(1)--> \|----------\|* * * * * * \|* * * * * \| \|
1083	!hitl(0)-->------------------------------------------------- \| ^
1084	! various snow from lower categories \| \|alfa(m_index)
1085
1086	! hitl - heights of the snow layers from thinner and current categories
1087	! aicetl - area of each snow depth in this layer
1088
1089	hitl(:) = z0
1090	aicetl(:) = z0
1091	DO n = 1, m_index
1092	hitl(n) = max(min(hsnon(n) + alfan(n) - alfan(m_index), &
1093	alfan(m_index+1) - alfan(m_index)), z0)
1094	aicetl(n) = asnon(n)
1095
1096	aicetl(0) = aicetl(0) + (aicen(n) - asnon(n))
1097	END DO
1098
1099	hitl(m_index+1) = alfan(m_index+1) - alfan(m_index)
1100	aicetl(m_index+1) = z0
1101
1102	!----------------------------------------------------------------
1103	! reorder array according to hitl
1104	! snow heights not necessarily in height order
1105	!----------------------------------------------------------------
1106
1107	DO ns = 1, m_index+1
1108	DO n = 0, m_index - ns + 1
1109	IF (hitl(n) > hitl(n+1)) THEN ! swap order
1110	tmp = hitl(n)
1111	hitl(n) = hitl(n+1)
1112	hitl(n+1) = tmp
1113	tmp = aicetl(n)
1114	aicetl(n) = aicetl(n+1)
1115	aicetl(n+1) = tmp
1116	END IF
1117	END DO
1118	END DO
1119
1120	!----------------------------------------------------------------
1121	! divide semi-filled layer into set of sublayers each vertically homogenous
1122	!----------------------------------------------------------------
1123
1124	!hitl(3)----------------------------------------------------------------
1125	! \| * * * * * * * *
1126	! \|* * * * * * * * *
1127	!hitl(2)----------------------------------------------------------------
1128	! \| * * * * * * * * \| * * * * * * * *
1129	! \|* * * * * * * * * \|* * * * * * * * *
1130	!hitl(1)----------------------------------------------------------------
1131	! \| * * * * * * * * \| * * * * * * * * \| * * * * * * * *
1132	! \|* * * * * * * * * \|* * * * * * * * * \|* * * * * * * * *
1133	!hitl(0)----------------------------------------------------------------
1134	! aicetl(0) aicetl(1) aicetl(2) aicetl(3)
1135
1136	! move up over layers incrementing volume
1137	DO n = 1, m_index+1
1138
1139	area = sum(aicetl(:)) - & ! total area of sub-layer
1140	(rhos(n)/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow
1141
1142	vol = (hitl(n) - hitl(n-1)) * area ! thickness of sub-layer times area
1143
1144	IF (vol >= rem_vol) THEN ! have reached the sub-layer with the depth within
1145	hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1)
1146
1147	exit
1148	ELSE ! still in sub-layer below the sub-layer with the depth
1149	rem_vol = rem_vol - vol
1150	END IF
1151
1152	END DO
1153
1154	END IF
1155
1156	END SUBROUTINE lim_mp_depth
1157
1158
1159	SUBROUTINE lim_mp_perm(ticen, salin, perm)
1160	!!-------------------------------------------------------------------
1161	!! * ROUTINE lim_mp_perm *
1162	!!
1163	!! ** Purpose : Determine the liquid fraction of brine in the ice
1164	!! and its permeability
1165	!!-------------------------------------------------------------------
1166
1167	REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: &
1168	ticen, & ! internal ice temperature (K)
1169	salin ! salinity (ppt) !js: ppt according to cice
1170
1171	REAL (wp), INTENT(OUT) :: &
1172	perm ! permeability
1173
1174	REAL (wp) :: &
1175	Sbr ! brine salinity
1176
1177	REAL (wp), DIMENSION(nlay_i) :: &
1178	Tin, & ! ice temperature
1179	phi ! liquid fraction
1180
1181	INTEGER :: k
1182
1183	!-----------------------------------------------------------------
1184	! Compute ice temperatures from enthalpies using quadratic formula
1185	!-----------------------------------------------------------------
1186
1187	DO k = 1,nlay_i
1188	Tin(k) = ticen(k) - rt0 !js: from K to degC
1189	END DO
1190
1191	!-----------------------------------------------------------------
1192	! brine salinity and liquid fraction
1193	!-----------------------------------------------------------------
1194
1195	IF (maxval(Tin) <= -2.0_wp) THEN
1196
1197	! Assur 1958
1198	DO k = 1,nlay_i
1199	Sbr = - 1.2_wp &
1200	- 21.8_wp * Tin(k) &
1201	- 0.919_wp * Tin(k)**2 &
1202	- 0.01878_wp * Tin(k)**3
1203	phi(k) = salin(k) / Sbr ! liquid fraction
1204
1205	!js: Sbr must not be zero
1206	END DO ! k
1207
1208	ELSE
1209
1210	! Notz 2005 thesis eq. 3.2
1211	!js: "interstitial brine Sbr (in ppt) as a function of temperature T (in degC)".
1212	DO k = 1,nlay_i
1213	Sbr = - 17.6_wp * Tin(k) &
1214	- 0.389_wp * Tin(k)**2 &
1215	- 0.00362_wp* Tin(k)**3
1216	phi(k) = salin(k) / Sbr ! liquid fraction
1217
1218	!js: Sbr must not be zero
1219	END DO
1220
1221	END IF
1222
1223	!-----------------------------------------------------------------
1224	! permeability
1225	!-----------------------------------------------------------------
1226
1227	perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007)
1228
1229	END SUBROUTINE lim_mp_perm
1230
1231
1232	!----------------------------------------------------------------------------------------------------------------------
1233
1234	SUBROUTINE ice_thd_pnd_init
1235	!!-------------------------------------------------------------------
1236	!! * ROUTINE ice_thd_pnd_init *
1237	!!
1238	!! ** Purpose : Physical constants and parameters linked to melt ponds
1239	!! over sea ice
1240	!!
1241	!! ** Method : Read the namthd_pnd namelist and check the melt pond
1242	!! parameter values called at the first timestep (nit000)
1243	!!
1244	!! ** input : Namelist namthd_pnd
1245	!!-------------------------------------------------------------------
1246	INTEGER :: ios, ioptio ! Local integer
1247	!!
1248	NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, &
1249	& ln_pnd_CST , rn_apnd, rn_hpnd, &
1250	& ln_pnd_TOPO , &
1251	& ln_pnd_lids, ln_pnd_alb
1252	!!-------------------------------------------------------------------
1253	!
1254	REWIND( numnam_ice_ref ) ! Namelist namthd_pnd in reference namelist : Melt Ponds
1255	READ ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901)
1256	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_pnd in reference namelist' )
1257	REWIND( numnam_ice_cfg ) ! Namelist namthd_pnd in configuration namelist : Melt Ponds
1258	READ ( numnam_ice_cfg, namthd_pnd, IOSTAT = ios, ERR = 902 )
1259	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namthd_pnd in configuration namelist' )
1260	IF(lwm) WRITE ( numoni, namthd_pnd )
1261	!
1262	IF(lwp) THEN ! control print
1263	WRITE(numout,*)
1264	WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
1265	WRITE(numout,*) '~~~~~~~~~~~~~~~~'
1266	WRITE(numout,*) ' Namelist namicethd_pnd:'
1267	WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd
1268	WRITE(numout,*) ' Topographic melt pond scheme ln_pnd_TOPO = ', ln_pnd_TOPO
1269	WRITE(numout,*) ' Level ice melt pond scheme ln_pnd_LEV = ', ln_pnd_LEV
1270	WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_apnd_min = ', rn_apnd_min
1271	WRITE(numout,*) ' Maximum ice fraction that contributes to melt ponds rn_apnd_max = ', rn_apnd_max
1272	WRITE(numout,*) ' Constant ice melt pond scheme ln_pnd_CST = ', ln_pnd_CST
1273	WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd
1274	WRITE(numout,*) ' Prescribed pond depth rn_hpnd = ', rn_hpnd
1275	WRITE(numout,*) ' Frozen lids on top of melt ponds ln_pnd_lids = ', ln_pnd_lids
1276	WRITE(numout,*) ' Melt ponds affect albedo or not ln_pnd_alb = ', ln_pnd_alb
1277	ENDIF
1278	!
1279	! !== set the choice of ice pond scheme ==!
1280	ioptio = 0
1281	IF( .NOT.ln_pnd ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndNO ; ENDIF
1282	IF( ln_pnd_CST ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndCST ; ENDIF
1283	IF( ln_pnd_LEV ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndLEV ; ENDIF
1284	IF( ioptio /= 1 ) &
1285	& CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_LEV or ln_pnd_CST)' )
1286	!
1287	SELECT CASE( nice_pnd )
1288	CASE( np_pndNO )
1289	IF( ln_pnd_alb ) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF
1290	IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when no ponds' ) ; ENDIF
1291	CASE( np_pndCST )
1292	IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when constant ponds' ) ; ENDIF
1293	END SELECT
1294	!
1295	END SUBROUTINE ice_thd_pnd_init
1296
1297	#else
1298	!!----------------------------------------------------------------------
1299	!! Default option Empty module NO SI3 sea-ice model
1300	!!----------------------------------------------------------------------
1301	#endif
1302
1303	!!======================================================================
1304	END MODULE icethd_pnd

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

Download in other formats: