Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

icevar.F90

Last change on this file was 14026, checked in by clem, 3 years ago
4.0-HEAD: fix bugs and defects related to tickets #2573 #2575 #2576 #2578. Sette passed and those fixes are now in the trunk, so unless there is a tricky trick somewhere, everything should be fine.
Property svn:keywords set to `Id`
File size: 65.6 KB

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

Download in other formats: