Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

icevar.F90

Last change on this file was 15385, checked in by clem, 3 years ago
cleaning some ice routines. No change in sette
Property svn:keywords set to `Id`
File size: 64.2 KB

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/trunk/src/ICE/icevar.F90

Download in other formats: