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.

icevar.F90 in NEMO/trunk/src/ICE – NEMO

Context Navigation

source: NEMO/trunk/src/ICE/icevar.F90 @ 15070

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

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

Download in other formats: