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

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

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

Download in other formats: