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.

limmp.F90 in branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90 @ 8142

Last change on this file since 8142 was 8142, checked in by vancop, 7 years ago
Melt pond interfaces practically operational
File size: 92.8 KB

Line
1	MODULE limmp
2	!!======================================================================
3	!! * MODULE limmp *
4	!! Melt ponds
5	!!======================================================================
6	!! history : ! Original code by Daniela Flocco and Adrian Turner
7	!! 1.0 ! 2012 (O. Lecomte) Adaptation for scientific tests (NEMO3.1)
8	!! 2.0 ! 2017 (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation in NEMO3.6
9	!! ! NB: Only lim_mp_cstt and lim_mp_cesm work in this
10	!! version
11	!!----------------------------------------------------------------------
12	#if defined key_lim3
13	!!----------------------------------------------------------------------
14	!! 'key_lim3' : LIM3 sea-ice model
15	!!----------------------------------------------------------------------
16	!! lim_mp_init : some initialization and namelist read
17	!! lim_mp : main calling routine
18	!! lim_mp_topo : main melt pond routine for the "topographic" formulation (FloccoFeltham)
19	!! lim_mp_area : ??? compute melt pond fraction per category
20	!! lim_mp_perm : computes permeability (should be a FUNCTION!)
21	!! calc_hpond : computes melt pond depth
22	!! permeability_phy : computes permeability
23
24	!!----------------------------------------------------------------------
25	USE phycst ! physical constants
26	USE dom_oce ! ocean space and time domain
27	! USE sbc_ice ! Surface boundary condition: ice fields
28	USE ice ! LIM-3 variables
29	USE lbclnk ! lateral boundary conditions - MPP exchanges
30	USE lib_mpp ! MPP library
31	USE wrk_nemo ! work arrays
32	USE in_out_manager ! I/O manager
33	USE lib_fortran ! glob_sum
34	USE timing ! Timing
35	! USE limcons ! conservation tests
36	! USE limctl ! control prints
37	! USE limvar
38
39	!OLI_CODE USE ice_oce, ONLY: rdt_ice, tatm_ice
40	!OLI_CODE USE phycst
41	!OLI_CODE USE dom_ice
42	!OLI_CODE USE dom_oce
43	!OLI_CODE USE sbc_oce
44	!OLI_CODE USE sbc_ice
45	!OLI_CODE USE par_ice
46	!OLI_CODE USE par_oce
47	!OLI_CODE USE ice
48	!OLI_CODE USE thd_ice
49	!OLI_CODE USE in_out_manager
50	!OLI_CODE USE lbclnk
51	!OLI_CODE USE lib_mpp
52	!OLI_CODE
53	!OLI_CODE IMPLICIT NONE
54	!OLI_CODE PRIVATE
55	!OLI_CODE
56	!OLI_CODE PUBLIC lim_mp_init
57	!OLI_CODE PUBLIC lim_mp
58
59	IMPLICIT NONE
60	PRIVATE
61
62	PUBLIC lim_mp_init ! routine called by sbcice_lim.F90
63	PUBLIC lim_mp ! routine called by sbcice_lim.F90
64
65	!! * Substitutions
66	# include "vectopt_loop_substitute.h90"
67	!!----------------------------------------------------------------------
68	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
69	!! $Id: limdyn.F90 6994 2016-10-05 13:07:10Z clem $
70	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
71	!!----------------------------------------------------------------------
72	CONTAINS
73
74	SUBROUTINE lim_mp_init
75	!!-------------------------------------------------------------------
76	!! * ROUTINE lim_mp_init *
77	!!
78	!! ** Purpose : Physical constants and parameters linked to melt ponds
79	!! over sea ice
80	!!
81	!! ** Method : Read the namicemp namelist and check the melt pond
82	!! parameter values called at the first timestep (nit000)
83	!!
84	!! ** input : Namelist namicemp
85	!!-------------------------------------------------------------------
86	INTEGER :: ios ! Local integer output status for namelist read
87	NAMELIST/namicemp/ ln_pnd, ln_pnd_rad, ln_pnd_fw, nn_pnd_scheme, rn_apnd, rn_hpnd
88	!!-------------------------------------------------------------------
89
90	REWIND( numnam_ice_ref ) ! Namelist namicemp in reference namelist : Melt Ponds
91	READ ( numnam_ice_ref, namicemp, IOSTAT = ios, ERR = 901)
92	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicemp in reference namelist', lwp )
93
94	REWIND( numnam_ice_cfg ) ! Namelist namicemp in configuration namelist : Melt Ponds
95	READ ( numnam_ice_cfg, namicemp, IOSTAT = ios, ERR = 902 )
96	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicemp in configuration namelist', lwp )
97	IF(lwm) WRITE ( numoni, namicemp )
98
99	IF(lwp) THEN ! control print
100	WRITE(numout,*)
101	WRITE(numout,*) 'lim_mp_init : ice parameters for melt ponds'
102	WRITE(numout,*) '~~~~~~~~~~~~'
103	WRITE(numout,*) ' Active melt ponds ln_pnd = ', ln_pnd
104	WRITE(numout,*) ' Active melt ponds radiative coupling ln_pnd_rad = ', ln_pnd_rad
105	WRITE(numout,*) ' Active melt ponds freshwater coupling ln_pnd_fw = ', ln_pnd_fw
106	WRITE(numout,*) ' Type of melt pond scheme =0 presc, =1 empirical = 2 topo nn_pnd_scheme = ', nn_pnd_scheme
107	WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd
108	WRITE(numout,*) ' Prescribed pond depth rn_hpnd = ', rn_hpnd
109	ENDIF
110
111	IF ( .NOT. ln_pnd ) THEN
112	WRITE(numout,*)
113	WRITE(numout,*) ' Melt ponds are not activated '
114	WRITE(numout,*) ' ln_pnd_rad and ln_pnd_fw set to .FALSE. '
115	WRITE(numout,*) ' nn_pnd_scheme, rn_apnd, rn_hpnd set to zero '
116	ln_pnd_rad = .FALSE.
117	ln_pnd_fw = .FALSE.
118	nn_pnd_scheme = 0
119	rn_apnd = 0._wp
120	rn_hpnd = 0._wp
121
122	IF(lwp) THEN ! control print
123	WRITE(numout,*) ' Active melt ponds radiative coupling ln_pnd_rad = ', ln_pnd_rad
124	WRITE(numout,*) ' Active melt ponds freshwater coupling ln_pnd_fw = ', ln_pnd_fw
125	WRITE(numout,*) ' Type of melt pond scheme =0 presc, =1 empirical = 2 topo nn_pnd_scheme = ', nn_pnd_scheme
126	WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd
127	WRITE(numout,*) ' Prescribed pond depth rn_hpnd = ', rn_hpnd
128	ENDIF
129	ENDIF
130
131	IF ( ln_pnd .AND. ( nn_pnd_scheme == 2 ) .AND. ( jpl == 1 ) ) THEN
132	WRITE(numout,*) ' Topographic melt ponds are incompatible with jpl = 1 '
133	WRITE(numout,*) ' Run aborted '
134	CALL ctl_stop( 'STOP', 'lim_mp_init: uncompatible options, reset namelist_ice_ref ' )
135	ENDIF
136
137	!
138	END SUBROUTINE lim_mp_init
139
140
141	SUBROUTINE lim_mp( kt )
142	!!-------------------------------------------------------------------
143	!! * ROUTINE lim_mp *
144	!!
145	!! ** Purpose : change melt pond fraction
146	!!
147	!! ** Method : brutal force
148	!!
149	!! ** Action : -
150	!! -
151	!!------------------------------------------------------------------------------------
152
153	INTEGER, INTENT(in) :: kt ! number of iteration
154	INTEGER :: ji, jj, jl ! dummy loop indices
155
156	!!-------------------------------------------------------------------
157
158	IF( nn_timing == 1 ) CALL timing_start('lim_mp')
159
160	SELECT CASE ( nn_pnd_scheme )
161
162	CASE (0)
163
164	CALL lim_mp_cstt ! staircase melt ponds
165
166	CASE (1)
167
168	CALL lim_mp_cesm ! empirical melt ponds
169
170	CASE (2)
171
172	CALL lim_mp_topo & ! topographic melt ponds
173	& (at_i, a_i, &
174	& vt_i, v_i, v_s, t_i, s_i, a_ip_frac, &
175	& h_ip, t_su)
176
177	END SELECT
178
179	IF( nn_timing == 1 ) CALL timing_stop('lim_mp')
180
181	END SUBROUTINE lim_mp
182
183	SUBROUTINE lim_mp_cstt
184	!!-------------------------------------------------------------------
185	!! * ROUTINE lim_mp_cstt *
186	!!
187	!! ** Purpose : Compute melt pond evolution
188	!!
189	!! ** Method : Melt pond fraction and thickness are prescribed
190	!! to non-zero values when t_su = 0C
191	!!
192	!! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd)
193	!!
194	!! ** Note : Coupling with such melt ponds is only radiative
195	!! Advection, ridging, rafting... are bypassed
196	!!
197	!! ** References : Bush, G.W., and Trump, D.J. (2017)
198	!!
199	!!-------------------------------------------------------------------
200	INTEGER :: ji, jj, jl
201	REAL(wp) :: z1_jpl ! 1/jpl
202	!!-------------------------------------------------------------------
203
204	z1_jpl = 1. / FLOAT(jpl)
205
206	WHERE ( ( a_i > epsi10 ) .AND. ( t_su >= rt0-epsi06 ) )
207	a_ip_frac = rn_apnd
208	h_ip = rn_hpnd
209	v_ip = a_ip_frac * a_i * h_ip
210	a_ip = a_ip_frac * a_i
211	ELSE WHERE
212	a_ip = 0._wp
213	h_ip = 0._wp
214	v_ip = 0._wp
215	a_ip_frac = 0._wp
216	END WHERE
217
218	wfx_pnd(:,:) = 0._wp
219
220	END SUBROUTINE lim_mp_cstt
221
222	SUBROUTINE lim_mp_cesm
223	!!-------------------------------------------------------------------
224	!! * ROUTINE lim_mp_cesm *
225	!!
226	!! ** Purpose : Compute melt pond evolution
227	!!
228	!! ** Method : Empirical method. A fraction of meltwater is accumulated
229	!! in pond volume. It is then released exponentially when
230	!! surface is freezing.
231	!!
232	!! ** Tunable parameters : (no real expertise yet, ideas?)
233	!!
234	!! ** Note : Stolen from CICE for quick test of the melt pond
235	!! radiation and freshwater interfaces
236	!! Coupling can be radiative AND freshwater
237	!! Advection, ridging, rafting are called
238	!!
239	!! ** References : Holland, M. M. et al (J Clim 2012)
240	!!
241	!!-------------------------------------------------------------------
242
243	INTEGER, POINTER, DIMENSION(:) :: indxi ! compressed indices for cells with ice melting
244	INTEGER, POINTER, DIMENSION(:) :: indxj !
245
246	REAL(wp), POINTER, DIMENSION(:,:) :: zwfx_mlw ! available meltwater for melt ponding
247	REAL(wp), POINTER, DIMENSION(:,:,:) :: zrfrac ! fraction of available meltwater retained for melt ponding
248
249	REAL(wp), PARAMETER :: zrmin = 0.15_wp ! minimum fraction of available meltwater retained for melt ponding
250	REAL(wp), PARAMETER :: zrmax = 0.70_wp ! maximum '' '' '' '' ''
251	REAL(wp), PARAMETER :: zrexp = 0.01_wp ! rate constant to refreeze melt ponds
252	REAL(wp), PARAMETER :: zpnd_aspect = 0.8_wp ! pond aspect ratio
253
254	REAL(wp) :: zhi ! dummy ice thickness
255	REAL(wp) :: zhs ! dummy snow depth
256	REAL(wp) :: zTp ! reference temperature
257	REAL(wp) :: zdTs ! dummy temperature difference
258	REAL(wp) :: z1_rhofw ! inverse freshwater density
259	REAL(wp) :: z1_zpnd_aspect ! inverse pond aspect ratio
260	REAL(wp) :: zvpold ! dummy pond volume
261
262	INTEGER :: ji, jj, jl, ij ! loop indices
263	INTEGER :: icells ! size of dummy array
264
265	!!-------------------------------------------------------------------
266
267	CALL wrk_alloc( jpi*jpj, indxi, indxj)
268	CALL wrk_alloc( jpi,jpj, zwfx_mlw )
269	CALL wrk_alloc( jpi,jpj,jpl, zrfrac )
270
271	z1_rhofw = 1. / rhofw
272	z1_zpnd_aspect = 1. / zpnd_aspect
273	zTp = -2.
274
275	a_ip_frac(:,:,:) = 0._wp
276	h_ip (:,:,:) = 0._wp
277
278	!------------------------------------------------------------------
279	! Available melt water for melt ponding and corresponding fraction
280	!------------------------------------------------------------------
281
282	zwfx_mlw(:,:) = MAX( wfx_sum(:,:) + wfx_snw(:,:), 0._wp ) ! available meltwater for melt ponding
283
284	! NB: zwfx_mlw can be slightly negative for very small values (why?)
285	! This can in some occasions give negative
286	! v_ip in the first category, which then gives crazy pond
287	! fractions and crashes the code as soon as the melt-pond
288	! radiative coupling is activated
289	! if we understand and remove why wfx_sum or wfx_snw could be
290	! negative, then, we can remove the MAX
291
292	zrfrac(:,:,:) = zrmin + ( zrmax - zrmin ) * a_i(:,:,:)
293
294	DO jl = 1, jpl
295
296	!------------------------------------------------------------------------------
297	! Identify grid cells where ponds should be updated (can probably be improved)
298	!------------------------------------------------------------------------------
299
300	indxi(:) = 0
301	indxj(:) = 0
302	icells = 0
303
304	DO jj = 1, jpj
305	DO ji = 1, jpi
306	IF ( a_i(ji,jj,jl) > epsi10 ) THEN
307	icells = icells + 1
308	indxi(icells) = ji
309	indxj(icells) = jj
310	ENDIF
311	END DO ! ji
312	END DO ! jj
313
314	DO ij = 1, icells
315
316	ji = indxi(ij)
317	jj = indxj(ij)
318
319	zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
320	zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl)
321
322	IF ( zhi < rn_himin) THEN !--- Remove ponds on thin ice if ice is too thin
323
324	a_ip(ji,jj,jl) = 0._wp !--- Dump ponds
325	v_ip(ji,jj,jl) = 0._wp
326	a_ip_frac(ji,jj,jl) = 0._wp
327	h_ip(ji,jj,jl) = 0._wp
328
329	IF ( ln_pnd_fw ) & !--- Give freshwater to the ocean
330	wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + v_ip(ji,jj,jl)
331
332
333	ELSE !--- Update pond characteristics
334
335	!--- Add retained melt water to melt ponds
336	! v_ip should never be positive, otherwise code crashes
337	! MV: as far as I saw, UM5 can create very small negative v_ip values
338	! hence I added the max, which was not required with Prather (1 yr run)
339	v_ip(ji,jj,jl) = MAX( v_ip(ji,jj,jl), 0._wp ) + zrfrac(ji,jj,jl) * z1_rhofw * zwfx_mlw(ji,jj) * a_i(ji,jj,jl) * rdt_ice
340
341	!--- Shrink pond due to refreezing
342	zdTs = MAX ( zTp - t_su(ji,jj,jl) + rt0 , 0. )
343
344	zvpold = v_ip(ji,jj,jl)
345
346	v_ip(ji,jj,jl) = v_ip(ji,jj,jl) * EXP( zrexp * zdTs / zTp )
347
348	!--- Dump meltwater due to refreezing ( of course this is wrong
349	!--- but this parameterization is too simple )
350	IF ( ln_pnd_fw ) &
351	wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + rhofw * ( v_ip(ji,jj,jl) - zvpold ) * r1_rdtice
352
353	a_ip_frac(ji,jj,jl) = MIN( 1._wp , SQRT( v_ip(ji,jj,jl) * z1_zpnd_aspect / a_i(ji,jj,jl) ) )
354
355	h_ip(ji,jj,jl) = zpnd_aspect * a_ip_frac(ji,jj,jl)
356
357	a_ip(ji,jj,jl) = a_ip_frac(ji,jj,jl) * a_i(ji,jj,jl)
358
359	!-----------------------------------------------------------
360	! Limit pond depth
361	!-----------------------------------------------------------
362	! The original version has pond depth limitation, which I did not
363	! keep here. Maybe we need it later on
364	!
365
366	ENDIF
367
368	END DO
369
370	END DO ! jpl
371
372	!--- Remove retained meltwater from surface fluxes
373
374	IF ( ln_pnd_fw ) THEN
375
376	wfx_snw(:,:) = wfx_snw(:,:) * ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) )
377
378	wfx_sum(:,:) = wfx_sum(:,:) * ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) )
379
380	ENDIF
381
382	CALL wrk_dealloc( jpi*jpj, indxi, indxj)
383	CALL wrk_dealloc( jpi,jpj, zwfx_mlw )
384	CALL wrk_dealloc( jpi,jpj,jpl, zrfrac )
385
386	END SUBROUTINE lim_mp_cesm
387
388	SUBROUTINE lim_mp_topo (aice, aicen, &
389	vice, vicen, &
390	vsnon, &
391	ticen, salin, &
392	a_ip_frac, h_ip, &
393	Tsfc )
394	!!-------------------------------------------------------------------
395	!! * ROUTINE lim_mp_topo *
396	!!
397	!! ** Purpose : Compute melt pond evolution based on the ice
398	!! topography as inferred from the ice thickness
399	!! distribution.
400	!!
401	!! ** Method : This code is initially based on Flocco and Feltham
402	!! (2007) and Flocco et al. (2010). More to come...
403	!!
404	!! ** Tunable parameters :
405	!!
406	!! ** Note :
407	!!
408	!! ** References
409	!! Flocco, D. and D. L. Feltham, 2007. A continuum model of melt pond
410	!! evolution on Arctic sea ice. J. Geophys. Res. 112, C08016, doi:
411	!! 10.1029/2006JC003836.
412	!! Flocco, D., D. L. Feltham and A. K. Turner, 2010. Incorporation of
413	!! a physically based melt pond scheme into the sea ice component of a
414	!! climate model. J. Geophys. Res. 115, C08012,
415	!! doi: 10.1029/2009JC005568.
416	!!
417	!!-------------------------------------------------------------------
418
419	REAL (wp), DIMENSION (jpi,jpj), &
420	INTENT(IN) :: &
421	aice, & ! total ice area fraction
422	vice ! total ice volume (m)
423
424	REAL (wp), DIMENSION (jpi,jpj,jpl), &
425	INTENT(IN) :: &
426	aicen, & ! ice area fraction, per category
427	vsnon, & ! snow volume, per category (m)
428	vicen ! ice volume, per category (m)
429
430	REAL (wp), DIMENSION (jpi,jpj,nlay_i,jpl), &
431	INTENT(IN) :: &
432	ticen, & ! ice enthalpy, per category
433	salin
434
435	REAL (wp), DIMENSION (jpi,jpj,jpl), &
436	INTENT(INOUT) :: &
437	a_ip_frac , & ! pond area fraction of ice, per ice category
438	h_ip ! pond depth, per ice category (m)
439
440	REAL (wp), DIMENSION (jpi,jpj,jpl), &
441	INTENT(IN) :: &
442	Tsfc ! snow/sea ice surface temperature
443
444	! local variables
445	REAL (wp), DIMENSION (jpi,jpj,jpl) :: &
446	zTsfcn, & ! ice/snow surface temperature (C)
447	zvolpn, & ! pond volume per unit area, per category (m)
448	zvuin ! water-equivalent volume of ice lid on melt pond ('upper ice', m)
449
450	REAL (wp), DIMENSION (jpi,jpj,jpl) :: &
451	zapondn,& ! pond area fraction, per category
452	zhpondn ! pond depth, per category (m)
453
454	REAL (wp), DIMENSION (jpi,jpj) :: &
455	zvolp ! total volume of pond, per unit area of pond (m)
456
457	REAL (wp) :: &
458	zhi, & ! ice thickness (m)
459	zdHui, & ! change in thickness of ice lid (m)
460	zomega, & ! conduction
461	zdTice, & ! temperature difference across ice lid (C)
462	zdvice, & ! change in ice volume (m)
463	zTavg, & ! mean surface temperature across categories (C)
464	zTp, & ! pond freezing temperature (C)
465	zdvn ! change in melt pond volume for fresh water budget
466	INTEGER, DIMENSION (jpi*jpj) :: &
467	indxi, indxj ! compressed indices for cells with ice melting
468
469	INTEGER :: n,k,i,j,ij,icells,indxij ! loop indices
470
471	INTEGER, DIMENSION (jpl) :: &
472	kcells ! cells where ice lid combines with vice
473
474	INTEGER, DIMENSION (jpi*jpj,jpl) :: &
475	indxii, indxjj ! i,j indices for kcells loop
476
477	REAL (wp), parameter :: &
478	zhicemin = 0.1_wp , & ! minimum ice thickness with ponds (m)
479	zTd = 0.15_wp, & ! temperature difference for freeze-up (C)
480	zr1_rlfus = 1._wp / 0.334e+6 / 917._wp , & ! (J/m^3)
481	zmin_volp = 1.e-4_wp, & ! minimum pond volume (m)
482	z0 = 0._wp, &
483	zTimelt = 0._wp, &
484	z01 = 0.01_wp, &
485	z25 = 0.25_wp, &
486	z5 = 0.5_wp
487
488	!---------------------------------------------------------------
489	! Initialization
490	!---------------------------------------------------------------
491
492	zhpondn(:,:,:) = 0._wp
493	zapondn(:,:,:) = 0._wp
494	indxii(:,:) = 0
495	indxjj(:,:) = 0
496	kcells(:) = 0
497
498	zvolp(:,:) = wfx_sum(:,:) + wfx_snw(:,:) + vt_ip(:,:) ! Total available melt water, to be distributed as melt ponds
499	zTsfcn(:,:,:) = zTsfcn(:,:,:) - rt0 ! Convert in Celsius
500
501	! The freezing temperature for meltponds is assumed slightly below 0C,
502	! as if meltponds had a little salt in them. The salt budget is not
503	! altered for meltponds, but if it were then an actual pond freezing
504	! temperature could be computed.
505
506	! zTp = zTimelt - zTd ---> for lids
507
508	!-----------------------------------------------------------------
509	! Identify grid cells with ponds
510	!-----------------------------------------------------------------
511
512	icells = 0
513	DO j = 1, jpj
514	DO i = 1, jpi
515	zhi = z0
516	IF (aice(i,j) > epsi10 ) zhi = vice(i,j)/aice(i,j)
517	IF ( aice(i,j) > z01 .and. zhi > zhicemin .and. &
518	zvolp(i,j) > zmin_volp*aice(i,j)) THEN
519	icells = icells + 1
520	indxi(icells) = i
521	indxj(icells) = j
522	ELSE ! remove ponds on thin ice
523	!fpond(i,j) = fpond(i,j) - zvolp(i,j)
524	zvolpn(i,j,:) = z0
525	zvuin (i,j,:) = z0
526	zvolp (i,j) = z0
527	END IF
528	END DO ! i
529	END DO ! j
530
531	DO ij = 1, icells
532	i = indxi(ij)
533	j = indxj(ij)
534
535	!--------------------------------------------------------------
536	! calculate pond area and depth
537	!--------------------------------------------------------------
538	CALL lim_mp_area(aice(i,j),vice(i,j), &
539	aicen(i,j,:), vicen(i,j,:), vsnon(i,j,:), &
540	ticen(i,j,:,:), salin(i,j,:,:), &
541	zvolpn(i,j,:), zvolp(i,j), &
542	zapondn(i,j,:),zhpondn(i,j,:), zdvn)
543	! outputs are
544	! - zdvn
545	! - zvolpn
546	! - zvolp
547	! - zapondn
548	! - zhpondn
549
550	wfx_pnd(i,j) = wfx_pnd(i,j) + zdvn ! update flux from ponds to ocean
551
552	! mean surface temperature MV - why do we need that ? --> for the lid
553
554	! zTavg = z0
555	! DO n = 1, jpl
556	! zTavg = zTavg + zTsfcn(i,j,n)*aicen(i,j,n)
557	! END DO
558	! zTavg = zTavg / aice(i,j)
559
560	END DO ! ij
561
562	!---------------------------------------------------------------
563	! Update pond volume and fraction
564	!---------------------------------------------------------------
565
566	a_ip(:,:,:) = zapondn(:,:,:)
567	v_ip(:,:,:) = zapondn(:,:,:) * zhpondn(:,:,:)
568	a_ip_frac(:,:,:) = 0._wp
569	h_ip (:,:,:) = 0._wp
570
571	END SUBROUTINE lim_mp_topo
572
573	SUBROUTINE lim_mp_area(aice,vice, &
574	aicen, vicen, vsnon, ticen, &
575	salin, zvolpn, zvolp, &
576	zapondn,zhpondn,dvolp)
577
578	!!-------------------------------------------------------------------
579	!! * ROUTINE lim_mp_area *
580	!!
581	!! ** Purpose : Given the total volume of meltwater, update
582	!! pond fraction (a_ip) and depth (should be volume)
583	!!
584	!! **
585	!!
586	!!------------------------------------------------------------------
587
588	REAL (wp), INTENT(IN) :: &
589	aice,vice
590
591	REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
592	aicen, vicen, vsnon
593
594	REAL (wp), DIMENSION(nlay_i,jpl), INTENT(IN) :: &
595	ticen, salin
596
597	REAL (wp), DIMENSION(jpl), INTENT(INOUT) :: &
598	zvolpn
599
600	REAL (wp), INTENT(INOUT) :: &
601	zvolp, dvolp
602
603	REAL (wp), DIMENSION(jpl), INTENT(OUT) :: &
604	zapondn, zhpondn
605
606	INTEGER :: &
607	n, ns, &
608	m_index, &
609	permflag
610
611	REAL (wp), DIMENSION(jpl) :: &
612	hicen, &
613	hsnon, &
614	asnon, &
615	alfan, &
616	betan, &
617	cum_max_vol, &
618	reduced_aicen
619
620	REAL (wp), DIMENSION(0:jpl) :: &
621	cum_max_vol_tmp
622
623	REAL (wp) :: &
624	hpond, &
625	drain, &
626	floe_weight, &
627	pressure_head, &
628	hsl_rel, &
629	deltah, &
630	perm, &
631	msno
632
633	REAL (wp), parameter :: &
634	viscosity = 1.79e-3_wp, & ! kinematic water viscosity in kg/m/s
635	z0 = 0.0_wp , &
636	c1 = 1.0_wp , &
637	p4 = 0.4_wp , &
638	p6 = 0.6_wp , &
639	epsi10 = 1.0e-11_wp
640
641	!-----------\|
642	! \|
643	! \|-----------\|
644	!___________\|___________\|______________________________________sea-level
645	! \| \|
646	! \| \|---^--------\|
647	! \| \| \| \|
648	! \| \| \| \|-----------\| \|-------
649	! \| \| \|alfan(n)\| \| \|
650	! \| \| \| \| \|--------------\|
651	! \| \| \| \| \| \|
652	!---------------------------v-------------------------------------------
653	! \| \| ^ \| \| \|
654	! \| \| \| \| \|--------------\|
655	! \| \| \|betan(n)\| \| \|
656	! \| \| \| \|-----------\| \|-------
657	! \| \| \| \|
658	! \| \|---v------- \|
659	! \| \|
660	! \|-----------\|
661	! \|
662	!-----------\|
663
664	!-------------------------------------------------------------------
665	! initialize
666	!-------------------------------------------------------------------
667
668	DO n = 1, jpl
669
670	zapondn(n) = z0
671	zhpondn(n) = z0
672
673	!----------------------------------------
674	! X) compute the effective snow fraction
675	!----------------------------------------
676	IF (aicen(n) < epsi10) THEN
677	hicen(n) = z0
678	hsnon(n) = z0
679	reduced_aicen(n) = z0
680	ELSE
681	hicen(n) = vicen(n) / aicen(n)
682	hsnon(n) = vsnon(n) / aicen(n)
683	reduced_aicen(n) = c1 ! n=jpl
684	IF (n < jpl) reduced_aicen(n) = aicen(n) &
685	* (-0.024_wp*hicen(n) + 0.832_wp)
686	asnon(n) = reduced_aicen(n) ! effective snow fraction (empirical)
687	! MV should check whether this makes sense to have the same effective snow fraction in here
688	END IF
689
690	! This choice for alfa and beta ignores hydrostatic equilibium of categories.
691	! Hydrostatic equilibium of the entire ITD is accounted for below, assuming
692	! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all
693	! categories. alfa and beta partition the ITD - they are areas not thicknesses!
694	! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area.
695	! Here, alfa = 60% of the ice area (and since hice is constant in a category,
696	! alfan = 60% of the ice volume) in each category lies above the reference line,
697	! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required.
698
699	! MV:
700	! Note that this choice is not in the original FF07 paper and has been adopted in CICE
701	! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe
702
703	! Where does that choice come from
704
705	alfan(n) = 0.6 * hicen(n)
706	betan(n) = 0.4 * hicen(n)
707
708	cum_max_vol(n) = z0
709	cum_max_vol_tmp(n) = z0
710
711	END DO ! jpl
712
713	cum_max_vol_tmp(0) = z0
714	drain = z0
715	dvolp = z0
716
717	!----------------------------------------------------------
718	! x) Drain overflow water, update pond fraction and volume
719	!----------------------------------------------------------
720
721	!--------------------------------------------------------------------------
722	! the maximum amount of water that can be contained up to each ice category
723	!--------------------------------------------------------------------------
724
725	! MV
726	! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW)
727	! Then the excess volume cum_max_vol(jl) drains out of the system
728	! It should be added to wfx_pnd
729	! END MV
730
731	DO n = 1, jpl-1 ! last category can not hold any volume
732
733	IF (alfan(n+1) >= alfan(n) .and. alfan(n+1) > z0) THEN
734
735	! total volume in level including snow
736	cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) + &
737	(alfan(n+1) - alfan(n)) * sum(reduced_aicen(1:n))
738
739	! subtract snow solid volumes from lower categories in current level
740	DO ns = 1, n
741	cum_max_vol_tmp(n) = cum_max_vol_tmp(n) &
742	- rhosn/rhofw * & ! free air fraction that can be filled by water
743	asnon(ns) * & ! effective areal fraction of snow in that category
744	max(min(hsnon(ns)+alfan(ns)-alfan(n), alfan(n+1)- &
745	alfan(n)), z0)
746	END DO
747
748	ELSE ! assume higher categories unoccupied
749	cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1)
750	END IF
751	!IF (cum_max_vol_tmp(n) < z0) THEN
752	! call abort_ice('negative melt pond volume')
753	!END IF
754	END DO
755	cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1) ! last category holds no volume
756	cum_max_vol (1:jpl) = cum_max_vol_tmp(1:jpl)
757
758	!----------------------------------------------------------------
759	! is there more meltwater than can be held in the floe?
760	!----------------------------------------------------------------
761	IF (zvolp >= cum_max_vol(jpl)) THEN
762	drain = zvolp - cum_max_vol(jpl) + epsi10
763	zvolp = zvolp - drain ! update meltwater volume available
764	dvolp = drain ! this is the drained water
765	IF (zvolp < epsi10) THEN
766	dvolp = dvolp + zvolp
767	zvolp = z0
768	END IF
769	END IF
770
771	! height and area corresponding to the remaining volume
772
773	! call calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, &
774	! zvolp, cum_max_vol, hpond, m_index)
775
776	DO n=1, m_index
777	zhpondn(n) = hpond - alfan(n) + alfan(1) ! here oui choulde update
778	! volume instead, no ?
779	zapondn(n) = reduced_aicen(n)
780	! in practise, pond fraction depends on the empirical snow fraction
781	! so in turn on ice thickness
782	END DO
783
784	!------------------------------------------------------------------------
785	! Drainage through brine network (permeability)
786	!------------------------------------------------------------------------
787	!!! drainage due to ice permeability - Darcy's law
788
789	! sea water level
790	msno = z0
791	DO n=1,jpl
792	msno = msno + vsnon(n) * rhosn
793	END DO
794	floe_weight = (msno + rhoicvice + rau0zvolp) / aice
795	hsl_rel = floe_weight / rau0 &
796	- ((sum(betan(:)*aicen(:))/aice) + alfan(1))
797
798	deltah = hpond - hsl_rel
799	pressure_head = grav * rau0 * max(deltah, z0)
800
801	! drain IF ice is permeable
802	permflag = 0
803	IF (pressure_head > z0) THEN
804	DO n = 1, jpl-1
805	IF (hicen(n) /= z0) THEN
806	perm = 0. ! MV ugly dummy patch
807	CALL lim_mp_perm(ticen(:,n), salin(:,n), vicen(n), perm)
808	IF (perm > z0) permflag = 1
809
810	drain = permzapondn(n)pressure_head*rdt_ice / &
811	(viscosity*hicen(n))
812	dvolp = dvolp + min(drain, zvolp)
813	zvolp = max(zvolp - drain, z0)
814	IF (zvolp < epsi10) THEN
815	dvolp = dvolp + zvolp
816	zvolp = z0
817	END IF
818	END IF
819	END DO
820
821	! adjust melt pond DIMENSIONs
822	IF (permflag > 0) THEN
823	! recompute pond depth
824	! CALL calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, &
825	! zvolp, cum_max_vol, hpond, m_index)
826	DO n=1, m_index
827	zhpondn(n) = hpond - alfan(n) + alfan(1)
828	zapondn(n) = reduced_aicen(n)
829	END DO
830	END IF
831	END IF ! pressure_head
832
833	!-------------------------------
834	! X) remove water from the snow
835	!-------------------------------
836	!------------------------------------------------------------------------
837	! total melt pond volume in category DOes not include snow volume
838	! snow in melt ponds is not melted
839	!------------------------------------------------------------------------
840
841	! Calculate pond volume for lower categories
842	DO n=1,m_index-1
843	zvolpn(n) = zapondn(n) * zhpondn(n) & ! what is not in the snow
844	- (rhosn/rhofw) * asnon(n) * min(hsnon(n), zhpondn(n))
845	END DO
846
847	! Calculate pond volume for highest category = remaining pond volume
848
849	! The following is completely unclear to Martin at least
850	! Could we redefine properly and recode in a more readable way ?
851
852	! m_index = last category with melt pond
853
854	IF (m_index == 1) zvolpn(m_index) = zvolp ! volume of mw in 1st category is the total volume of melt water
855
856	IF (m_index > 1) THEN
857	IF (zvolp > sum(zvolpn(1:m_index-1))) THEN
858	zvolpn(m_index) = zvolp - sum(zvolpn(1:m_index-1)) !
859	ELSE
860	zvolpn(m_index) = z0
861	zhpondn(m_index) = z0
862	zapondn(m_index) = z0
863	! If remaining pond volume is negative reduce pond volume of
864	! lower category
865	IF (zvolp+epsi10 < sum(zvolpn(1:m_index-1))) &
866	zvolpn(m_index-1) = zvolpn(m_index-1)-sum(zvolpn(1:m_index-1))&
867	+ zvolp
868	END IF
869	END IF
870
871	DO n=1,m_index
872	IF (zapondn(n) > epsi10) THEN
873	zhpondn(n) = zvolpn(n) / zapondn(n)
874	ELSE
875	dvolp = dvolp + zvolpn(n)
876	zhpondn(n) = z0
877	zvolpn(n) = z0
878	zapondn(n) = z0
879	end IF
880	END DO
881	DO n = m_index+1, jpl
882	zhpondn(n) = z0
883	zapondn(n) = z0
884	zvolpn (n) = z0
885	END DO
886
887	END SUBROUTINE lim_mp_area
888
889	!OLI_CODE
890	!OLI_CODE
891	!OLI_CODE SUBROUTINE calc_hpond(aicen, asnon, hsnon, rhos, alfan, &
892	!OLI_CODE zvolp, cum_max_vol, &
893	!OLI_CODE hpond, m_index)
894	!OLI_CODE !!-------------------------------------------------------------------
895	!OLI_CODE !! * ROUTINE calc_hpond *
896	!OLI_CODE !!
897	!OLI_CODE !! ** Purpose : Compute melt pond depth
898	!OLI_CODE !!-------------------------------------------------------------------
899	!OLI_CODE
900	!OLI_CODE REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
901	!OLI_CODE aicen, &
902	!OLI_CODE asnon, &
903	!OLI_CODE hsnon, &
904	!OLI_CODE rhos, &
905	!OLI_CODE alfan, &
906	!OLI_CODE cum_max_vol
907	!OLI_CODE
908	!OLI_CODE REAL (wp), INTENT(IN) :: &
909	!OLI_CODE zvolp
910	!OLI_CODE
911	!OLI_CODE REAL (wp), INTENT(OUT) :: &
912	!OLI_CODE hpond
913	!OLI_CODE
914	!OLI_CODE INTEGER, INTENT(OUT) :: &
915	!OLI_CODE m_index
916	!OLI_CODE
917	!OLI_CODE INTEGER :: n, ns
918	!OLI_CODE
919	!OLI_CODE REAL (wp), DIMENSION(0:jpl+1) :: &
920	!OLI_CODE hitl, &
921	!OLI_CODE aicetl
922	!OLI_CODE
923	!OLI_CODE REAL (wp) :: &
924	!OLI_CODE rem_vol, &
925	!OLI_CODE area, &
926	!OLI_CODE vol, &
927	!OLI_CODE tmp, &
928	!OLI_CODE z0 = 0.0_wp, &
929	!OLI_CODE epsi10 = 1.0e-11_wp
930	!OLI_CODE
931	!OLI_CODE !----------------------------------------------------------------
932	!OLI_CODE ! hpond is zero if zvolp is zero - have we fully drained?
933	!OLI_CODE !----------------------------------------------------------------
934	!OLI_CODE
935	!OLI_CODE IF (zvolp < epsi10) THEN
936	!OLI_CODE hpond = z0
937	!OLI_CODE m_index = 0
938	!OLI_CODE ELSE
939	!OLI_CODE
940	!OLI_CODE !----------------------------------------------------------------
941	!OLI_CODE ! Calculate the category where water fills up to
942	!OLI_CODE !----------------------------------------------------------------
943	!OLI_CODE
944	!OLI_CODE !----------\|
945	!OLI_CODE ! \|
946	!OLI_CODE ! \|
947	!OLI_CODE ! \|----------\| -- --
948	!OLI_CODE !__________\|__________\|_________________________________________ ^
949	!OLI_CODE ! \| \| rem_vol ^ \| Semi-filled
950	!OLI_CODE ! \| \|----------\|-- -- -- - ---\|-- ---- -- -- --v layer
951	!OLI_CODE ! \| \| \| \|
952	!OLI_CODE ! \| \| \| \|hpond
953	!OLI_CODE ! \| \| \|----------\| \| \|-------
954	!OLI_CODE ! \| \| \| \| \| \|
955	!OLI_CODE ! \| \| \| \|---v-----\|
956	!OLI_CODE ! \| \| m_index \| \| \|
957	!OLI_CODE !-------------------------------------------------------------
958	!OLI_CODE
959	!OLI_CODE m_index = 0 ! 1:m_index categories have water in them
960	!OLI_CODE DO n = 1, jpl
961	!OLI_CODE IF (zvolp <= cum_max_vol(n)) THEN
962	!OLI_CODE m_index = n
963	!OLI_CODE IF (n == 1) THEN
964	!OLI_CODE rem_vol = zvolp
965	!OLI_CODE ELSE
966	!OLI_CODE rem_vol = zvolp - cum_max_vol(n-1)
967	!OLI_CODE END IF
968	!OLI_CODE exit ! to break out of the loop
969	!OLI_CODE END IF
970	!OLI_CODE END DO
971	!OLI_CODE m_index = min(jpl-1, m_index)
972	!OLI_CODE
973	!OLI_CODE !----------------------------------------------------------------
974	!OLI_CODE ! semi-filled layer may have m_index different snow in it
975	!OLI_CODE !----------------------------------------------------------------
976	!OLI_CODE
977	!OLI_CODE !----------------------------------------------------------- ^
978	!OLI_CODE ! \| alfan(m_index+1)
979	!OLI_CODE ! \|
980	!OLI_CODE !hitl(3)--> \|----------\| \|
981	!OLI_CODE !hitl(2)--> \|------------\| * * * * *\| \|
982	!OLI_CODE !hitl(1)--> \|----------\|* * * * * * \|* * * * * \| \|
983	!OLI_CODE !hitl(0)-->------------------------------------------------- \| ^
984	!OLI_CODE ! various snow from lower categories \| \|alfa(m_index)
985	!OLI_CODE
986	!OLI_CODE ! hitl - heights of the snow layers from thinner and current categories
987	!OLI_CODE ! aicetl - area of each snow depth in this layer
988	!OLI_CODE
989	!OLI_CODE hitl(:) = z0
990	!OLI_CODE aicetl(:) = z0
991	!OLI_CODE DO n = 1, m_index
992	!OLI_CODE hitl(n) = max(min(hsnon(n) + alfan(n) - alfan(m_index), &
993	!OLI_CODE alfan(m_index+1) - alfan(m_index)), z0)
994	!OLI_CODE aicetl(n) = asnon(n)
995	!OLI_CODE
996	!OLI_CODE aicetl(0) = aicetl(0) + (aicen(n) - asnon(n))
997	!OLI_CODE END DO
998	!OLI_CODE hitl(m_index+1) = alfan(m_index+1) - alfan(m_index)
999	!OLI_CODE aicetl(m_index+1) = z0
1000	!OLI_CODE
1001	!OLI_CODE !----------------------------------------------------------------
1002	!OLI_CODE ! reorder array according to hitl
1003	!OLI_CODE ! snow heights not necessarily in height order
1004	!OLI_CODE !----------------------------------------------------------------
1005	!OLI_CODE
1006	!OLI_CODE DO ns = 1, m_index+1
1007	!OLI_CODE DO n = 0, m_index - ns + 1
1008	!OLI_CODE IF (hitl(n) > hitl(n+1)) THEN ! swap order
1009	!OLI_CODE tmp = hitl(n)
1010	!OLI_CODE hitl(n) = hitl(n+1)
1011	!OLI_CODE hitl(n+1) = tmp
1012	!OLI_CODE tmp = aicetl(n)
1013	!OLI_CODE aicetl(n) = aicetl(n+1)
1014	!OLI_CODE aicetl(n+1) = tmp
1015	!OLI_CODE END IF
1016	!OLI_CODE END DO
1017	!OLI_CODE END DO
1018	!OLI_CODE
1019	!OLI_CODE !----------------------------------------------------------------
1020	!OLI_CODE ! divide semi-filled layer into set of sublayers each vertically homogenous
1021	!OLI_CODE !----------------------------------------------------------------
1022	!OLI_CODE
1023	!OLI_CODE !hitl(3)----------------------------------------------------------------
1024	!OLI_CODE ! \| * * * * * * * *
1025	!OLI_CODE ! \|* * * * * * * * *
1026	!OLI_CODE !hitl(2)----------------------------------------------------------------
1027	!OLI_CODE ! \| * * * * * * * * \| * * * * * * * *
1028	!OLI_CODE ! \|* * * * * * * * * \|* * * * * * * * *
1029	!OLI_CODE !hitl(1)----------------------------------------------------------------
1030	!OLI_CODE ! \| * * * * * * * * \| * * * * * * * * \| * * * * * * * *
1031	!OLI_CODE ! \|* * * * * * * * * \|* * * * * * * * * \|* * * * * * * * *
1032	!OLI_CODE !hitl(0)----------------------------------------------------------------
1033	!OLI_CODE ! aicetl(0) aicetl(1) aicetl(2) aicetl(3)
1034	!OLI_CODE
1035	!OLI_CODE ! move up over layers incrementing volume
1036	!OLI_CODE DO n = 1, m_index+1
1037	!OLI_CODE
1038	!OLI_CODE area = sum(aicetl(:)) - & ! total area of sub-layer
1039	!OLI_CODE (rhos(n)/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow
1040	!OLI_CODE
1041	!OLI_CODE vol = (hitl(n) - hitl(n-1)) * area ! thickness of sub-layer times area
1042	!OLI_CODE
1043	!OLI_CODE IF (vol >= rem_vol) THEN ! have reached the sub-layer with the depth within
1044	!OLI_CODE hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - &
1045	!OLI_CODE alfan(1)
1046	!OLI_CODE exit
1047	!OLI_CODE ELSE ! still in sub-layer below the sub-layer with the depth
1048	!OLI_CODE rem_vol = rem_vol - vol
1049	!OLI_CODE END IF
1050	!OLI_CODE
1051	!OLI_CODE END DO
1052	!OLI_CODE
1053	!OLI_CODE END IF
1054	!OLI_CODE
1055	!OLI_CODE END SUBROUTINE calc_hpond
1056	!OLI_CODE
1057	!OLI_CODE
1058	SUBROUTINE lim_mp_perm(ticen, salin, vicen, perm)
1059	!!-------------------------------------------------------------------
1060	!! * ROUTINE lim_mp_perm *
1061	!!
1062	!! ** Purpose : Determine the liquid fraction of brine in the ice
1063	!! and its permeability
1064	!!-------------------------------------------------------------------
1065	REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: &
1066	ticen, & ! energy of melting for each ice layer (J/m2)
1067	salin
1068
1069	REAL (wp), INTENT(IN) :: &
1070	vicen ! ice volume
1071
1072	REAL (wp), INTENT(OUT) :: &
1073	perm ! permeability
1074
1075	REAL (wp) :: &
1076	Sbr ! brine salinity
1077
1078	REAL (wp), DIMENSION(nlay_i) :: &
1079	Tin, & ! ice temperature
1080	phi ! liquid fraction
1081
1082	INTEGER :: k
1083
1084	REAL (wp) :: &
1085	c2 = 2.0_wp
1086
1087	!-----------------------------------------------------------------
1088	! Compute ice temperatures from enthalpies using quadratic formula
1089	!-----------------------------------------------------------------
1090
1091	DO k = 1,nlay_i
1092	Tin(k) = ticen(k)
1093	END DO
1094
1095	!-----------------------------------------------------------------
1096	! brine salinity and liquid fraction
1097	!-----------------------------------------------------------------
1098
1099	IF (maxval(Tin-rtt) <= -c2) THEN
1100
1101	DO k = 1,nlay_i
1102	Sbr = - 1.2_wp &
1103	-21.8_wp * (Tin(k)-rtt) &
1104	- 0.919_wp * (Tin(k)-rtt)**2 &
1105	- 0.01878_wp * (Tin(k)-rtt)**3
1106	phi(k) = salin(k)/Sbr ! liquid fraction
1107	END DO ! k
1108
1109	ELSE
1110
1111	DO k = 1,nlay_i
1112	Sbr = -17.6_wp * (Tin(k)-rtt) &
1113	- 0.389_wp * (Tin(k)-rtt)**2 &
1114	- 0.00362_wp* (Tin(k)-rtt)**3
1115	phi(k) = salin(k)/Sbr ! liquid fraction
1116	END DO
1117
1118	END IF
1119
1120	!-----------------------------------------------------------------
1121	! permeability
1122	!-----------------------------------------------------------------
1123
1124	perm = 3.0e-08_wp * (minval(phi))**3 ! REFERENCE PLEASE (this fucking
1125	! bastard of Golden)
1126
1127	END SUBROUTINE lim_mp_perm
1128	!OLI_CODE
1129	!OLI_CODE #else
1130	!OLI_CODE !!----------------------------------------------------------------------
1131	!OLI_CODE !! Default option Dummy Module No LIM-3 sea-ice model
1132	!OLI_CODE !!----------------------------------------------------------------------
1133	!OLI_CODE CONTAINS
1134	!OLI_CODE SUBROUTINE lim_mp_init ! Empty routine
1135	!OLI_CODE END SUBROUTINE lim_mp_init
1136	!OLI_CODE SUBROUTINE lim_mp ! Empty routine
1137	!OLI_CODE END SUBROUTINE lim_mp
1138	!OLI_CODE SUBROUTINE compute_mp_topo ! Empty routine
1139	!OLI_CODE END SUBROUTINE compute_mp_topo
1140	!OLI_CODE SUBROUTINE pond_area ! Empty routine
1141	!OLI_CODE END SUBROUTINE pond_area
1142	!OLI_CODE SUBROUTINE calc_hpond ! Empty routine
1143	!OLI_CODE END SUBROUTINE calc_hpond
1144	!OLI_CODE SUBROUTINE permeability_phy ! Empty routine
1145	!OLI_CODE END SUBROUTINE permeability_phy
1146	!OLI_CODE #endif
1147	!OLI_CODE !!======================================================================
1148	!OLI_CODE END MODULE limmp_topo
1149
1150
1151	!OLI_CODE MODULE limmp_topo
1152	!OLI_CODE !!======================================================================
1153	!OLI_CODE !! * MODULE limmp_topo *
1154	!OLI_CODE !! LIM-3 sea-ice : computation of melt ponds' properties
1155	!OLI_CODE !!======================================================================
1156	!OLI_CODE !! History : Original code by Daniela Flocco and Adrian Turner
1157	!OLI_CODE !! ! 2012-09 (O. Lecomte) Adaptation for routine inclusion in
1158	!OLI_CODE !! NEMO-LIM3.1
1159	!OLI_CODE !! ! 2016-11 (O. Lecomte, C. Rousset, M. Vancoppenolle)
1160	!OLI_CODE !! Adaptation for merge with NEMO-LIM3.6
1161	!OLI_CODE !!----------------------------------------------------------------------
1162	!OLI_CODE #if defined key_lim3
1163	!OLI_CODE !!----------------------------------------------------------------------
1164	!OLI_CODE !! 'key_lim3' LIM-3 sea-ice model
1165	!OLI_CODE !!----------------------------------------------------------------------
1166	!OLI_CODE !! lim_mp_init : melt pond properties initialization
1167	!OLI_CODE !! lim_mp : melt pond routine caller
1168	!OLI_CODE !! compute_mp_topo : Actual melt pond routine
1169	!OLI_CODE !!----------------------------------------------------------------------
1170	!OLI_CODE USE ice_oce, ONLY: rdt_ice, tatm_ice
1171	!OLI_CODE USE phycst
1172	!OLI_CODE USE dom_ice
1173	!OLI_CODE USE dom_oce
1174	!OLI_CODE USE sbc_oce
1175	!OLI_CODE USE sbc_ice
1176	!OLI_CODE USE par_ice
1177	!OLI_CODE USE par_oce
1178	!OLI_CODE USE ice
1179	!OLI_CODE USE thd_ice
1180	!OLI_CODE USE in_out_manager
1181	!OLI_CODE USE lbclnk
1182	!OLI_CODE USE lib_mpp
1183	!OLI_CODE
1184	!OLI_CODE IMPLICIT NONE
1185	!OLI_CODE PRIVATE
1186	!OLI_CODE
1187	!OLI_CODE PUBLIC lim_mp_init
1188	!OLI_CODE PUBLIC lim_mp
1189	!OLI_CODE
1190	!OLI_CODE CONTAINS
1191	!OLI_CODE
1192	!OLI_CODE SUBROUTINE lim_mp_init
1193	!OLI_CODE !!-------------------------------------------------------------------
1194	!OLI_CODE !! * ROUTINE lim_mp_init *
1195	!OLI_CODE !!
1196	!OLI_CODE !! ** Purpose : Initialize melt ponds
1197	!OLI_CODE !!-------------------------------------------------------------------
1198	!OLI_CODE a_ip_frac(:,:,:) = 0._wp
1199	!OLI_CODE a_ip(:,:,:) = 0._wp
1200	!OLI_CODE h_ip(:,:,:) = 0._wp
1201	!OLI_CODE v_ip(:,:,:) = 0._wp
1202	!OLI_CODE h_il(:,:,:) = 0._wp
1203	!OLI_CODE v_il(:,:,:) = 0._wp
1204	!OLI_CODE
1205	!OLI_CODE END SUBROUTINE lim_mp_init
1206	!OLI_CODE
1207	!OLI_CODE
1208	!OLI_CODE SUBROUTINE lim_mp
1209	!OLI_CODE !!-------------------------------------------------------------------
1210	!OLI_CODE !! * ROUTINE lim_mp *
1211	!OLI_CODE !!
1212	!OLI_CODE !! ** Purpose : Compute surface heat flux and call main melt pond
1213	!OLI_CODE !! routine
1214	!OLI_CODE !!-------------------------------------------------------------------
1215	!OLI_CODE
1216	!OLI_CODE INTEGER :: ji, jj, jl ! dummy loop indices
1217	!OLI_CODE
1218	!OLI_CODE fsurf(:,:) = 0.e0
1219	!OLI_CODE DO jl = 1, jpl
1220	!OLI_CODE DO jj = 1, jpj
1221	!OLI_CODE DO ji = 1, jpi
1222	!OLI_CODE fsurf(ji,jj) = fsurf(ji,jj) + a_i(ji,jj,jl) * &
1223	!OLI_CODE (qns_ice(ji,jj,jl) + (1.0 - izero(ji,jj,jl)) &
1224	!OLI_CODE * qsr_ice(ji,jj,jl))
1225	!OLI_CODE END DO
1226	!OLI_CODE END DO
1227	!OLI_CODE END DO
1228	!OLI_CODE
1229	!OLI_CODE CALL compute_mp_topo(at_i, a_i, &
1230	!OLI_CODE vt_i, v_i, v_s, rhosn_glo, t_i, s_i, a_ip_frac, &
1231	!OLI_CODE h_ip, h_il, t_su, tatm_ice, diag_sur_me*rdt_ice, &
1232	!OLI_CODE fsurf, fwoc)
1233	!OLI_CODE
1234	!OLI_CODE at_ip(:,:) = 0.0
1235	!OLI_CODE vt_ip(:,:) = 0.0
1236	!OLI_CODE vt_il(:,:) = 0.0
1237	!OLI_CODE DO jl = 1, jpl
1238	!OLI_CODE DO jj = 1, jpj
1239	!OLI_CODE DO ji = 1, jpi
1240	!OLI_CODE a_ip(ji,jj,jl) = MAX(0.0_wp, a_ip_frac(ji,jj,jl) * &
1241	!OLI_CODE a_i(ji,jj,jl))
1242	!OLI_CODE v_ip(ji,jj,jl) = MAX(0.0_wp, a_ip_frac(ji,jj,jl) * &
1243	!OLI_CODE a_i(ji,jj,jl) * h_ip(ji,jj,jl))
1244	!OLI_CODE v_il(ji,jj,jl) = MAX(0.0_wp, a_ip_frac(ji,jj,jl) * &
1245	!OLI_CODE a_i(ji,jj,jl) * h_il(ji,jj,jl))
1246	!OLI_CODE at_ip(ji,jj) = at_ip(ji,jj) + a_ip(ji,jj,jl)
1247	!OLI_CODE vt_ip(ji,jj) = vt_ip(ji,jj) + v_ip(ji,jj,jl)
1248	!OLI_CODE vt_il(ji,jj) = vt_il(ji,jj) + v_il(ji,jj,jl)
1249	!OLI_CODE END DO
1250	!OLI_CODE END DO
1251	!OLI_CODE END DO
1252	!OLI_CODE
1253	!OLI_CODE END SUBROUTINE lim_mp
1254	!OLI_CODE
1255	!OLI_CODE
1256	!OLI_CODE SUBROUTINE compute_mp_topo(aice, aicen, &
1257	!OLI_CODE vice, vicen, &
1258	!OLI_CODE vsnon, rhos, &
1259	!OLI_CODE ticen, salin, &
1260	!OLI_CODE a_ip_frac, h_ip, &
1261	!OLI_CODE h_il, Tsfc, &
1262	!OLI_CODE potT, meltt, &
1263	!OLI_CODE fsurf, fwoc)
1264	!OLI_CODE !!-------------------------------------------------------------------
1265	!OLI_CODE !! * ROUTINE compute_mp_topo *
1266	!OLI_CODE !!
1267	!OLI_CODE !! ** Purpose : Compute melt pond evolution based on the ice
1268	!OLI_CODE !! topography as inferred from the ice thickness
1269	!OLI_CODE !! distribution.
1270	!OLI_CODE !!
1271	!OLI_CODE !! ** Method : This code is initially based on Flocco and Feltham
1272	!OLI_CODE !! (2007) and Flocco et al. (2010). More to come...
1273	!OLI_CODE !!
1274	!OLI_CODE !! ** Tunable parameters :
1275	!OLI_CODE !!
1276	!OLI_CODE !! ** Note :
1277	!OLI_CODE !!
1278	!OLI_CODE !! ** References
1279	!OLI_CODE !! Flocco, D. and D. L. Feltham, 2007. A continuum model of melt pond
1280	!OLI_CODE !! evolution on Arctic sea ice. J. Geophys. Res. 112, C08016, doi:
1281	!OLI_CODE !! 10.1029/2006JC003836.
1282	!OLI_CODE !! Flocco, D., D. L. Feltham and A. K. Turner, 2010. Incorporation of
1283	!OLI_CODE !! a physically based melt pond scheme into the sea ice component of a
1284	!OLI_CODE !! climate model. J. Geophys. Res. 115, C08012,
1285	!OLI_CODE !! doi: 10.1029/2009JC005568.
1286	!OLI_CODE !!
1287	!OLI_CODE !!-------------------------------------------------------------------
1288	!OLI_CODE
1289	!OLI_CODE REAL (wp), DIMENSION (jpi,jpj), &
1290	!OLI_CODE INTENT(IN) :: &
1291	!OLI_CODE aice, & ! total ice area fraction
1292	!OLI_CODE vice ! total ice volume (m)
1293	!OLI_CODE
1294	!OLI_CODE REAL (wp), DIMENSION (jpi,jpj,jpl), &
1295	!OLI_CODE INTENT(IN) :: &
1296	!OLI_CODE aicen, & ! ice area fraction, per category
1297	!OLI_CODE vsnon, & ! snow volume, per category (m)
1298	!OLI_CODE rhos, & ! equivalent snow density, per category (kg/m^3)
1299	!OLI_CODE vicen ! ice volume, per category (m)
1300	!OLI_CODE
1301	!OLI_CODE REAL (wp), DIMENSION (jpi,jpj,nlay_i,jpl), &
1302	!OLI_CODE INTENT(IN) :: &
1303	!OLI_CODE ticen, & ! ice enthalpy, per category
1304	!OLI_CODE salin
1305	!OLI_CODE
1306	!OLI_CODE REAL (wp), DIMENSION (jpi,jpj,jpl), &
1307	!OLI_CODE INTENT(INOUT) :: &
1308	!OLI_CODE a_ip_frac , & ! pond area fraction of ice, per ice category
1309	!OLI_CODE h_ip , & ! pond depth, per ice category (m)
1310	!OLI_CODE h_il ! Refrozen ice lid thickness, per ice category (m)
1311	!OLI_CODE
1312	!OLI_CODE REAL (wp), DIMENSION (jpi,jpj), &
1313	!OLI_CODE INTENT(IN) :: &
1314	!OLI_CODE potT, & ! air potential temperature
1315	!OLI_CODE meltt, & ! total surface meltwater flux
1316	!OLI_CODE fsurf ! thermodynamic heat flux at ice/snow surface (W/m^2)
1317	!OLI_CODE
1318	!OLI_CODE REAL (wp), DIMENSION (jpi,jpj), &
1319	!OLI_CODE INTENT(INOUT) :: &
1320	!OLI_CODE fwoc ! fresh water flux to the ocean (from draining and other pond volume adjustments)
1321	!OLI_CODE ! (m)
1322	!OLI_CODE
1323	!OLI_CODE REAL (wp), DIMENSION (jpi,jpj,jpl), &
1324	!OLI_CODE INTENT(IN) :: &
1325	!OLI_CODE Tsfc ! snow/sea ice surface temperature
1326	!OLI_CODE
1327	!OLI_CODE ! local variables
1328	!OLI_CODE REAL (wp), DIMENSION (jpi,jpj,jpl) :: &
1329	!OLI_CODE zTsfcn, & ! ice/snow surface temperature (C)
1330	!OLI_CODE zvolpn, & ! pond volume per unit area, per category (m)
1331	!OLI_CODE zvuin ! water-equivalent volume of ice lid on melt pond ('upper ice', m)
1332	!OLI_CODE
1333	!OLI_CODE REAL (wp), DIMENSION (jpi,jpj,jpl) :: &
1334	!OLI_CODE zapondn,& ! pond area fraction, per category
1335	!OLI_CODE zhpondn ! pond depth, per category (m)
1336	!OLI_CODE
1337	!OLI_CODE REAL (wp), DIMENSION (jpi,jpj) :: &
1338	!OLI_CODE zvolp ! total volume of pond, per unit area of pond (m)
1339	!OLI_CODE
1340	!OLI_CODE REAL (wp) :: &
1341	!OLI_CODE zhi, & ! ice thickness (m)
1342	!OLI_CODE zdHui, & ! change in thickness of ice lid (m)
1343	!OLI_CODE zomega, & ! conduction
1344	!OLI_CODE zdTice, & ! temperature difference across ice lid (C)
1345	!OLI_CODE zdvice, & ! change in ice volume (m)
1346	!OLI_CODE zTavg, & ! mean surface temperature across categories (C)
1347	!OLI_CODE zTp, & ! pond freezing temperature (C)
1348	!OLI_CODE zdvn ! change in melt pond volume for fresh water budget
1349	!OLI_CODE INTEGER, DIMENSION (jpi*jpj) :: &
1350	!OLI_CODE indxi, indxj ! compressed indices for cells with ice melting
1351	!OLI_CODE
1352	!OLI_CODE INTEGER :: n,k,i,j,ij,icells,indxij ! loop indices
1353	!OLI_CODE
1354	!OLI_CODE INTEGER, DIMENSION (jpl) :: &
1355	!OLI_CODE kcells ! cells where ice lid combines with vice
1356	!OLI_CODE
1357	!OLI_CODE INTEGER, DIMENSION (jpi*jpj,jpl) :: &
1358	!OLI_CODE indxii, indxjj ! i,j indices for kcells loop
1359	!OLI_CODE
1360	!OLI_CODE REAL (wp), parameter :: &
1361	!OLI_CODE zhicemin = 0.1_wp , & ! minimum ice thickness with ponds (m)
1362	!OLI_CODE zTd = 0.15_wp, & ! temperature difference for freeze-up (C)
1363	!OLI_CODE zr1_rlfus = 1._wp / 0.334e+6 / 917._wp , & ! (J/m^3)
1364	!OLI_CODE zmin_volp = 1.e-4_wp, & ! minimum pond volume (m)
1365	!OLI_CODE z0 = 0._wp, &
1366	!OLI_CODE zTimelt = 0._wp, &
1367	!OLI_CODE z01 = 0.01_wp, &
1368	!OLI_CODE z25 = 0.25_wp, &
1369	!OLI_CODE z5 = 0.5_wp, &
1370	!OLI_CODE epsi10 = 1.0e-11_wp
1371	!OLI_CODE !---------------------------------------------------------------
1372	!OLI_CODE ! initialize
1373	!OLI_CODE !---------------------------------------------------------------
1374	!OLI_CODE
1375	!OLI_CODE DO j = 1, jpj
1376	!OLI_CODE DO i = 1, jpi
1377	!OLI_CODE zvolp(i,j) = z0
1378	!OLI_CODE END DO
1379	!OLI_CODE END DO
1380	!OLI_CODE DO n = 1, jpl
1381	!OLI_CODE DO j = 1, jpj
1382	!OLI_CODE DO i = 1, jpi
1383	!OLI_CODE ! load tracers
1384	!OLI_CODE zvolp(i,j) = zvolp(i,j) + h_ip(i,j,n) &
1385	!OLI_CODE * a_ip_frac(i,j,n) * aicen(i,j,n)
1386	!OLI_CODE zTsfcn(i,j,n) = Tsfc(i,j,n) - rtt ! convert in Celsius - Oli
1387	!OLI_CODE zvuin (i,j,n) = h_il(i,j,n) &
1388	!OLI_CODE * a_ip_frac(i,j,n) * aicen(i,j,n)
1389	!OLI_CODE
1390	!OLI_CODE zhpondn(i,j,n) = z0 ! pond depth, per category
1391	!OLI_CODE zapondn(i,j,n) = z0 ! pond area, per category
1392	!OLI_CODE END DO
1393	!OLI_CODE END DO
1394	!OLI_CODE indxii(:,n) = 0
1395	!OLI_CODE indxjj(:,n) = 0
1396	!OLI_CODE kcells (n) = 0
1397	!OLI_CODE END DO
1398	!OLI_CODE
1399	!OLI_CODE ! The freezing temperature for meltponds is assumed slightly below 0C,
1400	!OLI_CODE ! as if meltponds had a little salt in them. The salt budget is not
1401	!OLI_CODE ! altered for meltponds, but if it were then an actual pond freezing
1402	!OLI_CODE ! temperature could be computed.
1403	!OLI_CODE
1404	!OLI_CODE zTp = zTimelt - zTd
1405	!OLI_CODE
1406	!OLI_CODE !-----------------------------------------------------------------
1407	!OLI_CODE ! Identify grid cells with ponds
1408	!OLI_CODE !-----------------------------------------------------------------
1409	!OLI_CODE
1410	!OLI_CODE icells = 0
1411	!OLI_CODE DO j = 1, jpj
1412	!OLI_CODE DO i = 1, jpi
1413	!OLI_CODE zhi = z0
1414	!OLI_CODE IF (aice(i,j) > epsi10) zhi = vice(i,j)/aice(i,j)
1415	!OLI_CODE IF ( aice(i,j) > z01 .and. zhi > zhicemin .and. &
1416	!OLI_CODE zvolp(i,j) > zmin_volp*aice(i,j)) THEN
1417	!OLI_CODE icells = icells + 1
1418	!OLI_CODE indxi(icells) = i
1419	!OLI_CODE indxj(icells) = j
1420	!OLI_CODE ELSE ! remove ponds on thin ice
1421	!OLI_CODE !fpond(i,j) = fpond(i,j) - zvolp(i,j)
1422	!OLI_CODE zvolpn(i,j,:) = z0
1423	!OLI_CODE zvuin (i,j,:) = z0
1424	!OLI_CODE zvolp (i,j) = z0
1425	!OLI_CODE END IF
1426	!OLI_CODE END DO ! i
1427	!OLI_CODE END DO ! j
1428	!OLI_CODE
1429	!OLI_CODE DO ij = 1, icells
1430	!OLI_CODE i = indxi(ij)
1431	!OLI_CODE j = indxj(ij)
1432	!OLI_CODE
1433	!OLI_CODE !--------------------------------------------------------------
1434	!OLI_CODE ! calculate pond area and depth
1435	!OLI_CODE !--------------------------------------------------------------
1436	!OLI_CODE CALL pond_area(aice(i,j),vice(i,j),rhos(i,j,:), &
1437	!OLI_CODE aicen(i,j,:), vicen(i,j,:), vsnon(i,j,:), &
1438	!OLI_CODE ticen(i,j,:,:), salin(i,j,:,:), &
1439	!OLI_CODE zvolpn(i,j,:), zvolp(i,j), &
1440	!OLI_CODE zapondn(i,j,:),zhpondn(i,j,:), zdvn)
1441	!OLI_CODE
1442	!OLI_CODE fwoc(i,j) = fwoc(i,j) + zdvn ! -> Goes to fresh water budget
1443	!OLI_CODE
1444	!OLI_CODE ! mean surface temperature
1445	!OLI_CODE zTavg = z0
1446	!OLI_CODE DO n = 1, jpl
1447	!OLI_CODE zTavg = zTavg + zTsfcn(i,j,n)*aicen(i,j,n)
1448	!OLI_CODE END DO
1449	!OLI_CODE zTavg = zTavg / aice(i,j)
1450	!OLI_CODE
1451	!OLI_CODE DO n = 1, jpl-1
1452	!OLI_CODE
1453	!OLI_CODE IF (zvuin(i,j,n) > epsi10) THEN
1454	!OLI_CODE
1455	!OLI_CODE !----------------------------------------------------------------
1456	!OLI_CODE ! melting: floating upper ice layer melts in whole or part
1457	!OLI_CODE !----------------------------------------------------------------
1458	!OLI_CODE ! IF (zTsfcn(i,j,n) > zTp) THEN
1459	!OLI_CODE IF (zTavg > zTp) THEN
1460	!OLI_CODE
1461	!OLI_CODE zdvice = min(meltt(i,j)*zapondn(i,j,n), zvuin(i,j,n))
1462	!OLI_CODE IF (zdvice > epsi10) THEN
1463	!OLI_CODE zvuin (i,j,n) = zvuin (i,j,n) - zdvice
1464	!OLI_CODE zvolpn(i,j,n) = zvolpn(i,j,n) + zdvice
1465	!OLI_CODE zvolp (i,j) = zvolp (i,j) + zdvice
1466	!OLI_CODE !fwoc(i,j) = fwoc(i,j) + zdvice
1467	!OLI_CODE
1468	!OLI_CODE IF (zvuin(i,j,n) < epsi10 .and. zvolpn(i,j,n) > puny) THEN
1469	!OLI_CODE ! ice lid melted and category is pond covered
1470	!OLI_CODE zvolpn(i,j,n) = zvolpn(i,j,n) + zvuin(i,j,n)
1471	!OLI_CODE !fwoc(i,j) = fwoc(i,j) + zvuin(i,j,n)
1472	!OLI_CODE zvuin(i,j,n) = z0
1473	!OLI_CODE END IF
1474	!OLI_CODE zhpondn(i,j,n) = zvolpn(i,j,n) / zapondn(i,j,n)
1475	!OLI_CODE END IF
1476	!OLI_CODE
1477	!OLI_CODE !----------------------------------------------------------------
1478	!OLI_CODE ! freezing: existing upper ice layer grows
1479	!OLI_CODE !----------------------------------------------------------------
1480	!OLI_CODE ELSE IF (zvolpn(i,j,n) > epsi10) THEN ! zTavg <= zTp
1481	!OLI_CODE
1482	!OLI_CODE ! dIFferential growth of base of surface floating ice layer
1483	!OLI_CODE zdTice = max(-zTavg, z0) ! > 0
1484	!OLI_CODE zomega = rcdiczdTice zr1_rlfus
1485	!OLI_CODE zdHui = sqrt(zomegardt_ice + z25(zvuin(i,j,n)/ &
1486	!OLI_CODE aicen(i,j,n))*2)- z5 zvuin(i,j,n)/aicen(i,j,n)
1487	!OLI_CODE
1488	!OLI_CODE zdvice = min(zdHui*zapondn(i,j,n), zvolpn(i,j,n))
1489	!OLI_CODE IF (zdvice > epsi10) THEN
1490	!OLI_CODE zvuin (i,j,n) = zvuin (i,j,n) + zdvice
1491	!OLI_CODE zvolpn(i,j,n) = zvolpn(i,j,n) - zdvice
1492	!OLI_CODE zvolp (i,j) = zvolp (i,j) - zdvice
1493	!OLI_CODE !fwoc(i,j) = fwoc(i,j) - zdvice
1494	!OLI_CODE zhpondn(i,j,n) = zvolpn(i,j,n) / zapondn(i,j,n)
1495	!OLI_CODE END IF
1496	!OLI_CODE
1497	!OLI_CODE END IF ! zTavg
1498	!OLI_CODE
1499	!OLI_CODE !----------------------------------------------------------------
1500	!OLI_CODE ! freezing: upper ice layer begins to form
1501	!OLI_CODE ! note: albedo does not change
1502	!OLI_CODE !----------------------------------------------------------------
1503	!OLI_CODE ELSE ! zvuin < epsi10
1504	!OLI_CODE
1505	!OLI_CODE ! thickness of newly formed ice
1506	!OLI_CODE ! the surface temperature of a meltpond is the same as that
1507	!OLI_CODE ! of the ice underneath (0C), and the thermodynamic surface
1508	!OLI_CODE ! flux is the same
1509	!OLI_CODE zdHui = max(-fsurf(i,j)rdt_icezr1_rlfus, z0)
1510	!OLI_CODE zdvice = min(zdHui*zapondn(i,j,n), zvolpn(i,j,n))
1511	!OLI_CODE IF (zdvice > epsi10) THEN
1512	!OLI_CODE zvuin (i,j,n) = zdvice
1513	!OLI_CODE zvolpn(i,j,n) = zvolpn(i,j,n) - zdvice
1514	!OLI_CODE zvolp (i,j) = zvolp (i,j) - zdvice
1515	!OLI_CODE !fwoc(i,j) = fwoc(i,j) - zdvice
1516	!OLI_CODE zhpondn(i,j,n)= zvolpn(i,j,n) / zapondn(i,j,n)
1517	!OLI_CODE END IF
1518	!OLI_CODE
1519	!OLI_CODE END IF ! zvuin
1520	!OLI_CODE
1521	!OLI_CODE END DO ! jpl
1522	!OLI_CODE
1523	!OLI_CODE END DO ! ij
1524	!OLI_CODE
1525	!OLI_CODE !---------------------------------------------------------------
1526	!OLI_CODE ! remove ice lid if there is no liquid pond
1527	!OLI_CODE ! zvuin may be nonzero on category jpl due to dynamics
1528	!OLI_CODE !---------------------------------------------------------------
1529	!OLI_CODE
1530	!OLI_CODE DO j = 1, jpj
1531	!OLI_CODE DO i = 1, jpi
1532	!OLI_CODE DO n = 1, jpl
1533	!OLI_CODE IF (aicen(i,j,n) > epsi10 .and. zvolpn(i,j,n) < puny &
1534	!OLI_CODE .and. zvuin (i,j,n) > epsi10) THEN
1535	!OLI_CODE kcells(n) = kcells(n) + 1
1536	!OLI_CODE indxij = kcells(n)
1537	!OLI_CODE indxii(indxij,n) = i
1538	!OLI_CODE indxjj(indxij,n) = j
1539	!OLI_CODE END IF
1540	!OLI_CODE END DO
1541	!OLI_CODE END DO ! i
1542	!OLI_CODE END DO ! j
1543	!OLI_CODE
1544	!OLI_CODE DO n = 1, jpl
1545	!OLI_CODE
1546	!OLI_CODE IF (kcells(n) > 0) THEN
1547	!OLI_CODE DO ij = 1, kcells(n)
1548	!OLI_CODE i = indxii(ij,n)
1549	!OLI_CODE j = indxjj(ij,n)
1550	!OLI_CODE fwoc(i,j) = fwoc(i,j) + rhoic/rauw * zvuin(i,j,n) ! Completely refrozen lid goes into ocean (to be changed)
1551	!OLI_CODE zvuin(i,j,n) = z0
1552	!OLI_CODE END DO ! ij
1553	!OLI_CODE END IF
1554	!OLI_CODE
1555	!OLI_CODE ! reload tracers
1556	!OLI_CODE DO j = 1, jpj
1557	!OLI_CODE DO i = 1, jpi
1558	!OLI_CODE IF (zapondn(i,j,n) > epsi10) THEN
1559	!OLI_CODE h_il(i,j,n) = zvuin(i,j,n) / zapondn(i,j,n)
1560	!OLI_CODE ELSE
1561	!OLI_CODE zvuin(i,j,n) = z0
1562	!OLI_CODE h_il(i,j,n) = z0
1563	!OLI_CODE END IF
1564	!OLI_CODE IF (aicen(i,j,n) > epsi10) THEN
1565	!OLI_CODE a_ip_frac(i,j,n) = zapondn(i,j,n) / aicen(i,j,n) * &
1566	!OLI_CODE (1.0_wp - MAX(z0, SIGN(1.0_wp, -zvolpn(i,j,n))))
1567	!OLI_CODE h_ip(i,j,n) = zhpondn(i,j,n)
1568	!OLI_CODE ELSE
1569	!OLI_CODE a_ip_frac(i,j,n) = z0
1570	!OLI_CODE h_ip(i,j,n) = z0
1571	!OLI_CODE h_il(i,j,n) = z0
1572	!OLI_CODE END IF
1573	!OLI_CODE END DO ! i
1574	!OLI_CODE END DO ! j
1575	!OLI_CODE
1576	!OLI_CODE END DO ! n
1577	!OLI_CODE
1578	!OLI_CODE END SUBROUTINE compute_mp_topo
1579	!OLI_CODE
1580	!OLI_CODE
1581	!OLI_CODE SUBROUTINE pond_area(aice,vice,rhos, &
1582	!OLI_CODE aicen, vicen, vsnon, ticen, &
1583	!OLI_CODE salin, zvolpn, zvolp, &
1584	!OLI_CODE zapondn,zhpondn,dvolp)
1585	!OLI_CODE !!-------------------------------------------------------------------
1586	!OLI_CODE !! * ROUTINE pond_area *
1587	!OLI_CODE !!
1588	!OLI_CODE !! ** Purpose : Compute melt pond area, depth and melting rates
1589	!OLI_CODE !!------------------------------------------------------------------
1590	!OLI_CODE REAL (wp), INTENT(IN) :: &
1591	!OLI_CODE aice,vice
1592	!OLI_CODE
1593	!OLI_CODE REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
1594	!OLI_CODE aicen, vicen, vsnon, rhos
1595	!OLI_CODE
1596	!OLI_CODE REAL (wp), DIMENSION(nlay_i,jpl), INTENT(IN) :: &
1597	!OLI_CODE ticen, salin
1598	!OLI_CODE
1599	!OLI_CODE REAL (wp), DIMENSION(jpl), INTENT(INOUT) :: &
1600	!OLI_CODE zvolpn
1601	!OLI_CODE
1602	!OLI_CODE REAL (wp), INTENT(INOUT) :: &
1603	!OLI_CODE zvolp, dvolp
1604	!OLI_CODE
1605	!OLI_CODE REAL (wp), DIMENSION(jpl), INTENT(OUT) :: &
1606	!OLI_CODE zapondn, zhpondn
1607	!OLI_CODE
1608	!OLI_CODE INTEGER :: &
1609	!OLI_CODE n, ns, &
1610	!OLI_CODE m_index, &
1611	!OLI_CODE permflag
1612	!OLI_CODE
1613	!OLI_CODE REAL (wp), DIMENSION(jpl) :: &
1614	!OLI_CODE hicen, &
1615	!OLI_CODE hsnon, &
1616	!OLI_CODE asnon, &
1617	!OLI_CODE alfan, &
1618	!OLI_CODE betan, &
1619	!OLI_CODE cum_max_vol, &
1620	!OLI_CODE reduced_aicen
1621	!OLI_CODE
1622	!OLI_CODE REAL (wp), DIMENSION(0:jpl) :: &
1623	!OLI_CODE cum_max_vol_tmp
1624	!OLI_CODE
1625	!OLI_CODE REAL (wp) :: &
1626	!OLI_CODE hpond, &
1627	!OLI_CODE drain, &
1628	!OLI_CODE floe_weight, &
1629	!OLI_CODE pressure_head, &
1630	!OLI_CODE hsl_rel, &
1631	!OLI_CODE deltah, &
1632	!OLI_CODE perm, &
1633	!OLI_CODE apond, &
1634	!OLI_CODE msno
1635	!OLI_CODE
1636	!OLI_CODE REAL (wp), parameter :: &
1637	!OLI_CODE viscosity = 1.79e-3_wp, & ! kinematic water viscosity in kg/m/s
1638	!OLI_CODE z0 = 0.0_wp , &
1639	!OLI_CODE c1 = 1.0_wp , &
1640	!OLI_CODE p4 = 0.4_wp , &
1641	!OLI_CODE p6 = 0.6_wp , &
1642	!OLI_CODE epsi10 = 1.0e-11_wp
1643	!OLI_CODE
1644	!OLI_CODE !-----------\|
1645	!OLI_CODE ! \|
1646	!OLI_CODE ! \|-----------\|
1647	!OLI_CODE !___________\|___________\|______________________________________sea-level
1648	!OLI_CODE ! \| \|
1649	!OLI_CODE ! \| \|---^--------\|
1650	!OLI_CODE ! \| \| \| \|
1651	!OLI_CODE ! \| \| \| \|-----------\| \|-------
1652	!OLI_CODE ! \| \| \|alfan(n)\| \| \|
1653	!OLI_CODE ! \| \| \| \| \|--------------\|
1654	!OLI_CODE ! \| \| \| \| \| \|
1655	!OLI_CODE !---------------------------v-------------------------------------------
1656	!OLI_CODE ! \| \| ^ \| \| \|
1657	!OLI_CODE ! \| \| \| \| \|--------------\|
1658	!OLI_CODE ! \| \| \|betan(n)\| \| \|
1659	!OLI_CODE ! \| \| \| \|-----------\| \|-------
1660	!OLI_CODE ! \| \| \| \|
1661	!OLI_CODE ! \| \|---v------- \|
1662	!OLI_CODE ! \| \|
1663	!OLI_CODE ! \|-----------\|
1664	!OLI_CODE ! \|
1665	!OLI_CODE !-----------\|
1666	!OLI_CODE
1667	!OLI_CODE !-------------------------------------------------------------------
1668	!OLI_CODE ! initialize
1669	!OLI_CODE !-------------------------------------------------------------------
1670	!OLI_CODE
1671	!OLI_CODE DO n = 1, jpl
1672	!OLI_CODE
1673	!OLI_CODE zapondn(n) = z0
1674	!OLI_CODE zhpondn(n) = z0
1675	!OLI_CODE
1676	!OLI_CODE IF (aicen(n) < epsi10) THEN
1677	!OLI_CODE hicen(n) = z0
1678	!OLI_CODE hsnon(n) = z0
1679	!OLI_CODE reduced_aicen(n) = z0
1680	!OLI_CODE ELSE
1681	!OLI_CODE hicen(n) = vicen(n) / aicen(n)
1682	!OLI_CODE hsnon(n) = vsnon(n) / aicen(n)
1683	!OLI_CODE reduced_aicen(n) = c1 ! n=jpl
1684	!OLI_CODE IF (n < jpl) reduced_aicen(n) = aicen(n) &
1685	!OLI_CODE * (-0.024_wp*hicen(n) + 0.832_wp)
1686	!OLI_CODE asnon(n) = reduced_aicen(n)
1687	!OLI_CODE END IF
1688	!OLI_CODE
1689	!OLI_CODE ! This choice for alfa and beta ignores hydrostatic equilibium of categories.
1690	!OLI_CODE ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming
1691	!OLI_CODE ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all
1692	!OLI_CODE ! categories. alfa and beta partition the ITD - they are areas not thicknesses!
1693	!OLI_CODE ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area.
1694	!OLI_CODE ! Here, alfa = 60% of the ice area (and since hice is constant in a category,
1695	!OLI_CODE ! alfan = 60% of the ice volume) in each category lies above the reference line,
1696	!OLI_CODE ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required.
1697	!OLI_CODE
1698	!OLI_CODE alfan(n) = p6 * hicen(n)
1699	!OLI_CODE betan(n) = p4 * hicen(n)
1700	!OLI_CODE
1701	!OLI_CODE cum_max_vol(n) = z0
1702	!OLI_CODE cum_max_vol_tmp(n) = z0
1703	!OLI_CODE
1704	!OLI_CODE END DO ! jpl
1705	!OLI_CODE
1706	!OLI_CODE cum_max_vol_tmp(0) = z0
1707	!OLI_CODE drain = z0
1708	!OLI_CODE dvolp = z0
1709	!OLI_CODE
1710	!OLI_CODE !--------------------------------------------------------------------------
1711	!OLI_CODE ! the maximum amount of water that can be contained up to each ice category
1712	!OLI_CODE !--------------------------------------------------------------------------
1713	!OLI_CODE
1714	!OLI_CODE DO n = 1, jpl-1 ! last category can not hold any volume
1715	!OLI_CODE
1716	!OLI_CODE IF (alfan(n+1) >= alfan(n) .and. alfan(n+1) > z0) THEN
1717	!OLI_CODE
1718	!OLI_CODE ! total volume in level including snow
1719	!OLI_CODE cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) + &
1720	!OLI_CODE (alfan(n+1) - alfan(n)) * sum(reduced_aicen(1:n))
1721	!OLI_CODE
1722	!OLI_CODE
1723	!OLI_CODE ! subtract snow solid volumes from lower categories in current level
1724	!OLI_CODE DO ns = 1, n
1725	!OLI_CODE cum_max_vol_tmp(n) = cum_max_vol_tmp(n) &
1726	!OLI_CODE - rhos(ns)/rauw * & ! fraction of snow that is occupied by solid ??rauw
1727	!OLI_CODE asnon(ns) * & ! area of snow from that category
1728	!OLI_CODE max(min(hsnon(ns)+alfan(ns)-alfan(n), alfan(n+1)- &
1729	!OLI_CODE alfan(n)), z0)
1730	!OLI_CODE ! thickness of snow from ns layer in n layer
1731	!OLI_CODE END DO
1732	!OLI_CODE
1733	!OLI_CODE ELSE ! assume higher categories unoccupied
1734	!OLI_CODE cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1)
1735	!OLI_CODE END IF
1736	!OLI_CODE !IF (cum_max_vol_tmp(n) < z0) THEN
1737	!OLI_CODE ! call abort_ice('negative melt pond volume')
1738	!OLI_CODE !END IF
1739	!OLI_CODE END DO
1740	!OLI_CODE cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1) ! last category holds no volume
1741	!OLI_CODE cum_max_vol (1:jpl) = cum_max_vol_tmp(1:jpl)
1742	!OLI_CODE
1743	!OLI_CODE !----------------------------------------------------------------
1744	!OLI_CODE ! is there more meltwater than can be held in the floe?
1745	!OLI_CODE !----------------------------------------------------------------
1746	!OLI_CODE IF (zvolp >= cum_max_vol(jpl)) THEN
1747	!OLI_CODE drain = zvolp - cum_max_vol(jpl) + epsi10
1748	!OLI_CODE zvolp = zvolp - drain
1749	!OLI_CODE dvolp = drain
1750	!OLI_CODE IF (zvolp < epsi10) THEN
1751	!OLI_CODE dvolp = dvolp + zvolp
1752	!OLI_CODE zvolp = z0
1753	!OLI_CODE END IF
1754	!OLI_CODE END IF
1755	!OLI_CODE
1756	!OLI_CODE ! height and area corresponding to the remaining volume
1757	!OLI_CODE
1758	!OLI_CODE call calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, &
1759	!OLI_CODE zvolp, cum_max_vol, hpond, m_index)
1760	!OLI_CODE
1761	!OLI_CODE DO n=1, m_index
1762	!OLI_CODE zhpondn(n) = hpond - alfan(n) + alfan(1)
1763	!OLI_CODE zapondn(n) = reduced_aicen(n)
1764	!OLI_CODE END DO
1765	!OLI_CODE apond = sum(zapondn(1:m_index))
1766	!OLI_CODE
1767	!OLI_CODE !------------------------------------------------------------------------
1768	!OLI_CODE ! drainage due to ice permeability - Darcy's law
1769	!OLI_CODE !------------------------------------------------------------------------
1770	!OLI_CODE
1771	!OLI_CODE ! sea water level
1772	!OLI_CODE msno = z0
1773	!OLI_CODE DO n=1,jpl
1774	!OLI_CODE msno = msno + vsnon(n) * rhos(n)
1775	!OLI_CODE END DO
1776	!OLI_CODE floe_weight = (msno + rhoicvice + rau0zvolp) / aice
1777	!OLI_CODE hsl_rel = floe_weight / rau0 &
1778	!OLI_CODE - ((sum(betan(:)*aicen(:))/aice) + alfan(1))
1779	!OLI_CODE
1780	!OLI_CODE deltah = hpond - hsl_rel
1781	!OLI_CODE pressure_head = grav * rau0 * max(deltah, z0)
1782	!OLI_CODE
1783	!OLI_CODE ! drain IF ice is permeable
1784	!OLI_CODE permflag = 0
1785	!OLI_CODE IF (pressure_head > z0) THEN
1786	!OLI_CODE DO n = 1, jpl-1
1787	!OLI_CODE IF (hicen(n) /= z0) THEN
1788	!OLI_CODE CALL permeability_phi(ticen(:,n), salin(:,n), vicen(n), perm)
1789	!OLI_CODE IF (perm > z0) permflag = 1
1790	!OLI_CODE drain = permzapondn(n)pressure_head*rdt_ice / &
1791	!OLI_CODE (viscosity*hicen(n))
1792	!OLI_CODE dvolp = dvolp + min(drain, zvolp)
1793	!OLI_CODE zvolp = max(zvolp - drain, z0)
1794	!OLI_CODE IF (zvolp < epsi10) THEN
1795	!OLI_CODE dvolp = dvolp + zvolp
1796	!OLI_CODE zvolp = z0
1797	!OLI_CODE END IF
1798	!OLI_CODE END IF
1799	!OLI_CODE END DO
1800	!OLI_CODE
1801	!OLI_CODE ! adjust melt pond DIMENSIONs
1802	!OLI_CODE IF (permflag > 0) THEN
1803	!OLI_CODE ! recompute pond depth
1804	!OLI_CODE CALL calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, &
1805	!OLI_CODE zvolp, cum_max_vol, hpond, m_index)
1806	!OLI_CODE DO n=1, m_index
1807	!OLI_CODE zhpondn(n) = hpond - alfan(n) + alfan(1)
1808	!OLI_CODE zapondn(n) = reduced_aicen(n)
1809	!OLI_CODE END DO
1810	!OLI_CODE apond = sum(zapondn(1:m_index))
1811	!OLI_CODE END IF
1812	!OLI_CODE END IF ! pressure_head
1813	!OLI_CODE
1814	!OLI_CODE !------------------------------------------------------------------------
1815	!OLI_CODE ! total melt pond volume in category DOes not include snow volume
1816	!OLI_CODE ! snow in melt ponds is not melted
1817	!OLI_CODE !------------------------------------------------------------------------
1818	!OLI_CODE
1819	!OLI_CODE ! Calculate pond volume for lower categories
1820	!OLI_CODE DO n=1,m_index-1
1821	!OLI_CODE zvolpn(n) = zapondn(n) * zhpondn(n) &
1822	!OLI_CODE - (rhos(n)/rauw) * asnon(n) * min(hsnon(n), zhpondn(n))!
1823	!OLI_CODE END DO
1824	!OLI_CODE
1825	!OLI_CODE ! Calculate pond volume for highest category = remaining pond volume
1826	!OLI_CODE IF (m_index == 1) zvolpn(m_index) = zvolp
1827	!OLI_CODE IF (m_index > 1) THEN
1828	!OLI_CODE IF (zvolp > sum(zvolpn(1:m_index-1))) THEN
1829	!OLI_CODE zvolpn(m_index) = zvolp - sum(zvolpn(1:m_index-1))
1830	!OLI_CODE ELSE
1831	!OLI_CODE zvolpn(m_index) = z0
1832	!OLI_CODE zhpondn(m_index) = z0
1833	!OLI_CODE zapondn(m_index) = z0
1834	!OLI_CODE ! If remaining pond volume is negative reduce pond volume of
1835	!OLI_CODE ! lower category
1836	!OLI_CODE IF (zvolp+epsi10 < sum(zvolpn(1:m_index-1))) &
1837	!OLI_CODE zvolpn(m_index-1) = zvolpn(m_index-1)-sum(zvolpn(1:m_index-1))&
1838	!OLI_CODE + zvolp
1839	!OLI_CODE END IF
1840	!OLI_CODE END IF
1841	!OLI_CODE
1842	!OLI_CODE DO n=1,m_index
1843	!OLI_CODE IF (zapondn(n) > epsi10) THEN
1844	!OLI_CODE zhpondn(n) = zvolpn(n) / zapondn(n)
1845	!OLI_CODE ELSE
1846	!OLI_CODE dvolp = dvolp + zvolpn(n)
1847	!OLI_CODE zhpondn(n) = z0
1848	!OLI_CODE zvolpn(n) = z0
1849	!OLI_CODE zapondn(n) = z0
1850	!OLI_CODE end IF
1851	!OLI_CODE END DO
1852	!OLI_CODE DO n = m_index+1, jpl
1853	!OLI_CODE zhpondn(n) = z0
1854	!OLI_CODE zapondn(n) = z0
1855	!OLI_CODE zvolpn (n) = z0
1856	!OLI_CODE END DO
1857	!OLI_CODE
1858	!OLI_CODE END SUBROUTINE pond_area
1859	!OLI_CODE
1860	!OLI_CODE
1861	!OLI_CODE SUBROUTINE calc_hpond(aicen, asnon, hsnon, rhos, alfan, &
1862	!OLI_CODE zvolp, cum_max_vol, &
1863	!OLI_CODE hpond, m_index)
1864	!OLI_CODE !!-------------------------------------------------------------------
1865	!OLI_CODE !! * ROUTINE calc_hpond *
1866	!OLI_CODE !!
1867	!OLI_CODE !! ** Purpose : Compute melt pond depth
1868	!OLI_CODE !!-------------------------------------------------------------------
1869	!OLI_CODE
1870	!OLI_CODE REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
1871	!OLI_CODE aicen, &
1872	!OLI_CODE asnon, &
1873	!OLI_CODE hsnon, &
1874	!OLI_CODE rhos, &
1875	!OLI_CODE alfan, &
1876	!OLI_CODE cum_max_vol
1877	!OLI_CODE
1878	!OLI_CODE REAL (wp), INTENT(IN) :: &
1879	!OLI_CODE zvolp
1880	!OLI_CODE
1881	!OLI_CODE REAL (wp), INTENT(OUT) :: &
1882	!OLI_CODE hpond
1883	!OLI_CODE
1884	!OLI_CODE INTEGER, INTENT(OUT) :: &
1885	!OLI_CODE m_index
1886	!OLI_CODE
1887	!OLI_CODE INTEGER :: n, ns
1888	!OLI_CODE
1889	!OLI_CODE REAL (wp), DIMENSION(0:jpl+1) :: &
1890	!OLI_CODE hitl, &
1891	!OLI_CODE aicetl
1892	!OLI_CODE
1893	!OLI_CODE REAL (wp) :: &
1894	!OLI_CODE rem_vol, &
1895	!OLI_CODE area, &
1896	!OLI_CODE vol, &
1897	!OLI_CODE tmp, &
1898	!OLI_CODE z0 = 0.0_wp, &
1899	!OLI_CODE epsi10 = 1.0e-11_wp
1900	!OLI_CODE
1901	!OLI_CODE !----------------------------------------------------------------
1902	!OLI_CODE ! hpond is zero if zvolp is zero - have we fully drained?
1903	!OLI_CODE !----------------------------------------------------------------
1904	!OLI_CODE
1905	!OLI_CODE IF (zvolp < epsi10) THEN
1906	!OLI_CODE hpond = z0
1907	!OLI_CODE m_index = 0
1908	!OLI_CODE ELSE
1909	!OLI_CODE
1910	!OLI_CODE !----------------------------------------------------------------
1911	!OLI_CODE ! Calculate the category where water fills up to
1912	!OLI_CODE !----------------------------------------------------------------
1913	!OLI_CODE
1914	!OLI_CODE !----------\|
1915	!OLI_CODE ! \|
1916	!OLI_CODE ! \|
1917	!OLI_CODE ! \|----------\| -- --
1918	!OLI_CODE !__________\|__________\|_________________________________________ ^
1919	!OLI_CODE ! \| \| rem_vol ^ \| Semi-filled
1920	!OLI_CODE ! \| \|----------\|-- -- -- - ---\|-- ---- -- -- --v layer
1921	!OLI_CODE ! \| \| \| \|
1922	!OLI_CODE ! \| \| \| \|hpond
1923	!OLI_CODE ! \| \| \|----------\| \| \|-------
1924	!OLI_CODE ! \| \| \| \| \| \|
1925	!OLI_CODE ! \| \| \| \|---v-----\|
1926	!OLI_CODE ! \| \| m_index \| \| \|
1927	!OLI_CODE !-------------------------------------------------------------
1928	!OLI_CODE
1929	!OLI_CODE m_index = 0 ! 1:m_index categories have water in them
1930	!OLI_CODE DO n = 1, jpl
1931	!OLI_CODE IF (zvolp <= cum_max_vol(n)) THEN
1932	!OLI_CODE m_index = n
1933	!OLI_CODE IF (n == 1) THEN
1934	!OLI_CODE rem_vol = zvolp
1935	!OLI_CODE ELSE
1936	!OLI_CODE rem_vol = zvolp - cum_max_vol(n-1)
1937	!OLI_CODE END IF
1938	!OLI_CODE exit ! to break out of the loop
1939	!OLI_CODE END IF
1940	!OLI_CODE END DO
1941	!OLI_CODE m_index = min(jpl-1, m_index)
1942	!OLI_CODE
1943	!OLI_CODE !----------------------------------------------------------------
1944	!OLI_CODE ! semi-filled layer may have m_index different snow in it
1945	!OLI_CODE !----------------------------------------------------------------
1946	!OLI_CODE
1947	!OLI_CODE !----------------------------------------------------------- ^
1948	!OLI_CODE ! \| alfan(m_index+1)
1949	!OLI_CODE ! \|
1950	!OLI_CODE !hitl(3)--> \|----------\| \|
1951	!OLI_CODE !hitl(2)--> \|------------\| * * * * *\| \|
1952	!OLI_CODE !hitl(1)--> \|----------\|* * * * * * \|* * * * * \| \|
1953	!OLI_CODE !hitl(0)-->------------------------------------------------- \| ^
1954	!OLI_CODE ! various snow from lower categories \| \|alfa(m_index)
1955	!OLI_CODE
1956	!OLI_CODE ! hitl - heights of the snow layers from thinner and current categories
1957	!OLI_CODE ! aicetl - area of each snow depth in this layer
1958	!OLI_CODE
1959	!OLI_CODE hitl(:) = z0
1960	!OLI_CODE aicetl(:) = z0
1961	!OLI_CODE DO n = 1, m_index
1962	!OLI_CODE hitl(n) = max(min(hsnon(n) + alfan(n) - alfan(m_index), &
1963	!OLI_CODE alfan(m_index+1) - alfan(m_index)), z0)
1964	!OLI_CODE aicetl(n) = asnon(n)
1965	!OLI_CODE
1966	!OLI_CODE aicetl(0) = aicetl(0) + (aicen(n) - asnon(n))
1967	!OLI_CODE END DO
1968	!OLI_CODE hitl(m_index+1) = alfan(m_index+1) - alfan(m_index)
1969	!OLI_CODE aicetl(m_index+1) = z0
1970	!OLI_CODE
1971	!OLI_CODE !----------------------------------------------------------------
1972	!OLI_CODE ! reorder array according to hitl
1973	!OLI_CODE ! snow heights not necessarily in height order
1974	!OLI_CODE !----------------------------------------------------------------
1975	!OLI_CODE
1976	!OLI_CODE DO ns = 1, m_index+1
1977	!OLI_CODE DO n = 0, m_index - ns + 1
1978	!OLI_CODE IF (hitl(n) > hitl(n+1)) THEN ! swap order
1979	!OLI_CODE tmp = hitl(n)
1980	!OLI_CODE hitl(n) = hitl(n+1)
1981	!OLI_CODE hitl(n+1) = tmp
1982	!OLI_CODE tmp = aicetl(n)
1983	!OLI_CODE aicetl(n) = aicetl(n+1)
1984	!OLI_CODE aicetl(n+1) = tmp
1985	!OLI_CODE END IF
1986	!OLI_CODE END DO
1987	!OLI_CODE END DO
1988	!OLI_CODE
1989	!OLI_CODE !----------------------------------------------------------------
1990	!OLI_CODE ! divide semi-filled layer into set of sublayers each vertically homogenous
1991	!OLI_CODE !----------------------------------------------------------------
1992	!OLI_CODE
1993	!OLI_CODE !hitl(3)----------------------------------------------------------------
1994	!OLI_CODE ! \| * * * * * * * *
1995	!OLI_CODE ! \|* * * * * * * * *
1996	!OLI_CODE !hitl(2)----------------------------------------------------------------
1997	!OLI_CODE ! \| * * * * * * * * \| * * * * * * * *
1998	!OLI_CODE ! \|* * * * * * * * * \|* * * * * * * * *
1999	!OLI_CODE !hitl(1)----------------------------------------------------------------
2000	!OLI_CODE ! \| * * * * * * * * \| * * * * * * * * \| * * * * * * * *
2001	!OLI_CODE ! \|* * * * * * * * * \|* * * * * * * * * \|* * * * * * * * *
2002	!OLI_CODE !hitl(0)----------------------------------------------------------------
2003	!OLI_CODE ! aicetl(0) aicetl(1) aicetl(2) aicetl(3)
2004	!OLI_CODE
2005	!OLI_CODE ! move up over layers incrementing volume
2006	!OLI_CODE DO n = 1, m_index+1
2007	!OLI_CODE
2008	!OLI_CODE area = sum(aicetl(:)) - & ! total area of sub-layer
2009	!OLI_CODE (rhos(n)/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow
2010	!OLI_CODE
2011	!OLI_CODE vol = (hitl(n) - hitl(n-1)) * area ! thickness of sub-layer times area
2012	!OLI_CODE
2013	!OLI_CODE IF (vol >= rem_vol) THEN ! have reached the sub-layer with the depth within
2014	!OLI_CODE hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - &
2015	!OLI_CODE alfan(1)
2016	!OLI_CODE exit
2017	!OLI_CODE ELSE ! still in sub-layer below the sub-layer with the depth
2018	!OLI_CODE rem_vol = rem_vol - vol
2019	!OLI_CODE END IF
2020	!OLI_CODE
2021	!OLI_CODE END DO
2022	!OLI_CODE
2023	!OLI_CODE END IF
2024	!OLI_CODE
2025	!OLI_CODE END SUBROUTINE calc_hpond
2026	!OLI_CODE
2027	!OLI_CODE
2028	!OLI_CODE SUBROUTINE permeability_phi(ticen, salin, vicen, perm)
2029	!OLI_CODE !!-------------------------------------------------------------------
2030	!OLI_CODE !! * ROUTINE permeability_phi *
2031	!OLI_CODE !!
2032	!OLI_CODE !! ** Purpose : Determine the liquid fraction of brine in the ice
2033	!OLI_CODE !! and its permeability
2034	!OLI_CODE !!-------------------------------------------------------------------
2035	!OLI_CODE REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: &
2036	!OLI_CODE ticen, & ! energy of melting for each ice layer (J/m2)
2037	!OLI_CODE salin
2038	!OLI_CODE
2039	!OLI_CODE REAL (wp), INTENT(IN) :: &
2040	!OLI_CODE vicen ! ice volume
2041	!OLI_CODE
2042	!OLI_CODE REAL (wp), INTENT(OUT) :: &
2043	!OLI_CODE perm ! permeability
2044	!OLI_CODE
2045	!OLI_CODE REAL (wp) :: &
2046	!OLI_CODE Sbr ! brine salinity
2047	!OLI_CODE
2048	!OLI_CODE REAL (wp), DIMENSION(nlay_i) :: &
2049	!OLI_CODE Tin, & ! ice temperature
2050	!OLI_CODE phi ! liquid fraction
2051	!OLI_CODE
2052	!OLI_CODE INTEGER :: k
2053	!OLI_CODE
2054	!OLI_CODE REAL (wp) :: &
2055	!OLI_CODE c2 = 2.0_wp
2056	!OLI_CODE
2057	!OLI_CODE !-----------------------------------------------------------------
2058	!OLI_CODE ! Compute ice temperatures from enthalpies using quadratic formula
2059	!OLI_CODE !-----------------------------------------------------------------
2060	!OLI_CODE
2061	!OLI_CODE DO k = 1,nlay_i
2062	!OLI_CODE Tin(k) = ticen(k)
2063	!OLI_CODE END DO
2064	!OLI_CODE
2065	!OLI_CODE !-----------------------------------------------------------------
2066	!OLI_CODE ! brine salinity and liquid fraction
2067	!OLI_CODE !-----------------------------------------------------------------
2068	!OLI_CODE
2069	!OLI_CODE IF (maxval(Tin-rtt) <= -c2) THEN
2070	!OLI_CODE
2071	!OLI_CODE DO k = 1,nlay_i
2072	!OLI_CODE Sbr = - 1.2_wp &
2073	!OLI_CODE -21.8_wp * (Tin(k)-rtt) &
2074	!OLI_CODE - 0.919_wp * (Tin(k)-rtt)**2 &
2075	!OLI_CODE - 0.01878_wp * (Tin(k)-rtt)**3
2076	!OLI_CODE phi(k) = salin(k)/Sbr ! liquid fraction
2077	!OLI_CODE END DO ! k
2078	!OLI_CODE
2079	!OLI_CODE ELSE
2080	!OLI_CODE
2081	!OLI_CODE DO k = 1,nlay_i
2082	!OLI_CODE Sbr = -17.6_wp * (Tin(k)-rtt) &
2083	!OLI_CODE - 0.389_wp * (Tin(k)-rtt)**2 &
2084	!OLI_CODE - 0.00362_wp* (Tin(k)-rtt)**3
2085	!OLI_CODE phi(k) = salin(k)/Sbr ! liquid fraction
2086	!OLI_CODE END DO
2087	!OLI_CODE
2088	!OLI_CODE END IF
2089	!OLI_CODE
2090	!OLI_CODE !-----------------------------------------------------------------
2091	!OLI_CODE ! permeability
2092	!OLI_CODE !-----------------------------------------------------------------
2093	!OLI_CODE
2094	!OLI_CODE perm = 3.0e-08_wp * (minval(phi))**3
2095	!OLI_CODE
2096	!OLI_CODE END SUBROUTINE permeability_phi
2097	!OLI_CODE
2098	!OLI_CODE #else
2099	!OLI_CODE !!----------------------------------------------------------------------
2100	!OLI_CODE !! Default option Dummy Module No LIM-3 sea-ice model
2101	!OLI_CODE !!----------------------------------------------------------------------
2102	!OLI_CODE CONTAINS
2103	!OLI_CODE SUBROUTINE lim_mp_init ! Empty routine
2104	!OLI_CODE END SUBROUTINE lim_mp_init
2105	!OLI_CODE SUBROUTINE lim_mp ! Empty routine
2106	!OLI_CODE END SUBROUTINE lim_mp
2107	!OLI_CODE SUBROUTINE compute_mp_topo ! Empty routine
2108	!OLI_CODE END SUBROUTINE compute_mp_topo
2109	!OLI_CODE SUBROUTINE pond_area ! Empty routine
2110	!OLI_CODE END SUBROUTINE pond_area
2111	!OLI_CODE SUBROUTINE calc_hpond ! Empty routine
2112	!OLI_CODE END SUBROUTINE calc_hpond
2113	!OLI_CODE SUBROUTINE permeability_phy ! Empty routine
2114	!OLI_CODE END SUBROUTINE permeability_phy
2115	!OLI_CODE #endif
2116	!OLI_CODE !!======================================================================
2117	!OLI_CODE END MODULE limmp_topo
2118	!OLI_CODE
2119	#else
2120	!!----------------------------------------------------------------------
2121	!! Default option Empty module NO LIM sea-ice model
2122	!!----------------------------------------------------------------------
2123	CONTAINS
2124	SUBROUTINE lim_mp_init ! Empty routine
2125	END SUBROUTINE lim_mp_init
2126	SUBROUTINE lim_mp ! Empty routine
2127	END SUBROUTINE lim_mp
2128	SUBROUTINE lim_mp_topo ! Empty routine
2129	END SUBROUTINE lim_mp_topo
2130	SUBROUTINE lim_mp_cesm ! Empty routine
2131	END SUBROUTINE lim_mp_cesm
2132	SUBROUTINE lim_mp_area ! Empty routine
2133	END SUBROUTINE lim_mp_area
2134	SUBROUTINE lim_mp_perm ! Empty routine
2135	END SUBROUTINE lim_mp_perm
2136	#endif
2137
2138	!!======================================================================
2139	END MODULE limmp

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

Download in other formats: