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

Last change on this file since 12725 was 12725, checked in by clem, 4 years ago
fix previous revision. Make sure a_ip and a_ip_eff are correct
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	WHERE( pv_il(1:npti,:) < 0._wp .AND. pv_il(1:npti,:) > -epsi10 ) pv_il(1:npti,:) = 0._wp ! v_il must be >= 0
686	ENDIF
687	!
688	END SUBROUTINE ice_var_roundoff
689
690
691	SUBROUTINE ice_var_bv
692	!!-------------------------------------------------------------------
693	!! * ROUTINE ice_var_bv *
694	!!
695	!! ** Purpose : computes mean brine volume (%) in sea ice
696	!!
697	!! ** Method : e = - 0.054 * S (ppt) / T (C)
698	!!
699	!! References : Vancoppenolle et al., JGR, 2007
700	!!-------------------------------------------------------------------
701	INTEGER :: ji, jj, jk, jl ! dummy loop indices
702	!!-------------------------------------------------------------------
703	!
704	!!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done
705	!! instead of setting everything to zero as just below
706	bv_i (:,:,:) = 0._wp
707	DO jl = 1, jpl
708	DO jk = 1, nlay_i
709	WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )
710	bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
711	END WHERE
712	END DO
713	END DO
714	WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
715	ELSEWHERE ; bvm_i(:,:) = 0._wp
716	END WHERE
717	!
718	END SUBROUTINE ice_var_bv
719
720
721	SUBROUTINE ice_var_enthalpy
722	!!-------------------------------------------------------------------
723	!! * ROUTINE ice_var_enthalpy *
724	!!
725	!! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature
726	!!
727	!! ** Method : Formula (Bitz and Lipscomb, 1999)
728	!!-------------------------------------------------------------------
729	INTEGER :: ji, jk ! dummy loop indices
730	REAL(wp) :: ztmelts ! local scalar
731	!!-------------------------------------------------------------------
732	!
733	DO jk = 1, nlay_i ! Sea ice energy of melting
734	DO ji = 1, npti
735	ztmelts = - rTmlt * sz_i_1d(ji,jk)
736	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
737	! (sometimes zdf scheme produces abnormally high temperatures)
738	e_i_1d(ji,jk) = rhoi * ( rcpi * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) ) &
739	& + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) ) &
740	& - rcp * ztmelts )
741	END DO
742	END DO
743	DO jk = 1, nlay_s ! Snow energy of melting
744	DO ji = 1, npti
745	e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
746	END DO
747	END DO
748	!
749	END SUBROUTINE ice_var_enthalpy
750
751
752	FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
753	!!---------------------------------------------------------------------
754	!! * ROUTINE ice_var_sshdyn *
755	!!
756	!! ** Purpose : compute the equivalent ssh in lead when sea ice is embedded
757	!!
758	!! ** Method : ssh_lead = ssh + (Mice + Msnow) / rau0
759	!!
760	!! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira,
761	!! Sea ice-ocean coupling using a rescaled vertical coordinate z*,
762	!! Ocean Modelling, Volume 24, Issues 1-2, 2008
763	!!----------------------------------------------------------------------
764	!
765	! input
766	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh !: ssh [m]
767	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass !: mass of snow and ice at current ice time step [Kg/m2]
768	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b !: mass of snow and ice at previous ice time step [Kg/m2]
769	!
770	! output
771	REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn ! equivalent ssh in lead [m]
772	!
773	! temporary
774	REAL(wp) :: zintn, zintb ! time interpolation weights []
775	REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload ! snow and ice load [m]
776	!
777	! compute ice load used to define the equivalent ssh in lead
778	IF( ln_ice_embd ) THEN
779	!
780	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
781	! = (1/nn_fsbc)^2 * {SUM[n] , n=0,nn_fsbc-1}
782	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
783	!
784	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
785	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
786	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
787	!
788	zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0
789	!
790	ELSE
791	zsnwiceload(:,:) = 0.0_wp
792	ENDIF
793	! compute equivalent ssh in lead
794	ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:)
795	!
796	END FUNCTION ice_var_sshdyn
797
798
799	!!-------------------------------------------------------------------
800	!! * INTERFACE ice_var_itd *
801	!!
802	!! ** Purpose : converting N-cat ice to jpl ice categories
803	!!-------------------------------------------------------------------
804	SUBROUTINE ice_var_itd_1c1c( phti, phts, pati , ph_i, ph_s, pa_i, &
805	& ptmi, ptms, ptmsu, psmi, patip, phtip, phtil, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
806	!!-------------------------------------------------------------------
807	!! ** Purpose : converting 1-cat ice to 1 ice category
808	!!-------------------------------------------------------------------
809	REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
810	REAL(wp), DIMENSION(:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
811	REAL(wp), DIMENSION(:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, phtil ! input ice/snow temp & sal & ponds
812	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
813	!!-------------------------------------------------------------------
814	! == thickness and concentration == !
815	ph_i(:) = phti(:)
816	ph_s(:) = phts(:)
817	pa_i(:) = pati(:)
818	!
819	! == temperature and salinity and ponds == !
820	pt_i (:) = ptmi (:)
821	pt_s (:) = ptms (:)
822	pt_su(:) = ptmsu(:)
823	ps_i (:) = psmi (:)
824	pa_ip(:) = patip(:)
825	ph_ip(:) = phtip(:)
826	ph_il(:) = phtil(:)
827
828	END SUBROUTINE ice_var_itd_1c1c
829
830	SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati , ph_i, ph_s, pa_i, &
831	& ptmi, ptms, ptmsu, psmi, patip, phtip, phtil, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
832	!!-------------------------------------------------------------------
833	!! ** Purpose : converting N-cat ice to 1 ice category
834	!!-------------------------------------------------------------------
835	REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
836	REAL(wp), DIMENSION(:) , INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
837	REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, phtil ! input ice/snow temp & sal & ponds
838	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
839	!
840	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs
841	!
842	INTEGER :: idim
843	!!-------------------------------------------------------------------
844	!
845	idim = SIZE( phti, 1 )
846	!
847	! == thickness and concentration == !
848	ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim) )
849	!
850	pa_i(:) = SUM( pati(:,:), dim=2 )
851
852	WHERE( ( pa_i(:) ) /= 0._wp ) ; z1_ai(:) = 1._wp / pa_i(:)
853	ELSEWHERE ; z1_ai(:) = 0._wp
854	END WHERE
855
856	ph_i(:) = SUM( phti(:,:) * pati(:,:), dim=2 ) * z1_ai(:)
857	ph_s(:) = SUM( phts(:,:) * pati(:,:), dim=2 ) * z1_ai(:)
858	!
859	! == temperature and salinity == !
860	WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp ) ; z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) )
861	ELSEWHERE ; z1_vi(:) = 0._wp
862	END WHERE
863	WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp ) ; z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) )
864	ELSEWHERE ; z1_vs(:) = 0._wp
865	END WHERE
866	pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
867	pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
868	pt_su(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:)
869	ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
870
871	! == ponds == !
872	pa_ip(:) = SUM( patip(:,:), dim=2 )
873	WHERE( pa_ip(:) /= 0._wp )
874	ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
875	ph_il(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
876	ELSEWHERE
877	ph_ip(:) = 0._wp
878	ph_il(:) = 0._wp
879	END WHERE
880	!
881	DEALLOCATE( z1_ai, z1_vi, z1_vs )
882	!
883	END SUBROUTINE ice_var_itd_Nc1c
884
885	SUBROUTINE ice_var_itd_1cMc( phti, phts, pati , ph_i, ph_s, pa_i, &
886	& ptmi, ptms, ptmsu, psmi, patip, phtip, phtil, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
887	!!-------------------------------------------------------------------
888	!!
889	!! ** Purpose : converting 1-cat ice to jpl ice categories
890	!!
891	!!
892	!! ** Method: ice thickness distribution follows a gamma function from Abraham et al. (2015)
893	!! it has the property of conserving total concentration and volume
894	!!
895	!!
896	!! ** Arguments : phti: 1-cat ice thickness
897	!! phts: 1-cat snow depth
898	!! pati: 1-cat ice concentration
899	!!
900	!! ** Output : jpl-cat
901	!!
902	!! Abraham, C., Steiner, N., Monahan, A. and Michel, C., 2015.
903	!! Effects of subgrid‐scale snow thickness variability on radiative transfer in sea ice.
904	!! Journal of Geophysical Research: Oceans, 120(8), pp.5597-5614
905	!!-------------------------------------------------------------------
906	REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
907	REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
908	REAL(wp), DIMENSION(:) , INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, phtil ! input ice/snow temp & sal & ponds
909	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
910	!
911	REAL(wp), ALLOCATABLE, DIMENSION(:) :: zfra, z1_hti
912	INTEGER :: ji, jk, jl
913	INTEGER :: idim
914	REAL(wp) :: zv, zdh
915	!!-------------------------------------------------------------------
916	!
917	idim = SIZE( phti , 1 )
918	!
919	ph_i(1:idim,1:jpl) = 0._wp
920	ph_s(1:idim,1:jpl) = 0._wp
921	pa_i(1:idim,1:jpl) = 0._wp
922	!
923	ALLOCATE( z1_hti(idim) )
924	WHERE( phti(:) /= 0._wp ) ; z1_hti(:) = 1._wp / phti(:)
925	ELSEWHERE ; z1_hti(:) = 0._wp
926	END WHERE
927	!
928	! == thickness and concentration == !
929	! for categories 1:jpl-1, integrate the gamma function from hi_max(jl-1) to hi_max(jl)
930	DO jl = 1, jpl-1
931	DO ji = 1, idim
932	!
933	IF( phti(ji) > 0._wp ) THEN
934	! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl)
935	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) ) &
936	& - ( phti(ji) + 2.hi_max(jl ) ) EXP( -2.hi_max(jl )z1_hti(ji) ) )
937	!
938	! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl)
939	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) ) &
940	& * EXP( -2.hi_max(jl-1)z1_hti(ji) ) &
941	& - ( phti(ji)phti(ji) + 2.phti(ji)hi_max(jl) + 2.hi_max(jl)*hi_max(jl) ) &
942	& * EXP(-2.hi_max(jl)z1_hti(ji)) )
943	! thickness
944	IF( pa_i(ji,jl) > epsi06 ) THEN
945	ph_i(ji,jl) = zv / pa_i(ji,jl)
946	ELSE
947	ph_i(ji,jl) = 0.
948	pa_i(ji,jl) = 0.
949	ENDIF
950	ENDIF
951	!
952	ENDDO
953	ENDDO
954	!
955	! for the last category (jpl), integrate the gamma function from hi_max(jpl-1) to infinity
956	DO ji = 1, idim
957	!
958	IF( phti(ji) > 0._wp ) THEN
959	! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jpl-1) to infinity
960	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) )
961
962	! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jpl-1) to infinity
963	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) ) &
964	& * EXP( -2.hi_max(jpl-1)z1_hti(ji) )
965	! thickness
966	IF( pa_i(ji,jpl) > epsi06 ) THEN
967	ph_i(ji,jpl) = zv / pa_i(ji,jpl)
968	else
969	ph_i(ji,jpl) = 0.
970	pa_i(ji,jpl) = 0.
971	ENDIF
972	ENDIF
973	!
974	ENDDO
975	!
976	! Add Snow in each category where pa_i is not 0
977	DO jl = 1, jpl
978	DO ji = 1, idim
979	IF( pa_i(ji,jl) > 0._wp ) THEN
980	ph_s(ji,jl) = ph_i(ji,jl) * phts(ji) * z1_hti(ji)
981	! In case snow load is in excess that would lead to transformation from snow to ice
982	! Then, transfer the snow excess into the ice (different from icethd_dh)
983	zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rau0 ) * ph_i(ji,jl) ) * r1_rau0 )
984	! recompute h_i, h_s avoiding out of bounds values
985	ph_i(ji,jl) = MIN( hi_max(jl), ph_i(ji,jl) + zdh )
986	ph_s(ji,jl) = MAX( 0._wp, ph_s(ji,jl) - zdh * rhoi * r1_rhos )
987	ENDIF
988	END DO
989	END DO
990	!
991	DEALLOCATE( z1_hti )
992	!
993	! == temperature and salinity == !
994	DO jl = 1, jpl
995	pt_i (:,jl) = ptmi (:)
996	pt_s (:,jl) = ptms (:)
997	pt_su(:,jl) = ptmsu(:)
998	ps_i (:,jl) = psmi (:)
999	ps_i (:,jl) = psmi (:)
1000	END DO
1001	!
1002	! == ponds == !
1003	ALLOCATE( zfra(idim) )
1004	! keep the same pond fraction atip/ati for each category
1005	WHERE( pati(:) /= 0._wp ) ; zfra(:) = patip(:) / pati(:)
1006	ELSEWHERE ; zfra(:) = 0._wp
1007	END WHERE
1008	DO jl = 1, jpl
1009	pa_ip(:,jl) = zfra(:) * pa_i(:,jl)
1010	END DO
1011	! keep the same v_ip/v_i ratio for each category
1012	WHERE( ( phti(:) * pati(:) ) /= 0._wp ) ; zfra(:) = ( phtip(:) * patip(:) ) / ( phti(:) * pati(:) )
1013	ELSEWHERE ; zfra(:) = 0._wp
1014	END WHERE
1015	DO jl = 1, jpl
1016	WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1017	ELSEWHERE ; ph_ip(:,jl) = 0._wp
1018	END WHERE
1019	END DO
1020	! keep the same v_il/v_i ratio for each category
1021	WHERE( ( phti(:) * pati(:) ) /= 0._wp ) ; zfra(:) = ( phtil(:) * patip(:) ) / ( phti(:) * pati(:) )
1022	ELSEWHERE ; zfra(:) = 0._wp
1023	END WHERE
1024	DO jl = 1, jpl
1025	WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1026	ELSEWHERE ; ph_il(:,jl) = 0._wp
1027	END WHERE
1028	END DO
1029	DEALLOCATE( zfra )
1030	!
1031	END SUBROUTINE ice_var_itd_1cMc
1032
1033	SUBROUTINE ice_var_itd_NcMc( phti, phts, pati , ph_i, ph_s, pa_i, &
1034	& ptmi, ptms, ptmsu, psmi, patip, phtip, phtil, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
1035	!!-------------------------------------------------------------------
1036	!!
1037	!! ** Purpose : converting N-cat ice to jpl ice categories
1038	!!
1039	!! ice thickness distribution follows a gaussian law
1040	!! around the concentration of the most likely ice thickness
1041	!! (similar as iceistate.F90)
1042	!!
1043	!! ** Method: Iterative procedure
1044	!!
1045	!! 1) Fill ice cat that correspond to input thicknesses
1046	!! Find the lowest(jlmin) and highest(jlmax) cat that are filled
1047	!!
1048	!! 2) Expand the filling to the cat jlmin-1 and jlmax+1
1049	!! by removing 25% ice area from jlmin and jlmax (resp.)
1050	!!
1051	!! 3) Expand the filling to the empty cat between jlmin and jlmax
1052	!! by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
1053	!! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
1054	!!
1055	!! ** Arguments : phti: N-cat ice thickness
1056	!! phts: N-cat snow depth
1057	!! pati: N-cat ice concentration
1058	!!
1059	!! ** Output : jpl-cat
1060	!!
1061	!! (Example of application: BDY forcings when inputs have N-cat /= jpl)
1062	!!-------------------------------------------------------------------
1063	REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
1064	REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
1065	REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, phtil ! input ice/snow temp & sal & ponds
1066	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
1067	!
1068	INTEGER , ALLOCATABLE, DIMENSION(:,:) :: jlfil, jlfil2
1069	INTEGER , ALLOCATABLE, DIMENSION(:) :: jlmax, jlmin
1070	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs, ztmp, zfra
1071	!
1072	REAL(wp), PARAMETER :: ztrans = 0.25_wp
1073	INTEGER :: ji, jl, jl1, jl2
1074	INTEGER :: idim, icat
1075	!!-------------------------------------------------------------------
1076	!
1077	idim = SIZE( phti, 1 )
1078	icat = SIZE( phti, 2 )
1079	!
1080	! == thickness and concentration == !
1081	! ! ---------------------- !
1082	IF( icat == jpl ) THEN ! input cat = output cat !
1083	! ! ---------------------- !
1084	ph_i(:,:) = phti(:,:)
1085	ph_s(:,:) = phts(:,:)
1086	pa_i(:,:) = pati(:,:)
1087	!
1088	! == temperature and salinity and ponds == !
1089	pt_i (:,:) = ptmi (:,:)
1090	pt_s (:,:) = ptms (:,:)
1091	pt_su(:,:) = ptmsu(:,:)
1092	ps_i (:,:) = psmi (:,:)
1093	pa_ip(:,:) = patip(:,:)
1094	ph_ip(:,:) = phtip(:,:)
1095	ph_il(:,:) = phtil(:,:)
1096	! ! ---------------------- !
1097	ELSEIF( icat == 1 ) THEN ! input cat = 1 !
1098	! ! ---------------------- !
1099	CALL ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), &
1100	& ph_i(:,:), ph_s(:,:), pa_i (:,:), &
1101	& ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), phtil(:,1), &
1102	& pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:), ph_il(:,:) )
1103	! ! ---------------------- !
1104	ELSEIF( jpl == 1 ) THEN ! output cat = 1 !
1105	! ! ---------------------- !
1106	CALL ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), &
1107	& ph_i(:,1), ph_s(:,1), pa_i (:,1), &
1108	& ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), phtil(:,:), &
1109	& pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1), ph_il(:,1) )
1110	! ! ----------------------- !
1111	ELSE ! input cat /= output cat !
1112	! ! ----------------------- !
1113
1114	ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays
1115	ALLOCATE( jlmin(idim), jlmax(idim) )
1116
1117	! --- initialize output fields to 0 --- !
1118	ph_i(1:idim,1:jpl) = 0._wp
1119	ph_s(1:idim,1:jpl) = 0._wp
1120	pa_i(1:idim,1:jpl) = 0._wp
1121	!
1122	! --- fill the categories --- !
1123	! find where cat-input = cat-output and fill cat-output fields
1124	jlmax(:) = 0
1125	jlmin(:) = 999
1126	jlfil(:,:) = 0
1127	DO jl1 = 1, jpl
1128	DO jl2 = 1, icat
1129	DO ji = 1, idim
1130	IF( hi_max(jl1-1) <= phti(ji,jl2) .AND. hi_max(jl1) > phti(ji,jl2) ) THEN
1131	! fill the right category
1132	ph_i(ji,jl1) = phti(ji,jl2)
1133	ph_s(ji,jl1) = phts(ji,jl2)
1134	pa_i(ji,jl1) = pati(ji,jl2)
1135	! record categories that are filled
1136	jlmax(ji) = MAX( jlmax(ji), jl1 )
1137	jlmin(ji) = MIN( jlmin(ji), jl1 )
1138	jlfil(ji,jl1) = jl1
1139	ENDIF
1140	END DO
1141	END DO
1142	END DO
1143	!
1144	! --- fill the gaps between categories --- !
1145	! transfer from categories filled at the previous step to the empty ones in between
1146	DO ji = 1, idim
1147	jl1 = jlmin(ji)
1148	jl2 = jlmax(ji)
1149	IF( jl1 > 1 ) THEN
1150	! fill the lower cat (jl1-1)
1151	pa_i(ji,jl1-1) = ztrans * pa_i(ji,jl1)
1152	ph_i(ji,jl1-1) = hi_mean(jl1-1)
1153	! remove from cat jl1
1154	pa_i(ji,jl1 ) = ( 1._wp - ztrans ) * pa_i(ji,jl1)
1155	ENDIF
1156	IF( jl2 < jpl ) THEN
1157	! fill the upper cat (jl2+1)
1158	pa_i(ji,jl2+1) = ztrans * pa_i(ji,jl2)
1159	ph_i(ji,jl2+1) = hi_mean(jl2+1)
1160	! remove from cat jl2
1161	pa_i(ji,jl2 ) = ( 1._wp - ztrans ) * pa_i(ji,jl2)
1162	ENDIF
1163	END DO
1164	!
1165	jlfil2(:,:) = jlfil(:,:)
1166	! fill categories from low to high
1167	DO jl = 2, jpl-1
1168	DO ji = 1, idim
1169	IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
1170	! fill high
1171	pa_i(ji,jl) = ztrans * pa_i(ji,jl-1)
1172	ph_i(ji,jl) = hi_mean(jl)
1173	jlfil(ji,jl) = jl
1174	! remove low
1175	pa_i(ji,jl-1) = ( 1._wp - ztrans ) * pa_i(ji,jl-1)
1176	ENDIF
1177	END DO
1178	END DO
1179	!
1180	! fill categories from high to low
1181	DO jl = jpl-1, 2, -1
1182	DO ji = 1, idim
1183	IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
1184	! fill low
1185	pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1)
1186	ph_i(ji,jl) = hi_mean(jl)
1187	jlfil2(ji,jl) = jl
1188	! remove high
1189	pa_i(ji,jl+1) = ( 1._wp - ztrans ) * pa_i(ji,jl+1)
1190	ENDIF
1191	END DO
1192	END DO
1193	!
1194	DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays
1195	DEALLOCATE( jlmin, jlmax )
1196	!
1197	! == temperature and salinity == !
1198	!
1199	ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) )
1200	!
1201	WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp ) ; z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 )
1202	ELSEWHERE ; z1_ai(:) = 0._wp
1203	END WHERE
1204	WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp ) ; z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 )
1205	ELSEWHERE ; z1_vi(:) = 0._wp
1206	END WHERE
1207	WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp ) ; z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 )
1208	ELSEWHERE ; z1_vs(:) = 0._wp
1209	END WHERE
1210	!
1211	! fill all the categories with the same value
1212	ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
1213	DO jl = 1, jpl
1214	pt_i (:,jl) = ztmp(:)
1215	END DO
1216	ztmp(:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
1217	DO jl = 1, jpl
1218	pt_s (:,jl) = ztmp(:)
1219	END DO
1220	ztmp(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:)
1221	DO jl = 1, jpl
1222	pt_su(:,jl) = ztmp(:)
1223	END DO
1224	ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
1225	DO jl = 1, jpl
1226	ps_i (:,jl) = ztmp(:)
1227	END DO
1228	!
1229	DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp )
1230	!
1231	! == ponds == !
1232	ALLOCATE( zfra(idim) )
1233	! keep the same pond fraction atip/ati for each category
1234	WHERE( SUM( pati(:,:), dim=2 ) /= 0._wp ) ; zfra(:) = SUM( patip(:,:), dim=2 ) / SUM( pati(:,:), dim=2 )
1235	ELSEWHERE ; zfra(:) = 0._wp
1236	END WHERE
1237	DO jl = 1, jpl
1238	pa_ip(:,jl) = zfra(:) * pa_i(:,jl)
1239	END DO
1240	! keep the same v_ip/v_i ratio for each category
1241	WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp )
1242	zfra(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 )
1243	ELSEWHERE
1244	zfra(:) = 0._wp
1245	END WHERE
1246	DO jl = 1, jpl
1247	WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1248	ELSEWHERE ; ph_ip(:,jl) = 0._wp
1249	END WHERE
1250	END DO
1251	! keep the same v_il/v_i ratio for each category
1252	WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp )
1253	zfra(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 )
1254	ELSEWHERE
1255	zfra(:) = 0._wp
1256	END WHERE
1257	DO jl = 1, jpl
1258	WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1259	ELSEWHERE ; ph_il(:,jl) = 0._wp
1260	END WHERE
1261	END DO
1262	DEALLOCATE( zfra )
1263	!
1264	ENDIF
1265	!
1266	END SUBROUTINE ice_var_itd_NcMc
1267
1268	#else
1269	!!----------------------------------------------------------------------
1270	!! Default option Dummy module NO SI3 sea-ice model
1271	!!----------------------------------------------------------------------
1272	#endif
1273
1274	!!======================================================================
1275	END MODULE icevar

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

Download in other formats: