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/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icevar.F90 @ 12854

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

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

Download in other formats: