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/branches/2020/r12377_ticket2386/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/r12377_ticket2386/src/ICE/icevar.F90 @ 13540

Last change on this file since 13540 was 13540, checked in by andmirek, 4 years ago
Ticket #2386: update to latest trunk
Property svn:keywords set to `Id`
File size: 64.0 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	& ( h_il(:,:,:) - zhl_min ) / ( 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( 1, 1, 1, 1, 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) :: zsal, z1_dS
344	REAL(wp) :: zargtemp , 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( 1, 1, 1, 1 )
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( 1, 1, 1, 1, 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	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
412	sz_i(:,:,jk,jl) = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) )
413	END DO
414	END DO
415	!
416	END SELECT
417	!
418	END SUBROUTINE ice_var_salprof
419
420
421	SUBROUTINE ice_var_salprof1d
422	!!-------------------------------------------------------------------
423	!! * ROUTINE ice_var_salprof1d *
424	!!
425	!! ** Purpose : 1d computation of the sea ice salinity profile
426	!! Works with 1d vectors and is used by thermodynamic modules
427	!!-------------------------------------------------------------------
428	INTEGER :: ji, jk ! dummy loop indices
429	REAL(wp) :: zargtemp, zsal, z1_dS ! local scalars
430	REAL(wp) :: zs, zs0 ! - -
431	!
432	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_slope_s, zalpha !
433	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
434	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
435	!!-------------------------------------------------------------------
436	!
437	SELECT CASE ( nn_icesal )
438	!
439	! !---------------------------------------!
440	CASE( 1 ) ! constant salinity in time and space !
441	! !---------------------------------------!
442	sz_i_1d(1:npti,:) = rn_icesal
443	!
444	! !---------------------------------------------!
445	CASE( 2 ) ! time varying salinity with linear profile !
446	! !---------------------------------------------!
447	!
448	ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
449	!
450	! ! Slope of the linear profile
451	WHERE( h_i_1d(1:npti) > epsi20 ) ; z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
452	ELSEWHERE ; z_slope_s(1:npti) = 0._wp
453	END WHERE
454
455	z1_dS = 1._wp / ( zsi1 - zsi0 )
456	DO ji = 1, npti
457	zalpha(ji) = MAX( 0._wp , MIN( ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp ) )
458	! ! force a constant profile when SSS too low (Baltic Sea)
459	IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) ) zalpha(ji) = 0._wp
460	END DO
461	!
462	! Computation of the profile
463	DO jk = 1, nlay_i
464	DO ji = 1, npti
465	! ! linear profile with 0 surface value
466	zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
467	zs = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * s_i_1d(ji)
468	sz_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
469	END DO
470	END DO
471	!
472	DEALLOCATE( z_slope_s, zalpha )
473
474	! !-------------------------------------------!
475	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
476	! !-------------------------------------------! (mean = 2.30)
477	!
478	s_i_1d(1:npti) = 2.30_wp
479	!
480	!!gm cf remark in ice_var_salprof routine, CASE( 3 )
481	DO jk = 1, nlay_i
482	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
483	zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
484	DO ji = 1, npti
485	sz_i_1d(ji,jk) = zsal
486	END DO
487	END DO
488	!
489	END SELECT
490	!
491	END SUBROUTINE ice_var_salprof1d
492
493
494	SUBROUTINE ice_var_zapsmall
495	!!-------------------------------------------------------------------
496	!! * ROUTINE ice_var_zapsmall *
497	!!
498	!! ** Purpose : Remove too small sea ice areas and correct fluxes
499	!!-------------------------------------------------------------------
500	INTEGER :: ji, jj, jl, jk ! dummy loop indices
501	REAL(wp), DIMENSION(jpi,jpj) :: zswitch
502	!!-------------------------------------------------------------------
503	!
504	DO jl = 1, jpl !== loop over the categories ==!
505	!
506	WHERE( a_i(:,:,jl) > epsi10 ) ; h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
507	ELSEWHERE ; h_i(:,:,jl) = 0._wp
508	END WHERE
509	!
510	WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 ) ; zswitch(:,:) = 0._wp
511	ELSEWHERE ; zswitch(:,:) = 1._wp
512	END WHERE
513	!
514	!-----------------------------------------------------------------
515	! Zap ice energy and use ocean heat to melt ice
516	!-----------------------------------------------------------------
517	DO_3D( 1, 1, 1, 1, 1, nlay_i )
518	! update exchanges with ocean
519	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
520	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
521	t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
522	END_3D
523	!
524	DO_3D( 1, 1, 1, 1, 1, nlay_s )
525	! update exchanges with ocean
526	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
527	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj)
528	t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
529	END_3D
530	!
531	!-----------------------------------------------------------------
532	! zap ice and snow volume, add water and salt to ocean
533	!-----------------------------------------------------------------
534	DO_2D( 1, 1, 1, 1 )
535	! update exchanges with ocean
536	sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoi * r1_Dt_ice
537	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoi * r1_Dt_ice
538	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhos * r1_Dt_ice
539	!
540	a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
541	v_i (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
542	v_s (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
543	t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
544	oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
545	sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
546	!
547	h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
548	h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
549	!
550	a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
551	v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
552	v_il (ji,jj,jl) = v_il (ji,jj,jl) * zswitch(ji,jj)
553	!
554	END_2D
555	!
556	END DO
557
558	! to be sure that at_i is the sum of a_i(jl)
559	at_i (:,:) = SUM( a_i (:,:,:), dim=3 )
560	vt_i (:,:) = SUM( v_i (:,:,:), dim=3 )
561	!!clem add?
562	! vt_s (:,:) = SUM( v_s (:,:,:), dim=3 )
563	! st_i (:,:) = SUM( sv_i(:,:,:), dim=3 )
564	! et_s(:,:) = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 )
565	! et_i(:,:) = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 )
566	!!clem
567
568	! open water = 1 if at_i=0
569	WHERE( at_i(:,:) == 0._wp ) ato_i(:,:) = 1._wp
570	!
571	END SUBROUTINE ice_var_zapsmall
572
573
574	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 )
575	!!-------------------------------------------------------------------
576	!! * ROUTINE ice_var_zapneg *
577	!!
578	!! ** Purpose : Remove negative sea ice fields and correct fluxes
579	!!-------------------------------------------------------------------
580	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
581	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area
582	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume
583	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_s ! snw volume
584	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psv_i ! salt content
585	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content
586	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration
587	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
588	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume
589	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_il ! melt pond lid volume
590	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content
591	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content
592	!
593	INTEGER :: ji, jj, jl, jk ! dummy loop indices
594	REAL(wp) :: z1_dt
595	!!-------------------------------------------------------------------
596	!
597	z1_dt = 1._wp / pdt
598	!
599	DO jl = 1, jpl !== loop over the categories ==!
600	!
601	! make sure a_i=0 where v_i<=0
602	WHERE( pv_i(:,:,:) <= 0._wp ) pa_i(:,:,:) = 0._wp
603
604	!----------------------------------------
605	! zap ice energy and send it to the ocean
606	!----------------------------------------
607	DO_3D( 1, 1, 1, 1, 1, nlay_i )
608	IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
609	hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
610	pe_i(ji,jj,jk,jl) = 0._wp
611	ENDIF
612	END_3D
613	!
614	DO_3D( 1, 1, 1, 1, 1, nlay_s )
615	IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
616	hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
617	pe_s(ji,jj,jk,jl) = 0._wp
618	ENDIF
619	END_3D
620	!
621	!-----------------------------------------------------
622	! zap ice and snow volume, add water and salt to ocean
623	!-----------------------------------------------------
624	DO_2D( 1, 1, 1, 1 )
625	IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
626	wfx_res(ji,jj) = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
627	pv_i (ji,jj,jl) = 0._wp
628	ENDIF
629	IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
630	wfx_res(ji,jj) = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt
631	pv_s (ji,jj,jl) = 0._wp
632	ENDIF
633	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
634	sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt
635	psv_i (ji,jj,jl) = 0._wp
636	ENDIF
637	END_2D
638	!
639	END DO
640	!
641	WHERE( pato_i(:,:) < 0._wp ) pato_i(:,:) = 0._wp
642	WHERE( poa_i (:,:,:) < 0._wp ) poa_i (:,:,:) = 0._wp
643	WHERE( pa_i (:,:,:) < 0._wp ) pa_i (:,:,:) = 0._wp
644	WHERE( pa_ip (:,:,:) < 0._wp ) pa_ip (:,:,:) = 0._wp
645	WHERE( pv_ip (:,:,:) < 0._wp ) pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
646	WHERE( pv_il (:,:,:) < 0._wp ) pv_il (:,:,:) = 0._wp ! but it does not change conservation, so keep it this way is ok
647	!
648	END SUBROUTINE ice_var_zapneg
649
650
651	SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
652	!!-------------------------------------------------------------------
653	!! * ROUTINE ice_var_roundoff *
654	!!
655	!! ** Purpose : Remove negative sea ice values arising from roundoff errors
656	!!-------------------------------------------------------------------
657	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_i ! ice concentration
658	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_i ! ice volume
659	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_s ! snw volume
660	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: psv_i ! salt content
661	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: poa_i ! age content
662	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
663	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_ip ! melt pond volume
664	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_il ! melt pond lid volume
665	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_s ! snw heat content
666	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_i ! ice heat content
667	!!-------------------------------------------------------------------
668	!
669
670	WHERE( pa_i (1:npti,:) < 0._wp ) pa_i (1:npti,:) = 0._wp ! a_i must be >= 0
671	WHERE( pv_i (1:npti,:) < 0._wp ) pv_i (1:npti,:) = 0._wp ! v_i must be >= 0
672	WHERE( pv_s (1:npti,:) < 0._wp ) pv_s (1:npti,:) = 0._wp ! v_s must be >= 0
673	WHERE( psv_i(1:npti,:) < 0._wp ) psv_i(1:npti,:) = 0._wp ! sv_i must be >= 0
674	WHERE( poa_i(1:npti,:) < 0._wp ) poa_i(1:npti,:) = 0._wp ! oa_i must be >= 0
675	WHERE( pe_i (1:npti,:,:) < 0._wp ) pe_i (1:npti,:,:) = 0._wp ! e_i must be >= 0
676	WHERE( pe_s (1:npti,:,:) < 0._wp ) pe_s (1:npti,:,:) = 0._wp ! e_s must be >= 0
677	IF( ln_pnd_LEV ) THEN
678	WHERE( pa_ip(1:npti,:) < 0._wp ) pa_ip(1:npti,:) = 0._wp ! a_ip must be >= 0
679	WHERE( pv_ip(1:npti,:) < 0._wp ) pv_ip(1:npti,:) = 0._wp ! v_ip must be >= 0
680	IF( ln_pnd_lids ) THEN
681	WHERE( pv_il(1:npti,:) < 0._wp .AND. pv_il(1:npti,:) > -epsi10 ) pv_il(1:npti,:) = 0._wp ! v_il must be >= 0
682	ENDIF
683	ENDIF
684	!
685	END SUBROUTINE ice_var_roundoff
686
687
688	SUBROUTINE ice_var_bv
689	!!-------------------------------------------------------------------
690	!! * ROUTINE ice_var_bv *
691	!!
692	!! ** Purpose : computes mean brine volume (%) in sea ice
693	!!
694	!! ** Method : e = - 0.054 * S (ppt) / T (C)
695	!!
696	!! References : Vancoppenolle et al., JGR, 2007
697	!!-------------------------------------------------------------------
698	INTEGER :: ji, jj, jk, jl ! dummy loop indices
699	!!-------------------------------------------------------------------
700	!
701	!!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done
702	!! instead of setting everything to zero as just below
703	bv_i (:,:,:) = 0._wp
704	DO jl = 1, jpl
705	DO jk = 1, nlay_i
706	WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )
707	bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
708	END WHERE
709	END DO
710	END DO
711	WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
712	ELSEWHERE ; bvm_i(:,:) = 0._wp
713	END WHERE
714	!
715	END SUBROUTINE ice_var_bv
716
717
718	SUBROUTINE ice_var_enthalpy
719	!!-------------------------------------------------------------------
720	!! * ROUTINE ice_var_enthalpy *
721	!!
722	!! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature
723	!!
724	!! ** Method : Formula (Bitz and Lipscomb, 1999)
725	!!-------------------------------------------------------------------
726	INTEGER :: ji, jk ! dummy loop indices
727	REAL(wp) :: ztmelts ! local scalar
728	!!-------------------------------------------------------------------
729	!
730	DO jk = 1, nlay_i ! Sea ice energy of melting
731	DO ji = 1, npti
732	ztmelts = - rTmlt * sz_i_1d(ji,jk)
733	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
734	! (sometimes zdf scheme produces abnormally high temperatures)
735	e_i_1d(ji,jk) = rhoi * ( rcpi * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) ) &
736	& + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) ) &
737	& - rcp * ztmelts )
738	END DO
739	END DO
740	DO jk = 1, nlay_s ! Snow energy of melting
741	DO ji = 1, npti
742	e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
743	END DO
744	END DO
745	!
746	END SUBROUTINE ice_var_enthalpy
747
748
749	FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
750	!!---------------------------------------------------------------------
751	!! * ROUTINE ice_var_sshdyn *
752	!!
753	!! ** Purpose : compute the equivalent ssh in lead when sea ice is embedded
754	!!
755	!! ** Method : ssh_lead = ssh + (Mice + Msnow) / rho0
756	!!
757	!! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira,
758	!! Sea ice-ocean coupling using a rescaled vertical coordinate z*,
759	!! Ocean Modelling, Volume 24, Issues 1-2, 2008
760	!!----------------------------------------------------------------------
761	!
762	! input
763	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh !: ssh [m]
764	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass !: mass of snow and ice at current ice time step [Kg/m2]
765	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b !: mass of snow and ice at previous ice time step [Kg/m2]
766	!
767	! output
768	REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn ! equivalent ssh in lead [m]
769	!
770	! temporary
771	REAL(wp) :: zintn, zintb ! time interpolation weights []
772	REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload ! snow and ice load [m]
773	!
774	! compute ice load used to define the equivalent ssh in lead
775	IF( ln_ice_embd ) THEN
776	!
777	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
778	! = (1/nn_fsbc)^2 * {SUM[n] , n=0,nn_fsbc-1}
779	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
780	!
781	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
782	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
783	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
784	!
785	zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rho0
786	!
787	ELSE
788	zsnwiceload(:,:) = 0.0_wp
789	ENDIF
790	! compute equivalent ssh in lead
791	ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:)
792	!
793	END FUNCTION ice_var_sshdyn
794
795
796	!!-------------------------------------------------------------------
797	!! * INTERFACE ice_var_itd *
798	!!
799	!! ** Purpose : converting N-cat ice to jpl ice categories
800	!!-------------------------------------------------------------------
801	SUBROUTINE ice_var_itd_1c1c( phti, phts, pati , ph_i, ph_s, pa_i, &
802	& ptmi, ptms, ptmsu, psmi, patip, phtip, phtil, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
803	!!-------------------------------------------------------------------
804	!! ** Purpose : converting 1-cat ice to 1 ice category
805	!!-------------------------------------------------------------------
806	REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
807	REAL(wp), DIMENSION(:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
808	REAL(wp), DIMENSION(:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, phtil ! input ice/snow temp & sal & ponds
809	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
810	!!-------------------------------------------------------------------
811	! == thickness and concentration == !
812	ph_i(:) = phti(:)
813	ph_s(:) = phts(:)
814	pa_i(:) = pati(:)
815	!
816	! == temperature and salinity and ponds == !
817	pt_i (:) = ptmi (:)
818	pt_s (:) = ptms (:)
819	pt_su(:) = ptmsu(:)
820	ps_i (:) = psmi (:)
821	pa_ip(:) = patip(:)
822	ph_ip(:) = phtip(:)
823	ph_il(:) = phtil(:)
824
825	END SUBROUTINE ice_var_itd_1c1c
826
827	SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati , ph_i, ph_s, pa_i, &
828	& ptmi, ptms, ptmsu, psmi, patip, phtip, phtil, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
829	!!-------------------------------------------------------------------
830	!! ** Purpose : converting N-cat ice to 1 ice category
831	!!-------------------------------------------------------------------
832	REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
833	REAL(wp), DIMENSION(:) , INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
834	REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, phtil ! input ice/snow temp & sal & ponds
835	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
836	!
837	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs
838	!
839	INTEGER :: idim
840	!!-------------------------------------------------------------------
841	!
842	idim = SIZE( phti, 1 )
843	!
844	! == thickness and concentration == !
845	ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim) )
846	!
847	pa_i(:) = SUM( pati(:,:), dim=2 )
848
849	WHERE( ( pa_i(:) ) /= 0._wp ) ; z1_ai(:) = 1._wp / pa_i(:)
850	ELSEWHERE ; z1_ai(:) = 0._wp
851	END WHERE
852
853	ph_i(:) = SUM( phti(:,:) * pati(:,:), dim=2 ) * z1_ai(:)
854	ph_s(:) = SUM( phts(:,:) * pati(:,:), dim=2 ) * z1_ai(:)
855	!
856	! == temperature and salinity == !
857	WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp ) ; z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) )
858	ELSEWHERE ; z1_vi(:) = 0._wp
859	END WHERE
860	WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp ) ; z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) )
861	ELSEWHERE ; z1_vs(:) = 0._wp
862	END WHERE
863	pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
864	pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
865	pt_su(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:)
866	ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
867
868	! == ponds == !
869	pa_ip(:) = SUM( patip(:,:), dim=2 )
870	WHERE( pa_ip(:) /= 0._wp )
871	ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
872	ph_il(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
873	ELSEWHERE
874	ph_ip(:) = 0._wp
875	ph_il(:) = 0._wp
876	END WHERE
877	!
878	DEALLOCATE( z1_ai, z1_vi, z1_vs )
879	!
880	END SUBROUTINE ice_var_itd_Nc1c
881
882	SUBROUTINE ice_var_itd_1cMc( phti, phts, pati , ph_i, ph_s, pa_i, &
883	& ptmi, ptms, ptmsu, psmi, patip, phtip, phtil, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
884	!!-------------------------------------------------------------------
885	!!
886	!! ** Purpose : converting 1-cat ice to jpl ice categories
887	!!
888	!!
889	!! ** Method: ice thickness distribution follows a gamma function from Abraham et al. (2015)
890	!! it has the property of conserving total concentration and volume
891	!!
892	!!
893	!! ** Arguments : phti: 1-cat ice thickness
894	!! phts: 1-cat snow depth
895	!! pati: 1-cat ice concentration
896	!!
897	!! ** Output : jpl-cat
898	!!
899	!! Abraham, C., Steiner, N., Monahan, A. and Michel, C., 2015.
900	!! Effects of subgrid‐scale snow thickness variability on radiative transfer in sea ice.
901	!! Journal of Geophysical Research: Oceans, 120(8), pp.5597-5614
902	!!-------------------------------------------------------------------
903	REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
904	REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
905	REAL(wp), DIMENSION(:) , INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, phtil ! input ice/snow temp & sal & ponds
906	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
907	!
908	REAL(wp), ALLOCATABLE, DIMENSION(:) :: zfra, z1_hti
909	INTEGER :: ji, jk, jl
910	INTEGER :: idim
911	REAL(wp) :: zv, zdh
912	!!-------------------------------------------------------------------
913	!
914	idim = SIZE( phti , 1 )
915	!
916	ph_i(1:idim,1:jpl) = 0._wp
917	ph_s(1:idim,1:jpl) = 0._wp
918	pa_i(1:idim,1:jpl) = 0._wp
919	!
920	ALLOCATE( z1_hti(idim) )
921	WHERE( phti(:) /= 0._wp ) ; z1_hti(:) = 1._wp / phti(:)
922	ELSEWHERE ; z1_hti(:) = 0._wp
923	END WHERE
924	!
925	! == thickness and concentration == !
926	! for categories 1:jpl-1, integrate the gamma function from hi_max(jl-1) to hi_max(jl)
927	DO jl = 1, jpl-1
928	DO ji = 1, idim
929	!
930	IF( phti(ji) > 0._wp ) THEN
931	! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl)
932	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) ) &
933	& - ( phti(ji) + 2.hi_max(jl ) ) EXP( -2.hi_max(jl )z1_hti(ji) ) )
934	!
935	! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl)
936	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) ) &
937	& * EXP( -2.hi_max(jl-1)z1_hti(ji) ) &
938	& - ( phti(ji)phti(ji) + 2.phti(ji)hi_max(jl) + 2.hi_max(jl)*hi_max(jl) ) &
939	& * EXP(-2.hi_max(jl)z1_hti(ji)) )
940	! thickness
941	IF( pa_i(ji,jl) > epsi06 ) THEN
942	ph_i(ji,jl) = zv / pa_i(ji,jl)
943	ELSE
944	ph_i(ji,jl) = 0.
945	pa_i(ji,jl) = 0.
946	ENDIF
947	ENDIF
948	!
949	ENDDO
950	ENDDO
951	!
952	! for the last category (jpl), integrate the gamma function from hi_max(jpl-1) to infinity
953	DO ji = 1, idim
954	!
955	IF( phti(ji) > 0._wp ) THEN
956	! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jpl-1) to infinity
957	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) )
958
959	! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jpl-1) to infinity
960	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) ) &
961	& * EXP( -2.hi_max(jpl-1)z1_hti(ji) )
962	! thickness
963	IF( pa_i(ji,jpl) > epsi06 ) THEN
964	ph_i(ji,jpl) = zv / pa_i(ji,jpl)
965	else
966	ph_i(ji,jpl) = 0.
967	pa_i(ji,jpl) = 0.
968	ENDIF
969	ENDIF
970	!
971	ENDDO
972	!
973	! Add Snow in each category where pa_i is not 0
974	DO jl = 1, jpl
975	DO ji = 1, idim
976	IF( pa_i(ji,jl) > 0._wp ) THEN
977	ph_s(ji,jl) = ph_i(ji,jl) * phts(ji) * z1_hti(ji)
978	! In case snow load is in excess that would lead to transformation from snow to ice
979	! Then, transfer the snow excess into the ice (different from icethd_dh)
980	zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rho0 ) * ph_i(ji,jl) ) * r1_rho0 )
981	! recompute h_i, h_s avoiding out of bounds values
982	ph_i(ji,jl) = MIN( hi_max(jl), ph_i(ji,jl) + zdh )
983	ph_s(ji,jl) = MAX( 0._wp, ph_s(ji,jl) - zdh * rhoi * r1_rhos )
984	ENDIF
985	END DO
986	END DO
987	!
988	DEALLOCATE( z1_hti )
989	!
990	! == temperature and salinity == !
991	DO jl = 1, jpl
992	pt_i (:,jl) = ptmi (:)
993	pt_s (:,jl) = ptms (:)
994	pt_su(:,jl) = ptmsu(:)
995	ps_i (:,jl) = psmi (:)
996	END DO
997	!
998	! == ponds == !
999	ALLOCATE( zfra(idim) )
1000	! keep the same pond fraction atip/ati for each category
1001	WHERE( pati(:) /= 0._wp ) ; zfra(:) = patip(:) / pati(:)
1002	ELSEWHERE ; zfra(:) = 0._wp
1003	END WHERE
1004	DO jl = 1, jpl
1005	pa_ip(:,jl) = zfra(:) * pa_i(:,jl)
1006	END DO
1007	! keep the same v_ip/v_i ratio for each category
1008	WHERE( ( phti(:) * pati(:) ) /= 0._wp ) ; zfra(:) = ( phtip(:) * patip(:) ) / ( phti(:) * pati(:) )
1009	ELSEWHERE ; zfra(:) = 0._wp
1010	END WHERE
1011	DO jl = 1, jpl
1012	WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1013	ELSEWHERE ; ph_ip(:,jl) = 0._wp
1014	END WHERE
1015	END DO
1016	! keep the same v_il/v_i ratio for each category
1017	WHERE( ( phti(:) * pati(:) ) /= 0._wp ) ; zfra(:) = ( phtil(:) * patip(:) ) / ( phti(:) * pati(:) )
1018	ELSEWHERE ; zfra(:) = 0._wp
1019	END WHERE
1020	DO jl = 1, jpl
1021	WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1022	ELSEWHERE ; ph_il(:,jl) = 0._wp
1023	END WHERE
1024	END DO
1025	DEALLOCATE( zfra )
1026	!
1027	END SUBROUTINE ice_var_itd_1cMc
1028
1029	SUBROUTINE ice_var_itd_NcMc( phti, phts, pati , ph_i, ph_s, pa_i, &
1030	& ptmi, ptms, ptmsu, psmi, patip, phtip, phtil, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
1031	!!-------------------------------------------------------------------
1032	!!
1033	!! ** Purpose : converting N-cat ice to jpl ice categories
1034	!!
1035	!! ice thickness distribution follows a gaussian law
1036	!! around the concentration of the most likely ice thickness
1037	!! (similar as iceistate.F90)
1038	!!
1039	!! ** Method: Iterative procedure
1040	!!
1041	!! 1) Fill ice cat that correspond to input thicknesses
1042	!! Find the lowest(jlmin) and highest(jlmax) cat that are filled
1043	!!
1044	!! 2) Expand the filling to the cat jlmin-1 and jlmax+1
1045	!! by removing 25% ice area from jlmin and jlmax (resp.)
1046	!!
1047	!! 3) Expand the filling to the empty cat between jlmin and jlmax
1048	!! by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
1049	!! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
1050	!!
1051	!! ** Arguments : phti: N-cat ice thickness
1052	!! phts: N-cat snow depth
1053	!! pati: N-cat ice concentration
1054	!!
1055	!! ** Output : jpl-cat
1056	!!
1057	!! (Example of application: BDY forcings when inputs have N-cat /= jpl)
1058	!!-------------------------------------------------------------------
1059	REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
1060	REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
1061	REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, phtil ! input ice/snow temp & sal & ponds
1062	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
1063	!
1064	INTEGER , ALLOCATABLE, DIMENSION(:,:) :: jlfil, jlfil2
1065	INTEGER , ALLOCATABLE, DIMENSION(:) :: jlmax, jlmin
1066	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs, ztmp, zfra
1067	!
1068	REAL(wp), PARAMETER :: ztrans = 0.25_wp
1069	INTEGER :: ji, jl, jl1, jl2
1070	INTEGER :: idim, icat
1071	!!-------------------------------------------------------------------
1072	!
1073	idim = SIZE( phti, 1 )
1074	icat = SIZE( phti, 2 )
1075	!
1076	! == thickness and concentration == !
1077	! ! ---------------------- !
1078	IF( icat == jpl ) THEN ! input cat = output cat !
1079	! ! ---------------------- !
1080	ph_i(:,:) = phti(:,:)
1081	ph_s(:,:) = phts(:,:)
1082	pa_i(:,:) = pati(:,:)
1083	!
1084	! == temperature and salinity and ponds == !
1085	pt_i (:,:) = ptmi (:,:)
1086	pt_s (:,:) = ptms (:,:)
1087	pt_su(:,:) = ptmsu(:,:)
1088	ps_i (:,:) = psmi (:,:)
1089	pa_ip(:,:) = patip(:,:)
1090	ph_ip(:,:) = phtip(:,:)
1091	ph_il(:,:) = phtil(:,:)
1092	! ! ---------------------- !
1093	ELSEIF( icat == 1 ) THEN ! input cat = 1 !
1094	! ! ---------------------- !
1095	CALL ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), &
1096	& ph_i(:,:), ph_s(:,:), pa_i (:,:), &
1097	& ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), phtil(:,1), &
1098	& pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:), ph_il(:,:) )
1099	! ! ---------------------- !
1100	ELSEIF( jpl == 1 ) THEN ! output cat = 1 !
1101	! ! ---------------------- !
1102	CALL ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), &
1103	& ph_i(:,1), ph_s(:,1), pa_i (:,1), &
1104	& ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), phtil(:,:), &
1105	& pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1), ph_il(:,1) )
1106	! ! ----------------------- !
1107	ELSE ! input cat /= output cat !
1108	! ! ----------------------- !
1109
1110	ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays
1111	ALLOCATE( jlmin(idim), jlmax(idim) )
1112
1113	! --- initialize output fields to 0 --- !
1114	ph_i(1:idim,1:jpl) = 0._wp
1115	ph_s(1:idim,1:jpl) = 0._wp
1116	pa_i(1:idim,1:jpl) = 0._wp
1117	!
1118	! --- fill the categories --- !
1119	! find where cat-input = cat-output and fill cat-output fields
1120	jlmax(:) = 0
1121	jlmin(:) = 999
1122	jlfil(:,:) = 0
1123	DO jl1 = 1, jpl
1124	DO jl2 = 1, icat
1125	DO ji = 1, idim
1126	IF( hi_max(jl1-1) <= phti(ji,jl2) .AND. hi_max(jl1) > phti(ji,jl2) ) THEN
1127	! fill the right category
1128	ph_i(ji,jl1) = phti(ji,jl2)
1129	ph_s(ji,jl1) = phts(ji,jl2)
1130	pa_i(ji,jl1) = pati(ji,jl2)
1131	! record categories that are filled
1132	jlmax(ji) = MAX( jlmax(ji), jl1 )
1133	jlmin(ji) = MIN( jlmin(ji), jl1 )
1134	jlfil(ji,jl1) = jl1
1135	ENDIF
1136	END DO
1137	END DO
1138	END DO
1139	!
1140	! --- fill the gaps between categories --- !
1141	! transfer from categories filled at the previous step to the empty ones in between
1142	DO ji = 1, idim
1143	jl1 = jlmin(ji)
1144	jl2 = jlmax(ji)
1145	IF( jl1 > 1 ) THEN
1146	! fill the lower cat (jl1-1)
1147	pa_i(ji,jl1-1) = ztrans * pa_i(ji,jl1)
1148	ph_i(ji,jl1-1) = hi_mean(jl1-1)
1149	! remove from cat jl1
1150	pa_i(ji,jl1 ) = ( 1._wp - ztrans ) * pa_i(ji,jl1)
1151	ENDIF
1152	IF( jl2 < jpl ) THEN
1153	! fill the upper cat (jl2+1)
1154	pa_i(ji,jl2+1) = ztrans * pa_i(ji,jl2)
1155	ph_i(ji,jl2+1) = hi_mean(jl2+1)
1156	! remove from cat jl2
1157	pa_i(ji,jl2 ) = ( 1._wp - ztrans ) * pa_i(ji,jl2)
1158	ENDIF
1159	END DO
1160	!
1161	jlfil2(:,:) = jlfil(:,:)
1162	! fill categories from low to high
1163	DO jl = 2, jpl-1
1164	DO ji = 1, idim
1165	IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
1166	! fill high
1167	pa_i(ji,jl) = ztrans * pa_i(ji,jl-1)
1168	ph_i(ji,jl) = hi_mean(jl)
1169	jlfil(ji,jl) = jl
1170	! remove low
1171	pa_i(ji,jl-1) = ( 1._wp - ztrans ) * pa_i(ji,jl-1)
1172	ENDIF
1173	END DO
1174	END DO
1175	!
1176	! fill categories from high to low
1177	DO jl = jpl-1, 2, -1
1178	DO ji = 1, idim
1179	IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
1180	! fill low
1181	pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1)
1182	ph_i(ji,jl) = hi_mean(jl)
1183	jlfil2(ji,jl) = jl
1184	! remove high
1185	pa_i(ji,jl+1) = ( 1._wp - ztrans ) * pa_i(ji,jl+1)
1186	ENDIF
1187	END DO
1188	END DO
1189	!
1190	DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays
1191	DEALLOCATE( jlmin, jlmax )
1192	!
1193	! == temperature and salinity == !
1194	!
1195	ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) )
1196	!
1197	WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp ) ; z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 )
1198	ELSEWHERE ; z1_ai(:) = 0._wp
1199	END WHERE
1200	WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp ) ; z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 )
1201	ELSEWHERE ; z1_vi(:) = 0._wp
1202	END WHERE
1203	WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp ) ; z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 )
1204	ELSEWHERE ; z1_vs(:) = 0._wp
1205	END WHERE
1206	!
1207	! fill all the categories with the same value
1208	ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
1209	DO jl = 1, jpl
1210	pt_i (:,jl) = ztmp(:)
1211	END DO
1212	ztmp(:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
1213	DO jl = 1, jpl
1214	pt_s (:,jl) = ztmp(:)
1215	END DO
1216	ztmp(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:)
1217	DO jl = 1, jpl
1218	pt_su(:,jl) = ztmp(:)
1219	END DO
1220	ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
1221	DO jl = 1, jpl
1222	ps_i (:,jl) = ztmp(:)
1223	END DO
1224	!
1225	DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp )
1226	!
1227	! == ponds == !
1228	ALLOCATE( zfra(idim) )
1229	! keep the same pond fraction atip/ati for each category
1230	WHERE( SUM( pati(:,:), dim=2 ) /= 0._wp ) ; zfra(:) = SUM( patip(:,:), dim=2 ) / SUM( pati(:,:), dim=2 )
1231	ELSEWHERE ; zfra(:) = 0._wp
1232	END WHERE
1233	DO jl = 1, jpl
1234	pa_ip(:,jl) = zfra(:) * pa_i(:,jl)
1235	END DO
1236	! keep the same v_ip/v_i ratio for each category
1237	WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp )
1238	zfra(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 )
1239	ELSEWHERE
1240	zfra(:) = 0._wp
1241	END WHERE
1242	DO jl = 1, jpl
1243	WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1244	ELSEWHERE ; ph_ip(:,jl) = 0._wp
1245	END WHERE
1246	END DO
1247	! keep the same v_il/v_i ratio for each category
1248	WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp )
1249	zfra(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 )
1250	ELSEWHERE
1251	zfra(:) = 0._wp
1252	END WHERE
1253	DO jl = 1, jpl
1254	WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1255	ELSEWHERE ; ph_il(:,jl) = 0._wp
1256	END WHERE
1257	END DO
1258	DEALLOCATE( zfra )
1259	!
1260	ENDIF
1261	!
1262	END SUBROUTINE ice_var_itd_NcMc
1263
1264	!!-------------------------------------------------------------------
1265	!! INTERFACE ice_var_snwfra
1266	!!
1267	!! ** Purpose : fraction of ice covered by snow
1268	!!
1269	!! ** Method : In absence of proper snow model on top of sea ice,
1270	!! we argue that snow does not cover the whole ice because
1271	!! of wind blowing...
1272	!!
1273	!! ** Arguments : ph_s: snow thickness
1274	!!
1275	!! ** Output : pa_s_fra: fraction of ice covered by snow
1276	!!
1277	!!-------------------------------------------------------------------
1278	SUBROUTINE ice_var_snwfra_3d( ph_s, pa_s_fra )
1279	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ph_s ! snow thickness
1280	REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pa_s_fra ! ice fraction covered by snow
1281	IF ( nn_snwfra == 0 ) THEN ! basic 0 or 1 snow cover
1282	WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp
1283	ELSEWHERE ; pa_s_fra = 0._wp
1284	END WHERE
1285	ELSEIF( nn_snwfra == 1 ) THEN ! snow cover depends on hsnow (met-office style)
1286	pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s )
1287	ELSEIF( nn_snwfra == 2 ) THEN ! snow cover depends on hsnow (cice style)
1288	pa_s_fra = ph_s / ( ph_s + 0.02_wp )
1289	ENDIF
1290	END SUBROUTINE ice_var_snwfra_3d
1291
1292	SUBROUTINE ice_var_snwfra_2d( ph_s, pa_s_fra )
1293	REAL(wp), DIMENSION(:,:), INTENT(in ) :: ph_s ! snow thickness
1294	REAL(wp), DIMENSION(:,:), INTENT( out) :: pa_s_fra ! ice fraction covered by snow
1295	IF ( nn_snwfra == 0 ) THEN ! basic 0 or 1 snow cover
1296	WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp
1297	ELSEWHERE ; pa_s_fra = 0._wp
1298	END WHERE
1299	ELSEIF( nn_snwfra == 1 ) THEN ! snow cover depends on hsnow (met-office style)
1300	pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s )
1301	ELSEIF( nn_snwfra == 2 ) THEN ! snow cover depends on hsnow (cice style)
1302	pa_s_fra = ph_s / ( ph_s + 0.02_wp )
1303	ENDIF
1304	END SUBROUTINE ice_var_snwfra_2d
1305
1306	SUBROUTINE ice_var_snwfra_1d( ph_s, pa_s_fra )
1307	REAL(wp), DIMENSION(:), INTENT(in ) :: ph_s ! snow thickness
1308	REAL(wp), DIMENSION(:), INTENT( out) :: pa_s_fra ! ice fraction covered by snow
1309	IF ( nn_snwfra == 0 ) THEN ! basic 0 or 1 snow cover
1310	WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp
1311	ELSEWHERE ; pa_s_fra = 0._wp
1312	END WHERE
1313	ELSEIF( nn_snwfra == 1 ) THEN ! snow cover depends on hsnow (met-office style)
1314	pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s )
1315	ELSEIF( nn_snwfra == 2 ) THEN ! snow cover depends on hsnow (cice style)
1316	pa_s_fra = ph_s / ( ph_s + 0.02_wp )
1317	ENDIF
1318	END SUBROUTINE ice_var_snwfra_1d
1319
1320	!!--------------------------------------------------------------------------
1321	!! INTERFACE ice_var_snwblow
1322	!!
1323	!! ** Purpose : Compute distribution of precip over the ice
1324	!!
1325	!! Snow accumulation in one thermodynamic time step
1326	!! snowfall is partitionned between leads and ice.
1327	!! If snow fall was uniform, a fraction (1-at_i) would fall into leads
1328	!! but because of the winds, more snow falls on leads than on sea ice
1329	!! and a greater fraction (1-at_i)^beta of the total mass of snow
1330	!! (beta < 1) falls in leads.
1331	!! In reality, beta depends on wind speed,
1332	!! and should decrease with increasing wind speed but here, it is
1333	!! considered as a constant. an average value is 0.66
1334	!!--------------------------------------------------------------------------
1335	!!gm I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE....
1336	SUBROUTINE ice_var_snwblow_2d( pin, pout )
1337	REAL(wp), DIMENSION(:,:), INTENT(in ) :: pin ! previous fraction lead ( 1. - a_i_b )
1338	REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
1339	pout = ( 1._wp - ( pin )**rn_snwblow )
1340	END SUBROUTINE ice_var_snwblow_2d
1341
1342	SUBROUTINE ice_var_snwblow_1d( pin, pout )
1343	REAL(wp), DIMENSION(:), INTENT(in ) :: pin
1344	REAL(wp), DIMENSION(:), INTENT(inout) :: pout
1345	pout = ( 1._wp - ( pin )**rn_snwblow )
1346	END SUBROUTINE ice_var_snwblow_1d
1347
1348	#else
1349	!!----------------------------------------------------------------------
1350	!! Default option Dummy module NO SI3 sea-ice model
1351	!!----------------------------------------------------------------------
1352	#endif
1353
1354	!!======================================================================
1355	END MODULE icevar

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

Download in other formats: