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

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

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

Download in other formats: