Context Navigation

icevar.F90

Last change on this file was 12636, checked in by dancopsey, 4 years ago

Add:

Melt pond lids
Melt pond maximum area and thickness
Melt pond vertical flushing
Area contributing to melt ponds depends on total ice fraction

File size: 57.9 KB

Line
1	MODULE icevar
2	!!======================================================================
3	!! * MODULE icevar *
4	!! sea-ice: series of functions to transform or compute ice variables
5	!!======================================================================
6	!! History : - ! 2006-01 (M. Vancoppenolle) Original code
7	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
8	!!----------------------------------------------------------------------
9	#if defined key_si3
10	!!----------------------------------------------------------------------
11	!! 'key_si3' SI3 sea-ice model
12	!!----------------------------------------------------------------------
13	!!
14	!! There are three sets of variables
15	!! VGLO : global variables of the model
16	!! - v_i (jpi,jpj,jpl)
17	!! - v_s (jpi,jpj,jpl)
18	!! - a_i (jpi,jpj,jpl)
19	!! - t_s (jpi,jpj,jpl)
20	!! - e_i (jpi,jpj,nlay_i,jpl)
21	!! - e_s (jpi,jpj,nlay_s,jpl)
22	!! - sv_i(jpi,jpj,jpl)
23	!! - oa_i(jpi,jpj,jpl)
24	!! VEQV : equivalent variables sometimes used in the model
25	!! - h_i(jpi,jpj,jpl)
26	!! - h_s(jpi,jpj,jpl)
27	!! - t_i(jpi,jpj,nlay_i,jpl)
28	!! ...
29	!! VAGG : aggregate variables, averaged/summed over all
30	!! thickness categories
31	!! - vt_i(jpi,jpj)
32	!! - vt_s(jpi,jpj)
33	!! - at_i(jpi,jpj)
34	!! - st_i(jpi,jpj)
35	!! - et_s(jpi,jpj) total snow heat content
36	!! - et_i(jpi,jpj) total ice thermal content
37	!! - sm_i(jpi,jpj) mean ice salinity
38	!! - tm_i(jpi,jpj) mean ice temperature
39	!! - tm_s(jpi,jpj) mean snw temperature
40	!!----------------------------------------------------------------------
41	!! ice_var_agg : integrate variables over layers and categories
42	!! ice_var_glo2eqv : transform from VGLO to VEQV
43	!! ice_var_eqv2glo : transform from VEQV to VGLO
44	!! ice_var_salprof : salinity profile in the ice
45	!! ice_var_salprof1d : salinity profile in the ice 1D
46	!! ice_var_zapsmall : remove very small area and volume
47	!! ice_var_zapneg : remove negative ice fields
48	!! ice_var_roundoff : remove negative values arising from roundoff erros
49	!! ice_var_bv : brine volume
50	!! ice_var_enthalpy : compute ice and snow enthalpies from temperature
51	!! ice_var_sshdyn : compute equivalent ssh in lead
52	!! ice_var_itd : convert N-cat to M-cat
53	!!----------------------------------------------------------------------
54	USE dom_oce ! ocean space and time domain
55	USE phycst ! physical constants (ocean directory)
56	USE sbc_oce , ONLY : sss_m, ln_ice_embd, nn_fsbc
57	USE ice ! sea-ice: variables
58	USE ice1D ! sea-ice: thermodynamics variables
59	!
60	USE in_out_manager ! I/O manager
61	USE lib_mpp ! MPP library
62	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
63
64	IMPLICIT NONE
65	PRIVATE
66
67	PUBLIC ice_var_agg
68	PUBLIC ice_var_glo2eqv
69	PUBLIC ice_var_eqv2glo
70	PUBLIC ice_var_salprof
71	PUBLIC ice_var_salprof1d
72	PUBLIC ice_var_zapsmall
73	PUBLIC ice_var_zapneg
74	PUBLIC ice_var_roundoff
75	PUBLIC ice_var_bv
76	PUBLIC ice_var_enthalpy
77	PUBLIC ice_var_sshdyn
78	PUBLIC ice_var_itd
79
80	INTERFACE ice_var_itd
81	MODULE PROCEDURE ice_var_itd_1c1c, ice_var_itd_Nc1c, ice_var_itd_1cMc, ice_var_itd_NcMc
82	END INTERFACE
83
84	!!----------------------------------------------------------------------
85	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
86	!! $Id$
87	!! Software governed by the CeCILL license (see ./LICENSE)
88	!!----------------------------------------------------------------------
89	CONTAINS
90
91	SUBROUTINE ice_var_agg( kn )
92	!!-------------------------------------------------------------------
93	!! * ROUTINE ice_var_agg *
94	!!
95	!! ** Purpose : aggregates ice-thickness-category variables to
96	!! all-ice variables, i.e. it turns VGLO into VAGG
97	!!-------------------------------------------------------------------
98	INTEGER, INTENT( in ) :: kn ! =1 state variables only
99	! ! >1 state variables + others
100	!
101	INTEGER :: ji, jj, jk, jl ! dummy loop indices
102	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z1_at_i, z1_vt_i, z1_vt_s
103	!!-------------------------------------------------------------------
104	!
105	! ! integrated values
106	vt_i(:,:) = SUM( v_i (:,:,:) , dim=3 )
107	vt_s(:,:) = SUM( v_s (:,:,:) , dim=3 )
108	st_i(:,:) = SUM( sv_i(:,:,:) , dim=3 )
109	at_i(:,:) = SUM( a_i (:,:,:) , dim=3 )
110	et_s(:,:) = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 )
111	et_i(:,:) = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 )
112	!
113	at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds
114	vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
115	!
116	ato_i(:,:) = 1._wp - at_i(:,:) ! open water fraction
117
118	! The following fields are calculated for diagnostics and outputs only
119	! ==> Do not use them for other purposes
120	IF( kn > 1 ) THEN
121	!
122	ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) )
123	WHERE( at_i(:,:) > epsi20 ) ; z1_at_i(:,:) = 1._wp / at_i(:,:)
124	ELSEWHERE ; z1_at_i(:,:) = 0._wp
125	END WHERE
126	WHERE( vt_i(:,:) > epsi20 ) ; z1_vt_i(:,:) = 1._wp / vt_i(:,:)
127	ELSEWHERE ; z1_vt_i(:,:) = 0._wp
128	END WHERE
129	WHERE( vt_s(:,:) > epsi20 ) ; z1_vt_s(:,:) = 1._wp / vt_s(:,:)
130	ELSEWHERE ; z1_vt_s(:,:) = 0._wp
131	END WHERE
132	!
133	! ! mean ice/snow thickness
134	hm_i(:,:) = vt_i(:,:) * z1_at_i(:,:)
135	hm_s(:,:) = vt_s(:,:) * z1_at_i(:,:)
136	!
137	! ! mean temperature (K), salinity and age
138	tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
139	tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
140	om_i (:,:) = SUM( oa_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
141	sm_i (:,:) = st_i(:,:) * z1_vt_i(:,:)
142	!
143	tm_i(:,:) = 0._wp
144	tm_s(:,:) = 0._wp
145	DO jl = 1, jpl
146	DO jk = 1, nlay_i
147	tm_i(:,:) = tm_i(:,:) + r1_nlay_i * t_i (:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
148	END DO
149	DO jk = 1, nlay_s
150	tm_s(:,:) = tm_s(:,:) + r1_nlay_s * t_s (:,:,jk,jl) * v_s(:,:,jl) * z1_vt_s(:,:)
151	END DO
152	END DO
153	!
154	! ! put rt0 where there is no ice
155	WHERE( at_i(:,:)<=epsi20 )
156	tm_su(:,:) = rt0
157	tm_si(:,:) = rt0
158	tm_i (:,:) = rt0
159	tm_s (:,:) = rt0
160	END WHERE
161	!
162	! ! mean melt pond depth
163	WHERE( at_ip(:,:) > epsi20 ) ; hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:)
164	ELSEWHERE ; hm_ip(:,:) = 0._wp
165	END WHERE
166	!
167	DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s )
168	!
169	ENDIF
170	!
171	END SUBROUTINE ice_var_agg
172
173
174	SUBROUTINE ice_var_glo2eqv
175	!!-------------------------------------------------------------------
176	!! * ROUTINE ice_var_glo2eqv *
177	!!
178	!! ** Purpose : computes equivalent variables as function of
179	!! global variables, i.e. it turns VGLO into VEQV
180	!!-------------------------------------------------------------------
181	INTEGER :: ji, jj, jk, jl ! dummy loop indices
182	REAL(wp) :: ze_i ! local scalars
183	REAL(wp) :: ze_s, ztmelts, zbbb, zccc ! - -
184	REAL(wp) :: zhmax, z1_zhmax ! - -
185	REAL(wp) :: zlay_i, zlay_s ! - -
186	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i, z1_v_i
187	!!-------------------------------------------------------------------
188
189	!!gm Question 2: It is possible to define existence of sea-ice in a common way between
190	!! ice area and ice volume ?
191	!! the idea is to be able to define one for all at the begining of this routine
192	!! a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
193
194	!---------------------------------------------------------------
195	! Ice thickness, snow thickness, ice salinity, ice age and ponds
196	!---------------------------------------------------------------
197	! !--- inverse of the ice area
198	WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
199	ELSEWHERE ; z1_a_i(:,:,:) = 0._wp
200	END WHERE
201	!
202	WHERE( v_i(:,:,:) > epsi20 ) ; z1_v_i(:,:,:) = 1._wp / v_i(:,:,:)
203	ELSEWHERE ; z1_v_i(:,:,:) = 0._wp
204	END WHERE
205	! !--- ice thickness
206	h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
207
208	zhmax = hi_max(jpl)
209	z1_zhmax = 1._wp / hi_max(jpl)
210	WHERE( h_i(:,:,jpl) > zhmax ) ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area
211	h_i (:,:,jpl) = zhmax
212	a_i (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax
213	z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl)
214	END WHERE
215	! !--- snow thickness
216	h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:)
217	! !--- ice age
218	o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:)
219	! !--- pond fraction and thickness
220	a_ip_frac(:,:,:) = a_ip(:,:,:) * z1_a_i(:,:,:)
221	WHERE( a_ip_frac(:,:,:) > epsi20 ) ; h_ip(:,:,:) = v_ip(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:)
222	ELSEWHERE ; h_ip(:,:,:) = 0._wp
223	END WHERE
224	!
225	! !--- salinity (with a minimum value imposed everywhere)
226	IF( nn_icesal == 2 ) THEN
227	WHERE( v_i(:,:,:) > epsi20 ) ; s_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, sv_i(:,:,:) * z1_v_i(:,:,:) ) )
228	ELSEWHERE ; s_i(:,:,:) = rn_simin
229	END WHERE
230	ENDIF
231	CALL ice_var_salprof ! salinity profile
232
233	!-------------------
234	! Ice temperature [K] (with a minimum value (rt0 - 100.))
235	!-------------------
236	zlay_i = REAL( nlay_i , wp ) ! number of layers
237	DO jl = 1, jpl
238	DO jk = 1, nlay_i
239	DO jj = 1, jpj
240	DO ji = 1, jpi
241	IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area
242	!
243	ze_i = e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3]
244	ztmelts = - sz_i(ji,jj,jk,jl) * rTmlt ! Ice layer melt temperature [C]
245	! Conversion q(S,T) -> T (second order equation)
246	zbbb = ( rcp - rcpi ) * ztmelts + ze_i * r1_rhoi - rLfus
247	zccc = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts , 0._wp) )
248	t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts
249	!
250	ELSE !--- no ice
251	t_i(ji,jj,jk,jl) = rt0
252	ENDIF
253	END DO
254	END DO
255	END DO
256	END DO
257
258	!--------------------
259	! Snow temperature [K] (with a minimum value (rt0 - 100.))
260	!--------------------
261	zlay_s = REAL( nlay_s , wp )
262	DO jk = 1, nlay_s
263	WHERE( v_s(:,:,:) > epsi20 ) !--- icy area
264	t_s(:,:,jk,:) = rt0 + MAX( -100._wp , &
265	& MIN( r1_rcpi * ( -r1_rhos * ( e_s(:,:,jk,:) / v_s(:,:,:) * zlay_s ) + rLfus ) , 0._wp ) )
266	ELSEWHERE !--- no ice
267	t_s(:,:,jk,:) = rt0
268	END WHERE
269	END DO
270	!
271	! integrated values
272	vt_i (:,:) = SUM( v_i , dim=3 )
273	vt_s (:,:) = SUM( v_s , dim=3 )
274	at_i (:,:) = SUM( a_i , dim=3 )
275	!
276	END SUBROUTINE ice_var_glo2eqv
277
278
279	SUBROUTINE ice_var_eqv2glo
280	!!-------------------------------------------------------------------
281	!! * ROUTINE ice_var_eqv2glo *
282	!!
283	!! ** Purpose : computes global variables as function of
284	!! equivalent variables, i.e. it turns VEQV into VGLO
285	!!-------------------------------------------------------------------
286	!
287	v_i (:,:,:) = h_i (:,:,:) * a_i (:,:,:)
288	v_s (:,:,:) = h_s (:,:,:) * a_i (:,:,:)
289	sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:)
290	v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
291	!
292	END SUBROUTINE ice_var_eqv2glo
293
294
295	SUBROUTINE ice_var_salprof
296	!!-------------------------------------------------------------------
297	!! * ROUTINE ice_var_salprof *
298	!!
299	!! ** Purpose : computes salinity profile in function of bulk salinity
300	!!
301	!! ** Method : If bulk salinity greater than zsi1,
302	!! the profile is assumed to be constant (S_inf)
303	!! If bulk salinity lower than zsi0,
304	!! the profile is linear with 0 at the surface (S_zero)
305	!! If it is between zsi0 and zsi1, it is a
306	!! alpha-weighted linear combination of s_inf and s_zero
307	!!
308	!! ** References : Vancoppenolle et al., 2007
309	!!-------------------------------------------------------------------
310	INTEGER :: ji, jj, jk, jl ! dummy loop index
311	REAL(wp) :: zsal, z1_dS
312	REAL(wp) :: zargtemp , zs0, zs
313	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_slope_s, zalpha ! case 2 only
314	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
315	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
316	!!-------------------------------------------------------------------
317
318	!!gm Question: Remove the option 3 ? How many years since it last use ?
319
320	SELECT CASE ( nn_icesal )
321	!
322	! !---------------------------------------!
323	CASE( 1 ) ! constant salinity in time and space !
324	! !---------------------------------------!
325	sz_i(:,:,:,:) = rn_icesal
326	s_i (:,:,:) = rn_icesal
327	!
328	! !---------------------------------------------!
329	CASE( 2 ) ! time varying salinity with linear profile !
330	! !---------------------------------------------!
331	!
332	ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
333	!
334	DO jl = 1, jpl
335	DO jk = 1, nlay_i
336	sz_i(:,:,jk,jl) = s_i(:,:,jl)
337	END DO
338	END DO
339	! ! Slope of the linear profile
340	WHERE( h_i(:,:,:) > epsi20 ) ; z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:)
341	ELSEWHERE ; z_slope_s(:,:,:) = 0._wp
342	END WHERE
343	!
344	z1_dS = 1._wp / ( zsi1 - zsi0 )
345	DO jl = 1, jpl
346	DO jj = 1, jpj
347	DO ji = 1, jpi
348	zalpha(ji,jj,jl) = MAX( 0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp ) )
349	! ! force a constant profile when SSS too low (Baltic Sea)
350	IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) ) zalpha(ji,jj,jl) = 0._wp
351	END DO
352	END DO
353	END DO
354	!
355	! Computation of the profile
356	DO jl = 1, jpl
357	DO jk = 1, nlay_i
358	DO jj = 1, jpj
359	DO ji = 1, jpi
360	! ! linear profile with 0 surface value
361	zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
362	zs = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl) ! weighting the profile
363	sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
364	END DO
365	END DO
366	END DO
367	END DO
368	!
369	DEALLOCATE( z_slope_s , zalpha )
370	!
371	! !-------------------------------------------!
372	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
373	! !-------------------------------------------! (mean = 2.30)
374	!
375	s_i(:,:,:) = 2.30_wp
376	!!gm Remark: if we keep the case 3, then compute an store one for all time-step
377	!! a array S_prof(1:nlay_i) containing the calculation and just do:
378	! DO jk = 1, nlay_i
379	! sz_i(:,:,jk,:) = S_prof(jk)
380	! END DO
381	!!gm end
382	!
383	DO jl = 1, jpl
384	DO jk = 1, nlay_i
385	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
386	sz_i(:,:,jk,jl) = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) )
387	END DO
388	END DO
389	!
390	END SELECT
391	!
392	END SUBROUTINE ice_var_salprof
393
394
395	SUBROUTINE ice_var_salprof1d
396	!!-------------------------------------------------------------------
397	!! * ROUTINE ice_var_salprof1d *
398	!!
399	!! ** Purpose : 1d computation of the sea ice salinity profile
400	!! Works with 1d vectors and is used by thermodynamic modules
401	!!-------------------------------------------------------------------
402	INTEGER :: ji, jk ! dummy loop indices
403	REAL(wp) :: zargtemp, zsal, z1_dS ! local scalars
404	REAL(wp) :: zs, zs0 ! - -
405	!
406	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_slope_s, zalpha !
407	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
408	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
409	!!-------------------------------------------------------------------
410	!
411	SELECT CASE ( nn_icesal )
412	!
413	! !---------------------------------------!
414	CASE( 1 ) ! constant salinity in time and space !
415	! !---------------------------------------!
416	sz_i_1d(1:npti,:) = rn_icesal
417	!
418	! !---------------------------------------------!
419	CASE( 2 ) ! time varying salinity with linear profile !
420	! !---------------------------------------------!
421	!
422	ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
423	!
424	! ! Slope of the linear profile
425	WHERE( h_i_1d(1:npti) > epsi20 ) ; z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
426	ELSEWHERE ; z_slope_s(1:npti) = 0._wp
427	END WHERE
428
429	z1_dS = 1._wp / ( zsi1 - zsi0 )
430	DO ji = 1, npti
431	zalpha(ji) = MAX( 0._wp , MIN( ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp ) )
432	! ! force a constant profile when SSS too low (Baltic Sea)
433	IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) ) zalpha(ji) = 0._wp
434	END DO
435	!
436	! Computation of the profile
437	DO jk = 1, nlay_i
438	DO ji = 1, npti
439	! ! linear profile with 0 surface value
440	zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
441	zs = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * s_i_1d(ji)
442	sz_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
443	END DO
444	END DO
445	!
446	DEALLOCATE( z_slope_s, zalpha )
447
448	! !-------------------------------------------!
449	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
450	! !-------------------------------------------! (mean = 2.30)
451	!
452	s_i_1d(1:npti) = 2.30_wp
453	!
454	!!gm cf remark in ice_var_salprof routine, CASE( 3 )
455	DO jk = 1, nlay_i
456	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
457	zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
458	DO ji = 1, npti
459	sz_i_1d(ji,jk) = zsal
460	END DO
461	END DO
462	!
463	END SELECT
464	!
465	END SUBROUTINE ice_var_salprof1d
466
467
468	SUBROUTINE ice_var_zapsmall
469	!!-------------------------------------------------------------------
470	!! * ROUTINE ice_var_zapsmall *
471	!!
472	!! ** Purpose : Remove too small sea ice areas and correct fluxes
473	!!-------------------------------------------------------------------
474	INTEGER :: ji, jj, jl, jk ! dummy loop indices
475	REAL(wp), DIMENSION(jpi,jpj) :: zswitch
476	!!-------------------------------------------------------------------
477	!
478	DO jl = 1, jpl !== loop over the categories ==!
479	!
480	WHERE( a_i(:,:,jl) > epsi10 ) ; h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
481	ELSEWHERE ; h_i(:,:,jl) = 0._wp
482	END WHERE
483	!
484	WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 ) ; zswitch(:,:) = 0._wp
485	ELSEWHERE ; zswitch(:,:) = 1._wp
486	END WHERE
487	!
488	!-----------------------------------------------------------------
489	! Zap ice energy and use ocean heat to melt ice
490	!-----------------------------------------------------------------
491	DO jk = 1, nlay_i
492	DO jj = 1 , jpj
493	DO ji = 1 , jpi
494	! update exchanges with ocean
495	hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
496	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
497	t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
498	END DO
499	END DO
500	END DO
501	!
502	DO jk = 1, nlay_s
503	DO jj = 1 , jpj
504	DO ji = 1 , jpi
505	! update exchanges with ocean
506	hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
507	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj)
508	t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
509	END DO
510	END DO
511	END DO
512	!
513	!-----------------------------------------------------------------
514	! zap ice and snow volume, add water and salt to ocean
515	!-----------------------------------------------------------------
516	DO jj = 1 , jpj
517	DO ji = 1 , jpi
518	! update exchanges with ocean
519	sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoi * r1_rdtice
520	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoi * r1_rdtice
521	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhos * r1_rdtice
522	!
523	a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
524	v_i (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
525	v_s (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
526	t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
527	oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
528	sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
529	!
530	h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
531	h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
532	!
533	a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
534	v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
535	lh_ip (ji,jj,jl) = lh_ip (ji,jj,jl) * zswitch(ji,jj)
536	!
537	END DO
538	END DO
539	!
540	END DO
541
542	! to be sure that at_i is the sum of a_i(jl)
543	at_i (:,:) = SUM( a_i (:,:,:), dim=3 )
544	vt_i (:,:) = SUM( v_i (:,:,:), dim=3 )
545	!!clem add?
546	! vt_s (:,:) = SUM( v_s (:,:,:), dim=3 )
547	! st_i (:,:) = SUM( sv_i(:,:,:), dim=3 )
548	! et_s(:,:) = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 )
549	! et_i(:,:) = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 )
550	!!clem
551
552	! open water = 1 if at_i=0
553	WHERE( at_i(:,:) == 0._wp ) ato_i(:,:) = 1._wp
554	!
555	END SUBROUTINE ice_var_zapsmall
556
557
558	SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i )
559	!!-------------------------------------------------------------------
560	!! * ROUTINE ice_var_zapneg *
561	!!
562	!! ** Purpose : Remove negative sea ice fields and correct fluxes
563	!!-------------------------------------------------------------------
564	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
565	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area
566	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume
567	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_s ! snw volume
568	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psv_i ! salt content
569	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content
570	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration
571	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
572	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume
573	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: plh_ip ! melt pond lid thickness
574	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content
575	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content
576	!
577	INTEGER :: ji, jj, jl, jk ! dummy loop indices
578	REAL(wp) :: z1_dt
579	!!-------------------------------------------------------------------
580	!
581	z1_dt = 1._wp / pdt
582	!
583	DO jl = 1, jpl !== loop over the categories ==!
584	!
585	! make sure a_i=0 where v_i<=0
586	WHERE( pv_i(:,:,:) <= 0._wp ) pa_i(:,:,:) = 0._wp
587
588	!----------------------------------------
589	! zap ice energy and send it to the ocean
590	!----------------------------------------
591	DO jk = 1, nlay_i
592	DO jj = 1 , jpj
593	DO ji = 1 , jpi
594	IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
595	hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
596	pe_i(ji,jj,jk,jl) = 0._wp
597	ENDIF
598	END DO
599	END DO
600	END DO
601	!
602	DO jk = 1, nlay_s
603	DO jj = 1 , jpj
604	DO ji = 1 , jpi
605	IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
606	hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
607	pe_s(ji,jj,jk,jl) = 0._wp
608	ENDIF
609	END DO
610	END DO
611	END DO
612	!
613	!-----------------------------------------------------
614	! zap ice and snow volume, add water and salt to ocean
615	!-----------------------------------------------------
616	DO jj = 1 , jpj
617	DO ji = 1 , jpi
618	IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
619	wfx_res(ji,jj) = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
620	pv_i (ji,jj,jl) = 0._wp
621	ENDIF
622	IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
623	wfx_res(ji,jj) = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt
624	pv_s (ji,jj,jl) = 0._wp
625	ENDIF
626	IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
627	sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt
628	psv_i (ji,jj,jl) = 0._wp
629	ENDIF
630	END DO
631	END DO
632	!
633	END DO
634	!
635	WHERE( pato_i(:,:) < 0._wp ) pato_i(:,:) = 0._wp
636	WHERE( poa_i (:,:,:) < 0._wp ) poa_i (:,:,:) = 0._wp
637	WHERE( pa_i (:,:,:) < 0._wp ) pa_i (:,:,:) = 0._wp
638	WHERE( pa_ip (:,:,:) < 0._wp ) pa_ip (:,:,:) = 0._wp
639	WHERE( pv_ip (:,:,:) < 0._wp ) pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
640	! but it does not change conservation, so keep it this way is ok
641	WHERE( plh_ip (:,:,:) < 0._wp ) plh_ip (:,:,:) = 0._wp
642	!
643	END SUBROUTINE ice_var_zapneg
644
645
646	SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i )
647	!!-------------------------------------------------------------------
648	!! * ROUTINE ice_var_roundoff *
649	!!
650	!! ** Purpose : Remove negative sea ice values arising from roundoff errors
651	!!-------------------------------------------------------------------
652	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_i ! ice concentration
653	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_i ! ice volume
654	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_s ! snw volume
655	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: psv_i ! salt content
656	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: poa_i ! age content
657	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
658	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_ip ! melt pond volume
659	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: plh_ip ! melt pond lid thickness
660	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_s ! snw heat content
661	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_i ! ice heat content
662	!!-------------------------------------------------------------------
663	!
664	WHERE( pa_i (1:npti,:) < 0._wp .AND. pa_i (1:npti,:) > -epsi10 ) pa_i (1:npti,:) = 0._wp ! a_i must be >= 0
665	WHERE( pv_i (1:npti,:) < 0._wp .AND. pv_i (1:npti,:) > -epsi10 ) pv_i (1:npti,:) = 0._wp ! v_i must be >= 0
666	WHERE( pv_s (1:npti,:) < 0._wp .AND. pv_s (1:npti,:) > -epsi10 ) pv_s (1:npti,:) = 0._wp ! v_s must be >= 0
667	WHERE( psv_i(1:npti,:) < 0._wp .AND. psv_i(1:npti,:) > -epsi10 ) psv_i(1:npti,:) = 0._wp ! sv_i must be >= 0
668	WHERE( poa_i(1:npti,:) < 0._wp .AND. poa_i(1:npti,:) > -epsi10 ) poa_i(1:npti,:) = 0._wp ! oa_i must be >= 0
669	WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 ) pe_i (1:npti,:,:) = 0._wp ! e_i must be >= 0
670	WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 ) pe_s (1:npti,:,:) = 0._wp ! e_s must be >= 0
671	IF( ln_pnd_H12 ) THEN
672	WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 ) pa_ip(1:npti,:) = 0._wp ! a_ip must be >= 0
673	WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 ) pv_ip(1:npti,:) = 0._wp ! v_ip must be >= 0
674	WHERE( plh_ip(1:npti,:) < 0._wp .AND. plh_ip(1:npti,:) > -epsi10 ) plh_ip(1:npti,:) = 0._wp ! lh_ip must be >= 0
675	ENDIF
676	!
677	END SUBROUTINE ice_var_roundoff
678
679
680	SUBROUTINE ice_var_bv
681	!!-------------------------------------------------------------------
682	!! * ROUTINE ice_var_bv *
683	!!
684	!! ** Purpose : computes mean brine volume (%) in sea ice
685	!!
686	!! ** Method : e = - 0.054 * S (ppt) / T (C)
687	!!
688	!! References : Vancoppenolle et al., JGR, 2007
689	!!-------------------------------------------------------------------
690	INTEGER :: ji, jj, jk, jl ! dummy loop indices
691	!!-------------------------------------------------------------------
692	!
693	!!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done
694	!! instead of setting everything to zero as just below
695	bv_i (:,:,:) = 0._wp
696	DO jl = 1, jpl
697	DO jk = 1, nlay_i
698	WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )
699	bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
700	END WHERE
701	END DO
702	END DO
703	WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
704	ELSEWHERE ; bvm_i(:,:) = 0._wp
705	END WHERE
706	!
707	END SUBROUTINE ice_var_bv
708
709
710	SUBROUTINE ice_var_enthalpy
711	!!-------------------------------------------------------------------
712	!! * ROUTINE ice_var_enthalpy *
713	!!
714	!! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature
715	!!
716	!! ** Method : Formula (Bitz and Lipscomb, 1999)
717	!!-------------------------------------------------------------------
718	INTEGER :: ji, jk ! dummy loop indices
719	REAL(wp) :: ztmelts ! local scalar
720	!!-------------------------------------------------------------------
721	!
722	DO jk = 1, nlay_i ! Sea ice energy of melting
723	DO ji = 1, npti
724	ztmelts = - rTmlt * sz_i_1d(ji,jk)
725	t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point => likely conservation issue
726	! (sometimes zdf scheme produces abnormally high temperatures)
727	e_i_1d(ji,jk) = rhoi * ( rcpi * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) ) &
728	& + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) ) &
729	& - rcp * ztmelts )
730	END DO
731	END DO
732	DO jk = 1, nlay_s ! Snow energy of melting
733	DO ji = 1, npti
734	e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
735	END DO
736	END DO
737	!
738	END SUBROUTINE ice_var_enthalpy
739
740
741	FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
742	!!---------------------------------------------------------------------
743	!! * ROUTINE ice_var_sshdyn *
744	!!
745	!! ** Purpose : compute the equivalent ssh in lead when sea ice is embedded
746	!!
747	!! ** Method : ssh_lead = ssh + (Mice + Msnow) / rau0
748	!!
749	!! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira,
750	!! Sea ice-ocean coupling using a rescaled vertical coordinate z*,
751	!! Ocean Modelling, Volume 24, Issues 1-2, 2008
752	!!----------------------------------------------------------------------
753	!
754	! input
755	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh !: ssh [m]
756	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass !: mass of snow and ice at current ice time step [Kg/m2]
757	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b !: mass of snow and ice at previous ice time step [Kg/m2]
758	!
759	! output
760	REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn ! equivalent ssh in lead [m]
761	!
762	! temporary
763	REAL(wp) :: zintn, zintb ! time interpolation weights []
764	REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload ! snow and ice load [m]
765	!
766	! compute ice load used to define the equivalent ssh in lead
767	IF( ln_ice_embd ) THEN
768	!
769	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
770	! = (1/nn_fsbc)^2 * {SUM[n] , n=0,nn_fsbc-1}
771	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
772	!
773	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
774	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
775	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
776	!
777	zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0
778	!
779	ELSE
780	zsnwiceload(:,:) = 0.0_wp
781	ENDIF
782	! compute equivalent ssh in lead
783	ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:)
784	!
785	END FUNCTION ice_var_sshdyn
786
787
788	!!-------------------------------------------------------------------
789	!! * INTERFACE ice_var_itd *
790	!!
791	!! ** Purpose : converting N-cat ice to jpl ice categories
792	!!-------------------------------------------------------------------
793	SUBROUTINE ice_var_itd_1c1c( phti, phts, pati , ph_i, ph_s, pa_i, &
794	& ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
795	!!-------------------------------------------------------------------
796	!! ** Purpose : converting 1-cat ice to 1 ice category
797	!!-------------------------------------------------------------------
798	REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
799	REAL(wp), DIMENSION(:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
800	REAL(wp), DIMENSION(:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds
801	REAL(wp), DIMENSION(:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds
802	!!-------------------------------------------------------------------
803	! == thickness and concentration == !
804	ph_i(:) = phti(:)
805	ph_s(:) = phts(:)
806	pa_i(:) = pati(:)
807	!
808	! == temperature and salinity and ponds == !
809	pt_i (:) = ptmi (:)
810	pt_s (:) = ptms (:)
811	pt_su(:) = ptmsu(:)
812	ps_i (:) = psmi (:)
813	pa_ip(:) = patip(:)
814	ph_ip(:) = phtip(:)
815
816	END SUBROUTINE ice_var_itd_1c1c
817
818	SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati , ph_i, ph_s, pa_i, &
819	& ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
820	!!-------------------------------------------------------------------
821	!! ** Purpose : converting N-cat ice to 1 ice category
822	!!-------------------------------------------------------------------
823	REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
824	REAL(wp), DIMENSION(:) , INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
825	REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds
826	REAL(wp), DIMENSION(:) , INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds
827	!
828	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs
829	!
830	INTEGER :: idim
831	!!-------------------------------------------------------------------
832	!
833	idim = SIZE( phti, 1 )
834	!
835	! == thickness and concentration == !
836	ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim) )
837	!
838	pa_i(:) = SUM( pati(:,:), dim=2 )
839
840	WHERE( ( pa_i(:) ) /= 0._wp ) ; z1_ai(:) = 1._wp / pa_i(:)
841	ELSEWHERE ; z1_ai(:) = 0._wp
842	END WHERE
843
844	ph_i(:) = SUM( phti(:,:) * pati(:,:), dim=2 ) * z1_ai(:)
845	ph_s(:) = SUM( phts(:,:) * pati(:,:), dim=2 ) * z1_ai(:)
846	!
847	! == temperature and salinity == !
848	WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp ) ; z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) )
849	ELSEWHERE ; z1_vi(:) = 0._wp
850	END WHERE
851	WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp ) ; z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) )
852	ELSEWHERE ; z1_vs(:) = 0._wp
853	END WHERE
854	pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
855	pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
856	pt_su(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:)
857	ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
858
859	! == ponds == !
860	pa_ip(:) = SUM( patip(:,:), dim=2 )
861	WHERE( pa_ip(:) /= 0._wp ) ; ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
862	ELSEWHERE ; ph_ip(:) = 0._wp
863	END WHERE
864	!
865	DEALLOCATE( z1_ai, z1_vi, z1_vs )
866	!
867	END SUBROUTINE ice_var_itd_Nc1c
868
869	SUBROUTINE ice_var_itd_1cMc( phti, phts, pati , ph_i, ph_s, pa_i, &
870	& ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
871	!!-------------------------------------------------------------------
872	!!
873	!! ** Purpose : converting 1-cat ice to jpl ice categories
874	!!
875	!!
876	!! ** Method: ice thickness distribution follows a gamma function from Abraham et al. (2015)
877	!! it has the property of conserving total concentration and volume
878	!!
879	!!
880	!! ** Arguments : phti: 1-cat ice thickness
881	!! phts: 1-cat snow depth
882	!! pati: 1-cat ice concentration
883	!!
884	!! ** Output : jpl-cat
885	!!
886	!! Abraham, C., Steiner, N., Monahan, A. and Michel, C., 2015.
887	!! Effects of subgrid‐scale snow thickness variability on radiative transfer in sea ice.
888	!! Journal of Geophysical Research: Oceans, 120(8), pp.5597-5614
889	!!-------------------------------------------------------------------
890	REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
891	REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
892	REAL(wp), DIMENSION(:) , INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds
893	REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds
894	!
895	REAL(wp), ALLOCATABLE, DIMENSION(:) :: zfra, z1_hti
896	INTEGER :: ji, jk, jl
897	INTEGER :: idim
898	REAL(wp) :: zv, zdh
899	!!-------------------------------------------------------------------
900	!
901	idim = SIZE( phti , 1 )
902	!
903	ph_i(1:idim,1:jpl) = 0._wp
904	ph_s(1:idim,1:jpl) = 0._wp
905	pa_i(1:idim,1:jpl) = 0._wp
906	!
907	ALLOCATE( z1_hti(idim) )
908	WHERE( phti(:) /= 0._wp ) ; z1_hti(:) = 1._wp / phti(:)
909	ELSEWHERE ; z1_hti(:) = 0._wp
910	END WHERE
911	!
912	! == thickness and concentration == !
913	! for categories 1:jpl-1, integrate the gamma function from hi_max(jl-1) to hi_max(jl)
914	DO jl = 1, jpl-1
915	DO ji = 1, idim
916	!
917	IF( phti(ji) > 0._wp ) THEN
918	! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl)
919	pa_i(ji,jl) = pati(ji) * z1_hti(ji) * ( ( phti(ji) + 2.hi_max(jl-1) ) EXP( -2.hi_max(jl-1)z1_hti(ji) ) &
920	& - ( phti(ji) + 2.hi_max(jl ) ) EXP( -2.hi_max(jl )z1_hti(ji) ) )
921	!
922	! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl)
923	zv = pati(ji) * z1_hti(ji) * ( ( phti(ji)phti(ji) + 2.phti(ji)hi_max(jl-1) + 2.hi_max(jl-1)*hi_max(jl-1) ) &
924	& * EXP( -2.hi_max(jl-1)z1_hti(ji) ) &
925	& - ( phti(ji)phti(ji) + 2.phti(ji)hi_max(jl) + 2.hi_max(jl)*hi_max(jl) ) &
926	& * EXP(-2.hi_max(jl)z1_hti(ji)) )
927	! thickness
928	IF( pa_i(ji,jl) > epsi06 ) THEN
929	ph_i(ji,jl) = zv / pa_i(ji,jl)
930	ELSE
931	ph_i(ji,jl) = 0.
932	pa_i(ji,jl) = 0.
933	ENDIF
934	ENDIF
935	!
936	ENDDO
937	ENDDO
938	!
939	! for the last category (jpl), integrate the gamma function from hi_max(jpl-1) to infinity
940	DO ji = 1, idim
941	!
942	IF( phti(ji) > 0._wp ) THEN
943	! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jpl-1) to infinity
944	pa_i(ji,jpl) = pati(ji) * z1_hti(ji) * ( phti(ji) + 2.hi_max(jpl-1) ) EXP( -2.hi_max(jpl-1)z1_hti(ji) )
945
946	! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jpl-1) to infinity
947	zv = pati(ji) * z1_hti(ji) * ( phti(ji)phti(ji) + 2.phti(ji)hi_max(jpl-1) + 2.hi_max(jpl-1)*hi_max(jpl-1) ) &
948	& * EXP( -2.hi_max(jpl-1)z1_hti(ji) )
949	! thickness
950	IF( pa_i(ji,jpl) > epsi06 ) THEN
951	ph_i(ji,jpl) = zv / pa_i(ji,jpl)
952	else
953	ph_i(ji,jpl) = 0.
954	pa_i(ji,jpl) = 0.
955	ENDIF
956	ENDIF
957	!
958	ENDDO
959	!
960	! Add Snow in each category where pa_i is not 0
961	DO jl = 1, jpl
962	DO ji = 1, idim
963	IF( pa_i(ji,jl) > 0._wp ) THEN
964	ph_s(ji,jl) = ph_i(ji,jl) * phts(ji) * z1_hti(ji)
965	! In case snow load is in excess that would lead to transformation from snow to ice
966	! Then, transfer the snow excess into the ice (different from icethd_dh)
967	zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rau0 ) * ph_i(ji,jl) ) * r1_rau0 )
968	! recompute h_i, h_s avoiding out of bounds values
969	ph_i(ji,jl) = MIN( hi_max(jl), ph_i(ji,jl) + zdh )
970	ph_s(ji,jl) = MAX( 0._wp, ph_s(ji,jl) - zdh * rhoi * r1_rhos )
971	ENDIF
972	END DO
973	END DO
974	!
975	DEALLOCATE( z1_hti )
976	!
977	! == temperature and salinity == !
978	DO jl = 1, jpl
979	pt_i (:,jl) = ptmi (:)
980	pt_s (:,jl) = ptms (:)
981	pt_su(:,jl) = ptmsu(:)
982	ps_i (:,jl) = psmi (:)
983	ps_i (:,jl) = psmi (:)
984	END DO
985	!
986	! == ponds == !
987	ALLOCATE( zfra(idim) )
988	! keep the same pond fraction atip/ati for each category
989	WHERE( pati(:) /= 0._wp ) ; zfra(:) = patip(:) / pati(:)
990	ELSEWHERE ; zfra(:) = 0._wp
991	END WHERE
992	DO jl = 1, jpl
993	pa_ip(:,jl) = zfra(:) * pa_i(:,jl)
994	END DO
995	! keep the same v_ip/v_i ratio for each category
996	WHERE( ( phti(:) * pati(:) ) /= 0._wp ) ; zfra(:) = ( phtip(:) * patip(:) ) / ( phti(:) * pati(:) )
997	ELSEWHERE ; zfra(:) = 0._wp
998	END WHERE
999	DO jl = 1, jpl
1000	WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1001	ELSEWHERE ; ph_ip(:,jl) = 0._wp
1002	END WHERE
1003	END DO
1004	DEALLOCATE( zfra )
1005	!
1006	END SUBROUTINE ice_var_itd_1cMc
1007
1008	SUBROUTINE ice_var_itd_NcMc( phti, phts, pati , ph_i, ph_s, pa_i, &
1009	& ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
1010	!!-------------------------------------------------------------------
1011	!!
1012	!! ** Purpose : converting N-cat ice to jpl ice categories
1013	!!
1014	!! ice thickness distribution follows a gaussian law
1015	!! around the concentration of the most likely ice thickness
1016	!! (similar as iceistate.F90)
1017	!!
1018	!! ** Method: Iterative procedure
1019	!!
1020	!! 1) Fill ice cat that correspond to input thicknesses
1021	!! Find the lowest(jlmin) and highest(jlmax) cat that are filled
1022	!!
1023	!! 2) Expand the filling to the cat jlmin-1 and jlmax+1
1024	!! by removing 25% ice area from jlmin and jlmax (resp.)
1025	!!
1026	!! 3) Expand the filling to the empty cat between jlmin and jlmax
1027	!! by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
1028	!! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
1029	!!
1030	!! ** Arguments : phti: N-cat ice thickness
1031	!! phts: N-cat snow depth
1032	!! pati: N-cat ice concentration
1033	!!
1034	!! ** Output : jpl-cat
1035	!!
1036	!! (Example of application: BDY forcings when inputs have N-cat /= jpl)
1037	!!-------------------------------------------------------------------
1038	REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
1039	REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
1040	REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds
1041	REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds
1042	!
1043	INTEGER , ALLOCATABLE, DIMENSION(:,:) :: jlfil, jlfil2
1044	INTEGER , ALLOCATABLE, DIMENSION(:) :: jlmax, jlmin
1045	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs, ztmp, zfra
1046	!
1047	REAL(wp), PARAMETER :: ztrans = 0.25_wp
1048	INTEGER :: ji, jl, jl1, jl2
1049	INTEGER :: idim, icat
1050	!!-------------------------------------------------------------------
1051	!
1052	idim = SIZE( phti, 1 )
1053	icat = SIZE( phti, 2 )
1054	!
1055	! == thickness and concentration == !
1056	! ! ---------------------- !
1057	IF( icat == jpl ) THEN ! input cat = output cat !
1058	! ! ---------------------- !
1059	ph_i(:,:) = phti(:,:)
1060	ph_s(:,:) = phts(:,:)
1061	pa_i(:,:) = pati(:,:)
1062	!
1063	! == temperature and salinity and ponds == !
1064	pt_i (:,:) = ptmi (:,:)
1065	pt_s (:,:) = ptms (:,:)
1066	pt_su(:,:) = ptmsu(:,:)
1067	ps_i (:,:) = psmi (:,:)
1068	pa_ip(:,:) = patip(:,:)
1069	ph_ip(:,:) = phtip(:,:)
1070	! ! ---------------------- !
1071	ELSEIF( icat == 1 ) THEN ! input cat = 1 !
1072	! ! ---------------------- !
1073	CALL ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), &
1074	& ph_i(:,:), ph_s(:,:), pa_i (:,:), &
1075	& ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), &
1076	& pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:) )
1077	! ! ---------------------- !
1078	ELSEIF( jpl == 1 ) THEN ! output cat = 1 !
1079	! ! ---------------------- !
1080	CALL ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), &
1081	& ph_i(:,1), ph_s(:,1), pa_i (:,1), &
1082	& ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), &
1083	& pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1) )
1084	! ! ----------------------- !
1085	ELSE ! input cat /= output cat !
1086	! ! ----------------------- !
1087
1088	ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays
1089	ALLOCATE( jlmin(idim), jlmax(idim) )
1090
1091	! --- initialize output fields to 0 --- !
1092	ph_i(1:idim,1:jpl) = 0._wp
1093	ph_s(1:idim,1:jpl) = 0._wp
1094	pa_i(1:idim,1:jpl) = 0._wp
1095	!
1096	! --- fill the categories --- !
1097	! find where cat-input = cat-output and fill cat-output fields
1098	jlmax(:) = 0
1099	jlmin(:) = 999
1100	jlfil(:,:) = 0
1101	DO jl1 = 1, jpl
1102	DO jl2 = 1, icat
1103	DO ji = 1, idim
1104	IF( hi_max(jl1-1) <= phti(ji,jl2) .AND. hi_max(jl1) > phti(ji,jl2) ) THEN
1105	! fill the right category
1106	ph_i(ji,jl1) = phti(ji,jl2)
1107	ph_s(ji,jl1) = phts(ji,jl2)
1108	pa_i(ji,jl1) = pati(ji,jl2)
1109	! record categories that are filled
1110	jlmax(ji) = MAX( jlmax(ji), jl1 )
1111	jlmin(ji) = MIN( jlmin(ji), jl1 )
1112	jlfil(ji,jl1) = jl1
1113	ENDIF
1114	END DO
1115	END DO
1116	END DO
1117	!
1118	! --- fill the gaps between categories --- !
1119	! transfer from categories filled at the previous step to the empty ones in between
1120	DO ji = 1, idim
1121	jl1 = jlmin(ji)
1122	jl2 = jlmax(ji)
1123	IF( jl1 > 1 ) THEN
1124	! fill the lower cat (jl1-1)
1125	pa_i(ji,jl1-1) = ztrans * pa_i(ji,jl1)
1126	ph_i(ji,jl1-1) = hi_mean(jl1-1)
1127	! remove from cat jl1
1128	pa_i(ji,jl1 ) = ( 1._wp - ztrans ) * pa_i(ji,jl1)
1129	ENDIF
1130	IF( jl2 < jpl ) THEN
1131	! fill the upper cat (jl2+1)
1132	pa_i(ji,jl2+1) = ztrans * pa_i(ji,jl2)
1133	ph_i(ji,jl2+1) = hi_mean(jl2+1)
1134	! remove from cat jl2
1135	pa_i(ji,jl2 ) = ( 1._wp - ztrans ) * pa_i(ji,jl2)
1136	ENDIF
1137	END DO
1138	!
1139	jlfil2(:,:) = jlfil(:,:)
1140	! fill categories from low to high
1141	DO jl = 2, jpl-1
1142	DO ji = 1, idim
1143	IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
1144	! fill high
1145	pa_i(ji,jl) = ztrans * pa_i(ji,jl-1)
1146	ph_i(ji,jl) = hi_mean(jl)
1147	jlfil(ji,jl) = jl
1148	! remove low
1149	pa_i(ji,jl-1) = ( 1._wp - ztrans ) * pa_i(ji,jl-1)
1150	ENDIF
1151	END DO
1152	END DO
1153	!
1154	! fill categories from high to low
1155	DO jl = jpl-1, 2, -1
1156	DO ji = 1, idim
1157	IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
1158	! fill low
1159	pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1)
1160	ph_i(ji,jl) = hi_mean(jl)
1161	jlfil2(ji,jl) = jl
1162	! remove high
1163	pa_i(ji,jl+1) = ( 1._wp - ztrans ) * pa_i(ji,jl+1)
1164	ENDIF
1165	END DO
1166	END DO
1167	!
1168	DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays
1169	DEALLOCATE( jlmin, jlmax )
1170	!
1171	! == temperature and salinity == !
1172	!
1173	ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) )
1174	!
1175	WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp ) ; z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 )
1176	ELSEWHERE ; z1_ai(:) = 0._wp
1177	END WHERE
1178	WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp ) ; z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 )
1179	ELSEWHERE ; z1_vi(:) = 0._wp
1180	END WHERE
1181	WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp ) ; z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 )
1182	ELSEWHERE ; z1_vs(:) = 0._wp
1183	END WHERE
1184	!
1185	! fill all the categories with the same value
1186	ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
1187	DO jl = 1, jpl
1188	pt_i (:,jl) = ztmp(:)
1189	END DO
1190	ztmp(:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
1191	DO jl = 1, jpl
1192	pt_s (:,jl) = ztmp(:)
1193	END DO
1194	ztmp(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:)
1195	DO jl = 1, jpl
1196	pt_su(:,jl) = ztmp(:)
1197	END DO
1198	ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
1199	DO jl = 1, jpl
1200	ps_i (:,jl) = ztmp(:)
1201	END DO
1202	!
1203	DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp )
1204	!
1205	! == ponds == !
1206	ALLOCATE( zfra(idim) )
1207	! keep the same pond fraction atip/ati for each category
1208	WHERE( SUM( pati(:,:), dim=2 ) /= 0._wp ) ; zfra(:) = SUM( patip(:,:), dim=2 ) / SUM( pati(:,:), dim=2 )
1209	ELSEWHERE ; zfra(:) = 0._wp
1210	END WHERE
1211	DO jl = 1, jpl
1212	pa_ip(:,jl) = zfra(:) * pa_i(:,jl)
1213	END DO
1214	! keep the same v_ip/v_i ratio for each category
1215	WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp )
1216	zfra(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 )
1217	ELSEWHERE
1218	zfra(:) = 0._wp
1219	END WHERE
1220	DO jl = 1, jpl
1221	WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1222	ELSEWHERE ; ph_ip(:,jl) = 0._wp
1223	END WHERE
1224	END DO
1225	DEALLOCATE( zfra )
1226	!
1227	ENDIF
1228	!
1229	END SUBROUTINE ice_var_itd_NcMc
1230
1231	#else
1232	!!----------------------------------------------------------------------
1233	!! Default option Dummy module NO SI3 sea-ice model
1234	!!----------------------------------------------------------------------
1235	#endif
1236
1237	!!======================================================================
1238	END MODULE icevar

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icevar.F90

Download in other formats: