New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

icethd_pnd.F90 in NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icethd_pnd.F90 @ 13850

Last change on this file since 13850 was 13850, checked in by vancop, 4 years ago
Running but slow topographic ponds
Property svn:keywords set to `Id`
File size: 59.2 KB

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

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

Download in other formats: