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 @ 12744

Last change on this file since 12744 was 12744, checked in by clem, 4 years ago
make sure all pond lids are set to 0 when not using this option
Property svn:keywords set to `Id`
File size: 60.3 KB

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

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

Download in other formats: