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

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

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

Download in other formats: