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/temporary_r4_trunk/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/temporary_r4_trunk/src/ICE/icevar.F90 @ 13470

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

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

Download in other formats: