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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

icevar.F90 in NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/ICE – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/ICE/icevar.F90 @ 14451

Last change on this file since 14451 was 14302, checked in by edblockley, 3 years ago
Porting topographic melt-ponds from Martin's branch branches/2020/SI3-05_MeltPonds_topo@14301
File size: 65.7 KB

Line
1	MODULE icevar
2	!!======================================================================
3	!! * MODULE icevar *
4	!! sea-ice: series of functions to transform or compute ice variables
5	!!======================================================================
6	!! History : - ! 2006-01 (M. Vancoppenolle) Original code
7	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
8	!!----------------------------------------------------------------------
9	#if defined key_si3
10	!!----------------------------------------------------------------------
11	!! 'key_si3' SI3 sea-ice model
12	!!----------------------------------------------------------------------
13	!!
14	!! There are three sets of variables
15	!! VGLO : global variables of the model
16	!! - v_i (jpi,jpj,jpl)
17	!! - v_s (jpi,jpj,jpl)
18	!! - a_i (jpi,jpj,jpl)
19	!! - t_s (jpi,jpj,jpl)
20	!! - e_i (jpi,jpj,nlay_i,jpl)
21	!! - e_s (jpi,jpj,nlay_s,jpl)
22	!! - sv_i(jpi,jpj,jpl)
23	!! - oa_i(jpi,jpj,jpl)
24	!! VEQV : equivalent variables sometimes used in the model
25	!! - h_i(jpi,jpj,jpl)
26	!! - h_s(jpi,jpj,jpl)
27	!! - t_i(jpi,jpj,nlay_i,jpl)
28	!! ...
29	!! VAGG : aggregate variables, averaged/summed over all
30	!! thickness categories
31	!! - vt_i(jpi,jpj)
32	!! - vt_s(jpi,jpj)
33	!! - at_i(jpi,jpj)
34	!! - st_i(jpi,jpj)
35	!! - et_s(jpi,jpj) total snow heat content
36	!! - et_i(jpi,jpj) total ice thermal content
37	!! - sm_i(jpi,jpj) mean ice salinity
38	!! - tm_i(jpi,jpj) mean ice temperature
39	!! - tm_s(jpi,jpj) mean snw temperature
40	!!----------------------------------------------------------------------
41	!! ice_var_agg : integrate variables over layers and categories
42	!! ice_var_glo2eqv : transform from VGLO to VEQV
43	!! ice_var_eqv2glo : transform from VEQV to VGLO
44	!! ice_var_salprof : salinity profile in the ice
45	!! ice_var_salprof1d : salinity profile in the ice 1D
46	!! ice_var_zapsmall : remove very small area and volume
47	!! ice_var_zapneg : remove negative ice fields
48	!! ice_var_roundoff : remove negative values arising from roundoff erros
49	!! ice_var_bv : brine volume
50	!! ice_var_enthalpy : compute ice and snow enthalpies from temperature
51	!! ice_var_sshdyn : compute equivalent ssh in lead
52	!! ice_var_itd : convert N-cat to M-cat
53	!! ice_var_snwfra : fraction of ice covered by snow
54	!! 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	& ( h_il(:,:,:) - zhl_min ) / ( 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	!
551	a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
552	v_i (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
553	v_s (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
554	t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
555	oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
556	sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
557	!
558	h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
559	h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
560	!
561	a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
562	v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
563	v_il (ji,jj,jl) = v_il (ji,jj,jl) * zswitch(ji,jj)
564	!
565	END DO
566	END DO
567
568	!-----------------------------------------------------------------
569	! zap small ponds
570	!-----------------------------------------------------------------
571	DO jj = 1, jpj
572	DO ji = 1, jpi
573	IF ( v_ip(ji,jj,jl) <= epsi10 ) THEN
574	a_ip(ji,jj,jl) = 0._wp
575	a_ip_frac(ji,jj,jl) = 0._wp
576	v_ip(ji,jj,jl) = 0._wp
577	h_ip(ji,jj,jl) = 0._wp
578	v_il(ji,jj,jl) = 0._wp
579	h_il(ji,jj,jl) = 0._wp
580	ENDIF
581	END DO
582	END DO
583	!
584	END DO
585
586	! to be sure that at_i is the sum of a_i(jl)
587	at_i (:,:) = SUM( a_i (:,:,:), dim=3 )
588	vt_i (:,:) = SUM( v_i (:,:,:), dim=3 )
589	!!clem add?
590	! vt_s (:,:) = SUM( v_s (:,:,:), dim=3 )
591	! st_i (:,:) = SUM( sv_i(:,:,:), dim=3 )
592	! et_s(:,:) = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 )
593	! et_i(:,:) = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 )
594	!!clem
595
596	! open water = 1 if at_i=0
597	WHERE( at_i(:,:) == 0._wp ) ato_i(:,:) = 1._wp
598	!
599	END SUBROUTINE ice_var_zapsmall
600
601
602	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 )
603	!!-------------------------------------------------------------------
604	!! * ROUTINE ice_var_zapneg *
605	!!
606	!! ** Purpose : Remove negative sea ice fields and correct fluxes
607	!!-------------------------------------------------------------------
608	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
609	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area
610	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume
611	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_s ! snw volume
612	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psv_i ! salt content
613	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content
614	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration
615	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
616	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume
617	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_il ! melt pond lid volume
618	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content
619	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content
620	!
621	INTEGER :: ji, jj, jl, jk ! dummy loop indices
622	REAL(wp) :: z1_dt
623	!!-------------------------------------------------------------------
624	!
625	z1_dt = 1._wp / pdt
626	!
627	DO jl = 1, jpl !== loop over the categories ==!
628	!
629	! make sure a_i=0 where v_i<=0
630	WHERE( pv_i(:,:,:) <= 0._wp ) pa_i(:,:,:) = 0._wp
631
632	!----------------------------------------
633	! zap ice energy and send it to the ocean
634	!----------------------------------------
635	DO jk = 1, nlay_i
636	DO jj = 1 , jpj
637	DO ji = 1 , jpi
638	IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
639	hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
640	pe_i(ji,jj,jk,jl) = 0._wp
641	ENDIF
642	END DO
643	END DO
644	END DO
645	!
646	DO jk = 1, nlay_s
647	DO jj = 1 , jpj
648	DO ji = 1 , jpi
649	IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
650	hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
651	pe_s(ji,jj,jk,jl) = 0._wp
652	ENDIF
653	END DO
654	END DO
655	END DO
656	!
657	!-----------------------------------------------------
658	! zap ice and snow volume, add water and salt to ocean
659	!-----------------------------------------------------
660	DO jj = 1 , jpj
661	DO ji = 1 , jpi
662	IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
663	wfx_res(ji,jj) = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
664	pv_i (ji,jj,jl) = 0._wp
665	ENDIF
666	IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
667	wfx_res(ji,jj) = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt
668	pv_s (ji,jj,jl) = 0._wp
669	ENDIF
670	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
671	sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt
672	psv_i (ji,jj,jl) = 0._wp
673	ENDIF
674	END DO
675	END DO
676	!
677	END DO
678	!
679	WHERE( pato_i(:,:) < 0._wp ) pato_i(:,:) = 0._wp
680	WHERE( poa_i (:,:,:) < 0._wp ) poa_i (:,:,:) = 0._wp
681	WHERE( pa_i (:,:,:) < 0._wp ) pa_i (:,:,:) = 0._wp
682	WHERE( pa_ip (:,:,:) < 0._wp ) pa_ip (:,:,:) = 0._wp
683	WHERE( pv_ip (:,:,:) < 0._wp ) pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
684	WHERE( pv_il (:,:,:) < 0._wp ) pv_il (:,:,:) = 0._wp ! but it does not change conservation, so keep it this way is ok
685	!
686	END SUBROUTINE ice_var_zapneg
687
688
689	SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
690	!!-------------------------------------------------------------------
691	!! * ROUTINE ice_var_roundoff *
692	!!
693	!! ** Purpose : Remove negative sea ice values arising from roundoff errors
694	!!-------------------------------------------------------------------
695	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_i ! ice concentration
696	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_i ! ice volume
697	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_s ! snw volume
698	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: psv_i ! salt content
699	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: poa_i ! age content
700	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
701	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_ip ! melt pond volume
702	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_il ! melt pond lid volume
703	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_s ! snw heat content
704	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_i ! ice heat content
705	!!-------------------------------------------------------------------
706	!
707	WHERE( pa_i (1:npti,:) < 0._wp .AND. pa_i (1:npti,:) > -epsi10 ) pa_i (1:npti,:) = 0._wp ! a_i must be >= 0
708	WHERE( pv_i (1:npti,:) < 0._wp .AND. pv_i (1:npti,:) > -epsi10 ) pv_i (1:npti,:) = 0._wp ! v_i must be >= 0
709	WHERE( pv_s (1:npti,:) < 0._wp .AND. pv_s (1:npti,:) > -epsi10 ) pv_s (1:npti,:) = 0._wp ! v_s must be >= 0
710	WHERE( psv_i(1:npti,:) < 0._wp .AND. psv_i(1:npti,:) > -epsi10 ) psv_i(1:npti,:) = 0._wp ! sv_i must be >= 0
711	WHERE( poa_i(1:npti,:) < 0._wp .AND. poa_i(1:npti,:) > -epsi10 ) poa_i(1:npti,:) = 0._wp ! oa_i must be >= 0
712	WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 ) pe_i (1:npti,:,:) = 0._wp ! e_i must be >= 0
713	WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 ) pe_s (1:npti,:,:) = 0._wp ! e_s must be >= 0
714	IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
715	WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 ) pa_ip(1:npti,:) = 0._wp ! a_ip must be >= 0
716	WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 ) pv_ip(1:npti,:) = 0._wp ! v_ip must be >= 0
717	IF( ln_pnd_lids ) THEN
718	WHERE( pv_il(1:npti,:) < 0._wp .AND. pv_il(1:npti,:) > -epsi10 ) pv_il(1:npti,:) = 0._wp ! v_il must be >= 0
719	ENDIF
720	ENDIF
721	!
722	END SUBROUTINE ice_var_roundoff
723
724
725	SUBROUTINE ice_var_bv
726	!!-------------------------------------------------------------------
727	!! * ROUTINE ice_var_bv *
728	!!
729	!! ** Purpose : computes mean brine volume (%) in sea ice
730	!!
731	!! ** Method : e = - 0.054 * S (ppt) / T (C)
732	!!
733	!! References : Vancoppenolle et al., JGR, 2007
734	!!-------------------------------------------------------------------
735	INTEGER :: ji, jj, jk, jl ! dummy loop indices
736	!!-------------------------------------------------------------------
737	!
738	!!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done
739	!! instead of setting everything to zero as just below
740	bv_i (:,:,:) = 0._wp
741	DO jl = 1, jpl
742	DO jk = 1, nlay_i
743	WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )
744	bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
745	END WHERE
746	END DO
747	END DO
748	WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
749	ELSEWHERE ; bvm_i(:,:) = 0._wp
750	END WHERE
751	!
752	END SUBROUTINE ice_var_bv
753
754
755	SUBROUTINE ice_var_enthalpy
756	!!-------------------------------------------------------------------
757	!! * ROUTINE ice_var_enthalpy *
758	!!
759	!! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature
760	!!
761	!! ** Method : Formula (Bitz and Lipscomb, 1999)
762	!!-------------------------------------------------------------------
763	INTEGER :: ji, jk ! dummy loop indices
764	REAL(wp) :: ztmelts ! local scalar
765	!!-------------------------------------------------------------------
766	!
767	DO jk = 1, nlay_i ! Sea ice energy of melting
768	DO ji = 1, npti
769	ztmelts = - rTmlt * sz_i_1d(ji,jk)
770	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
771	! (sometimes zdf scheme produces abnormally high temperatures)
772	e_i_1d(ji,jk) = rhoi * ( rcpi * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) ) &
773	& + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) ) &
774	& - rcp * ztmelts )
775	END DO
776	END DO
777	DO jk = 1, nlay_s ! Snow energy of melting
778	DO ji = 1, npti
779	e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
780	END DO
781	END DO
782	!
783	END SUBROUTINE ice_var_enthalpy
784
785
786	FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
787	!!---------------------------------------------------------------------
788	!! * ROUTINE ice_var_sshdyn *
789	!!
790	!! ** Purpose : compute the equivalent ssh in lead when sea ice is embedded
791	!!
792	!! ** Method : ssh_lead = ssh + (Mice + Msnow) / rau0
793	!!
794	!! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira,
795	!! Sea ice-ocean coupling using a rescaled vertical coordinate z*,
796	!! Ocean Modelling, Volume 24, Issues 1-2, 2008
797	!!----------------------------------------------------------------------
798	!
799	! input
800	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh !: ssh [m]
801	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass !: mass of snow and ice at current ice time step [Kg/m2]
802	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b !: mass of snow and ice at previous ice time step [Kg/m2]
803	!
804	! output
805	REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn ! equivalent ssh in lead [m]
806	!
807	! temporary
808	REAL(wp) :: zintn, zintb ! time interpolation weights []
809	REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload ! snow and ice load [m]
810	!
811	! compute ice load used to define the equivalent ssh in lead
812	IF( ln_ice_embd ) THEN
813	!
814	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
815	! = (1/nn_fsbc)^2 * {SUM[n] , n=0,nn_fsbc-1}
816	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
817	!
818	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
819	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
820	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
821	!
822	zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0
823	!
824	ELSE
825	zsnwiceload(:,:) = 0.0_wp
826	ENDIF
827	! compute equivalent ssh in lead
828	ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:)
829	!
830	END FUNCTION ice_var_sshdyn
831
832
833	!!-------------------------------------------------------------------
834	!! * INTERFACE ice_var_itd *
835	!!
836	!! ** Purpose : converting N-cat ice to jpl ice categories
837	!!-------------------------------------------------------------------
838	SUBROUTINE ice_var_itd_1c1c( phti, phts, pati , ph_i, ph_s, pa_i, &
839	& ptmi, ptms, ptmsu, psmi, patip, phtip, phtil, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
840	!!-------------------------------------------------------------------
841	!! ** Purpose : converting 1-cat ice to 1 ice category
842	!!-------------------------------------------------------------------
843	REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
844	REAL(wp), DIMENSION(:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
845	REAL(wp), DIMENSION(:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, phtil ! input ice/snow temp & sal & ponds
846	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
847	!!-------------------------------------------------------------------
848	! == thickness and concentration == !
849	ph_i(:) = phti(:)
850	ph_s(:) = phts(:)
851	pa_i(:) = pati(:)
852	!
853	! == temperature and salinity and ponds == !
854	pt_i (:) = ptmi (:)
855	pt_s (:) = ptms (:)
856	pt_su(:) = ptmsu(:)
857	ps_i (:) = psmi (:)
858	pa_ip(:) = patip(:)
859	ph_ip(:) = phtip(:)
860	ph_il(:) = phtil(:)
861
862	END SUBROUTINE ice_var_itd_1c1c
863
864	SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati , ph_i, ph_s, pa_i, &
865	& ptmi, ptms, ptmsu, psmi, patip, phtip, phtil, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
866	!!-------------------------------------------------------------------
867	!! ** Purpose : converting N-cat ice to 1 ice category
868	!!-------------------------------------------------------------------
869	REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
870	REAL(wp), DIMENSION(:) , INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
871	REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, phtil ! input ice/snow temp & sal & ponds
872	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
873	!
874	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs
875	!
876	INTEGER :: idim
877	!!-------------------------------------------------------------------
878	!
879	idim = SIZE( phti, 1 )
880	!
881	! == thickness and concentration == !
882	ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim) )
883	!
884	pa_i(:) = SUM( pati(:,:), dim=2 )
885
886	WHERE( ( pa_i(:) ) /= 0._wp ) ; z1_ai(:) = 1._wp / pa_i(:)
887	ELSEWHERE ; z1_ai(:) = 0._wp
888	END WHERE
889
890	ph_i(:) = SUM( phti(:,:) * pati(:,:), dim=2 ) * z1_ai(:)
891	ph_s(:) = SUM( phts(:,:) * pati(:,:), dim=2 ) * z1_ai(:)
892	!
893	! == temperature and salinity == !
894	WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp ) ; z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) )
895	ELSEWHERE ; z1_vi(:) = 0._wp
896	END WHERE
897	WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp ) ; z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) )
898	ELSEWHERE ; z1_vs(:) = 0._wp
899	END WHERE
900	pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
901	pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
902	pt_su(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:)
903	ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
904
905	! == ponds == !
906	pa_ip(:) = SUM( patip(:,:), dim=2 )
907	WHERE( pa_ip(:) /= 0._wp )
908	ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
909	ph_il(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
910	ELSEWHERE
911	ph_ip(:) = 0._wp
912	ph_il(:) = 0._wp
913	END WHERE
914	!
915	DEALLOCATE( z1_ai, z1_vi, z1_vs )
916	!
917	END SUBROUTINE ice_var_itd_Nc1c
918
919	SUBROUTINE ice_var_itd_1cMc( phti, phts, pati , ph_i, ph_s, pa_i, &
920	& ptmi, ptms, ptmsu, psmi, patip, phtip, phtil, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
921	!!-------------------------------------------------------------------
922	!!
923	!! ** Purpose : converting 1-cat ice to jpl ice categories
924	!!
925	!!
926	!! ** Method: ice thickness distribution follows a gamma function from Abraham et al. (2015)
927	!! it has the property of conserving total concentration and volume
928	!!
929	!!
930	!! ** Arguments : phti: 1-cat ice thickness
931	!! phts: 1-cat snow depth
932	!! pati: 1-cat ice concentration
933	!!
934	!! ** Output : jpl-cat
935	!!
936	!! Abraham, C., Steiner, N., Monahan, A. and Michel, C., 2015.
937	!! Effects of subgrid‐scale snow thickness variability on radiative transfer in sea ice.
938	!! Journal of Geophysical Research: Oceans, 120(8), pp.5597-5614
939	!!-------------------------------------------------------------------
940	REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
941	REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
942	REAL(wp), DIMENSION(:) , INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, phtil ! input ice/snow temp & sal & ponds
943	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
944	!
945	REAL(wp), ALLOCATABLE, DIMENSION(:) :: zfra, z1_hti
946	INTEGER :: ji, jk, jl
947	INTEGER :: idim
948	REAL(wp) :: zv, zdh
949	!!-------------------------------------------------------------------
950	!
951	idim = SIZE( phti , 1 )
952	!
953	ph_i(1:idim,1:jpl) = 0._wp
954	ph_s(1:idim,1:jpl) = 0._wp
955	pa_i(1:idim,1:jpl) = 0._wp
956	!
957	ALLOCATE( z1_hti(idim) )
958	WHERE( phti(:) /= 0._wp ) ; z1_hti(:) = 1._wp / phti(:)
959	ELSEWHERE ; z1_hti(:) = 0._wp
960	END WHERE
961	!
962	! == thickness and concentration == !
963	! for categories 1:jpl-1, integrate the gamma function from hi_max(jl-1) to hi_max(jl)
964	DO jl = 1, jpl-1
965	DO ji = 1, idim
966	!
967	IF( phti(ji) > 0._wp ) THEN
968	! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl)
969	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) ) &
970	& - ( phti(ji) + 2.hi_max(jl ) ) EXP( -2.hi_max(jl )z1_hti(ji) ) )
971	!
972	! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl)
973	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) ) &
974	& * EXP( -2.hi_max(jl-1)z1_hti(ji) ) &
975	& - ( phti(ji)phti(ji) + 2.phti(ji)hi_max(jl) + 2.hi_max(jl)*hi_max(jl) ) &
976	& * EXP(-2.hi_max(jl)z1_hti(ji)) )
977	! thickness
978	IF( pa_i(ji,jl) > epsi06 ) THEN
979	ph_i(ji,jl) = zv / pa_i(ji,jl)
980	ELSE
981	ph_i(ji,jl) = 0.
982	pa_i(ji,jl) = 0.
983	ENDIF
984	ENDIF
985	!
986	ENDDO
987	ENDDO
988	!
989	! for the last category (jpl), integrate the gamma function from hi_max(jpl-1) to infinity
990	DO ji = 1, idim
991	!
992	IF( phti(ji) > 0._wp ) THEN
993	! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jpl-1) to infinity
994	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) )
995
996	! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jpl-1) to infinity
997	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) ) &
998	& * EXP( -2.hi_max(jpl-1)z1_hti(ji) )
999	! thickness
1000	IF( pa_i(ji,jpl) > epsi06 ) THEN
1001	ph_i(ji,jpl) = zv / pa_i(ji,jpl)
1002	else
1003	ph_i(ji,jpl) = 0.
1004	pa_i(ji,jpl) = 0.
1005	ENDIF
1006	ENDIF
1007	!
1008	ENDDO
1009	!
1010	! Add Snow in each category where pa_i is not 0
1011	DO jl = 1, jpl
1012	DO ji = 1, idim
1013	IF( pa_i(ji,jl) > 0._wp ) THEN
1014	ph_s(ji,jl) = ph_i(ji,jl) * phts(ji) * z1_hti(ji)
1015	! In case snow load is in excess that would lead to transformation from snow to ice
1016	! Then, transfer the snow excess into the ice (different from icethd_dh)
1017	zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rau0 ) * ph_i(ji,jl) ) * r1_rau0 )
1018	! recompute h_i, h_s avoiding out of bounds values
1019	ph_i(ji,jl) = MIN( hi_max(jl), ph_i(ji,jl) + zdh )
1020	ph_s(ji,jl) = MAX( 0._wp, ph_s(ji,jl) - zdh * rhoi * r1_rhos )
1021	ENDIF
1022	END DO
1023	END DO
1024	!
1025	DEALLOCATE( z1_hti )
1026	!
1027	! == temperature and salinity == !
1028	DO jl = 1, jpl
1029	pt_i (:,jl) = ptmi (:)
1030	pt_s (:,jl) = ptms (:)
1031	pt_su(:,jl) = ptmsu(:)
1032	ps_i (:,jl) = psmi (:)
1033	END DO
1034	!
1035	! == ponds == !
1036	ALLOCATE( zfra(idim) )
1037	! keep the same pond fraction atip/ati for each category
1038	WHERE( pati(:) /= 0._wp ) ; zfra(:) = patip(:) / pati(:)
1039	ELSEWHERE ; zfra(:) = 0._wp
1040	END WHERE
1041	DO jl = 1, jpl
1042	pa_ip(:,jl) = zfra(:) * pa_i(:,jl)
1043	END DO
1044	! keep the same v_ip/v_i ratio for each category
1045	WHERE( ( phti(:) * pati(:) ) /= 0._wp ) ; zfra(:) = ( phtip(:) * patip(:) ) / ( phti(:) * pati(:) )
1046	ELSEWHERE ; zfra(:) = 0._wp
1047	END WHERE
1048	DO jl = 1, jpl
1049	WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1050	ELSEWHERE ; ph_ip(:,jl) = 0._wp
1051	END WHERE
1052	END DO
1053	! keep the same v_il/v_i ratio for each category
1054	WHERE( ( phti(:) * pati(:) ) /= 0._wp ) ; zfra(:) = ( phtil(:) * patip(:) ) / ( phti(:) * pati(:) )
1055	ELSEWHERE ; zfra(:) = 0._wp
1056	END WHERE
1057	DO jl = 1, jpl
1058	WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1059	ELSEWHERE ; ph_il(:,jl) = 0._wp
1060	END WHERE
1061	END DO
1062	DEALLOCATE( zfra )
1063	!
1064	END SUBROUTINE ice_var_itd_1cMc
1065
1066	SUBROUTINE ice_var_itd_NcMc( phti, phts, pati , ph_i, ph_s, pa_i, &
1067	& ptmi, ptms, ptmsu, psmi, patip, phtip, phtil, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
1068	!!-------------------------------------------------------------------
1069	!!
1070	!! ** Purpose : converting N-cat ice to jpl ice categories
1071	!!
1072	!! ice thickness distribution follows a gaussian law
1073	!! around the concentration of the most likely ice thickness
1074	!! (similar as iceistate.F90)
1075	!!
1076	!! ** Method: Iterative procedure
1077	!!
1078	!! 1) Fill ice cat that correspond to input thicknesses
1079	!! Find the lowest(jlmin) and highest(jlmax) cat that are filled
1080	!!
1081	!! 2) Expand the filling to the cat jlmin-1 and jlmax+1
1082	!! by removing 25% ice area from jlmin and jlmax (resp.)
1083	!!
1084	!! 3) Expand the filling to the empty cat between jlmin and jlmax
1085	!! by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
1086	!! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
1087	!!
1088	!! ** Arguments : phti: N-cat ice thickness
1089	!! phts: N-cat snow depth
1090	!! pati: N-cat ice concentration
1091	!!
1092	!! ** Output : jpl-cat
1093	!!
1094	!! (Example of application: BDY forcings when inputs have N-cat /= jpl)
1095	!!-------------------------------------------------------------------
1096	REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
1097	REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
1098	REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, phtil ! input ice/snow temp & sal & ponds
1099	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
1100	!
1101	INTEGER , ALLOCATABLE, DIMENSION(:,:) :: jlfil, jlfil2
1102	INTEGER , ALLOCATABLE, DIMENSION(:) :: jlmax, jlmin
1103	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs, ztmp, zfra
1104	!
1105	REAL(wp), PARAMETER :: ztrans = 0.25_wp
1106	INTEGER :: ji, jl, jl1, jl2
1107	INTEGER :: idim, icat
1108	!!-------------------------------------------------------------------
1109	!
1110	idim = SIZE( phti, 1 )
1111	icat = SIZE( phti, 2 )
1112	!
1113	! == thickness and concentration == !
1114	! ! ---------------------- !
1115	IF( icat == jpl ) THEN ! input cat = output cat !
1116	! ! ---------------------- !
1117	ph_i(:,:) = phti(:,:)
1118	ph_s(:,:) = phts(:,:)
1119	pa_i(:,:) = pati(:,:)
1120	!
1121	! == temperature and salinity and ponds == !
1122	pt_i (:,:) = ptmi (:,:)
1123	pt_s (:,:) = ptms (:,:)
1124	pt_su(:,:) = ptmsu(:,:)
1125	ps_i (:,:) = psmi (:,:)
1126	pa_ip(:,:) = patip(:,:)
1127	ph_ip(:,:) = phtip(:,:)
1128	ph_il(:,:) = phtil(:,:)
1129	! ! ---------------------- !
1130	ELSEIF( icat == 1 ) THEN ! input cat = 1 !
1131	! ! ---------------------- !
1132	CALL ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), &
1133	& ph_i(:,:), ph_s(:,:), pa_i (:,:), &
1134	& ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), phtil(:,1), &
1135	& pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:), ph_il(:,:) )
1136	! ! ---------------------- !
1137	ELSEIF( jpl == 1 ) THEN ! output cat = 1 !
1138	! ! ---------------------- !
1139	CALL ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), &
1140	& ph_i(:,1), ph_s(:,1), pa_i (:,1), &
1141	& ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), phtil(:,:), &
1142	& pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1), ph_il(:,1) )
1143	! ! ----------------------- !
1144	ELSE ! input cat /= output cat !
1145	! ! ----------------------- !
1146
1147	ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays
1148	ALLOCATE( jlmin(idim), jlmax(idim) )
1149
1150	! --- initialize output fields to 0 --- !
1151	ph_i(1:idim,1:jpl) = 0._wp
1152	ph_s(1:idim,1:jpl) = 0._wp
1153	pa_i(1:idim,1:jpl) = 0._wp
1154	!
1155	! --- fill the categories --- !
1156	! find where cat-input = cat-output and fill cat-output fields
1157	jlmax(:) = 0
1158	jlmin(:) = 999
1159	jlfil(:,:) = 0
1160	DO jl1 = 1, jpl
1161	DO jl2 = 1, icat
1162	DO ji = 1, idim
1163	IF( hi_max(jl1-1) <= phti(ji,jl2) .AND. hi_max(jl1) > phti(ji,jl2) ) THEN
1164	! fill the right category
1165	ph_i(ji,jl1) = phti(ji,jl2)
1166	ph_s(ji,jl1) = phts(ji,jl2)
1167	pa_i(ji,jl1) = pati(ji,jl2)
1168	! record categories that are filled
1169	jlmax(ji) = MAX( jlmax(ji), jl1 )
1170	jlmin(ji) = MIN( jlmin(ji), jl1 )
1171	jlfil(ji,jl1) = jl1
1172	ENDIF
1173	END DO
1174	END DO
1175	END DO
1176	!
1177	! --- fill the gaps between categories --- !
1178	! transfer from categories filled at the previous step to the empty ones in between
1179	DO ji = 1, idim
1180	jl1 = jlmin(ji)
1181	jl2 = jlmax(ji)
1182	IF( jl1 > 1 ) THEN
1183	! fill the lower cat (jl1-1)
1184	pa_i(ji,jl1-1) = ztrans * pa_i(ji,jl1)
1185	ph_i(ji,jl1-1) = hi_mean(jl1-1)
1186	! remove from cat jl1
1187	pa_i(ji,jl1 ) = ( 1._wp - ztrans ) * pa_i(ji,jl1)
1188	ENDIF
1189	IF( jl2 < jpl ) THEN
1190	! fill the upper cat (jl2+1)
1191	pa_i(ji,jl2+1) = ztrans * pa_i(ji,jl2)
1192	ph_i(ji,jl2+1) = hi_mean(jl2+1)
1193	! remove from cat jl2
1194	pa_i(ji,jl2 ) = ( 1._wp - ztrans ) * pa_i(ji,jl2)
1195	ENDIF
1196	END DO
1197	!
1198	jlfil2(:,:) = jlfil(:,:)
1199	! fill categories from low to high
1200	DO jl = 2, jpl-1
1201	DO ji = 1, idim
1202	IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
1203	! fill high
1204	pa_i(ji,jl) = ztrans * pa_i(ji,jl-1)
1205	ph_i(ji,jl) = hi_mean(jl)
1206	jlfil(ji,jl) = jl
1207	! remove low
1208	pa_i(ji,jl-1) = ( 1._wp - ztrans ) * pa_i(ji,jl-1)
1209	ENDIF
1210	END DO
1211	END DO
1212	!
1213	! fill categories from high to low
1214	DO jl = jpl-1, 2, -1
1215	DO ji = 1, idim
1216	IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
1217	! fill low
1218	pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1)
1219	ph_i(ji,jl) = hi_mean(jl)
1220	jlfil2(ji,jl) = jl
1221	! remove high
1222	pa_i(ji,jl+1) = ( 1._wp - ztrans ) * pa_i(ji,jl+1)
1223	ENDIF
1224	END DO
1225	END DO
1226	!
1227	DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays
1228	DEALLOCATE( jlmin, jlmax )
1229	!
1230	! == temperature and salinity == !
1231	!
1232	ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) )
1233	!
1234	WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp ) ; z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 )
1235	ELSEWHERE ; z1_ai(:) = 0._wp
1236	END WHERE
1237	WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp ) ; z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 )
1238	ELSEWHERE ; z1_vi(:) = 0._wp
1239	END WHERE
1240	WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp ) ; z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 )
1241	ELSEWHERE ; z1_vs(:) = 0._wp
1242	END WHERE
1243	!
1244	! fill all the categories with the same value
1245	ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
1246	DO jl = 1, jpl
1247	pt_i (:,jl) = ztmp(:)
1248	END DO
1249	ztmp(:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
1250	DO jl = 1, jpl
1251	pt_s (:,jl) = ztmp(:)
1252	END DO
1253	ztmp(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:)
1254	DO jl = 1, jpl
1255	pt_su(:,jl) = ztmp(:)
1256	END DO
1257	ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
1258	DO jl = 1, jpl
1259	ps_i (:,jl) = ztmp(:)
1260	END DO
1261	!
1262	DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp )
1263	!
1264	! == ponds == !
1265	ALLOCATE( zfra(idim) )
1266	! keep the same pond fraction atip/ati for each category
1267	WHERE( SUM( pati(:,:), dim=2 ) /= 0._wp ) ; zfra(:) = SUM( patip(:,:), dim=2 ) / SUM( pati(:,:), dim=2 )
1268	ELSEWHERE ; zfra(:) = 0._wp
1269	END WHERE
1270	DO jl = 1, jpl
1271	pa_ip(:,jl) = zfra(:) * pa_i(:,jl)
1272	END DO
1273	! keep the same v_ip/v_i ratio for each category
1274	WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp )
1275	zfra(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 )
1276	ELSEWHERE
1277	zfra(:) = 0._wp
1278	END WHERE
1279	DO jl = 1, jpl
1280	WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1281	ELSEWHERE ; ph_ip(:,jl) = 0._wp
1282	END WHERE
1283	END DO
1284	! keep the same v_il/v_i ratio for each category
1285	WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp )
1286	zfra(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 )
1287	ELSEWHERE
1288	zfra(:) = 0._wp
1289	END WHERE
1290	DO jl = 1, jpl
1291	WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1292	ELSEWHERE ; ph_il(:,jl) = 0._wp
1293	END WHERE
1294	END DO
1295	DEALLOCATE( zfra )
1296	!
1297	ENDIF
1298	!
1299	END SUBROUTINE ice_var_itd_NcMc
1300
1301	!!-------------------------------------------------------------------
1302	!! INTERFACE ice_var_snwfra
1303	!!
1304	!! ** Purpose : fraction of ice covered by snow
1305	!!
1306	!! ** Method : In absence of proper snow model on top of sea ice,
1307	!! we argue that snow does not cover the whole ice because
1308	!! of wind blowing...
1309	!!
1310	!! ** Arguments : ph_s: snow thickness
1311	!!
1312	!! ** Output : pa_s_fra: fraction of ice covered by snow
1313	!!
1314	!!-------------------------------------------------------------------
1315	SUBROUTINE ice_var_snwfra_3d( ph_s, pa_s_fra )
1316	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ph_s ! snow thickness
1317	REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pa_s_fra ! ice fraction covered by snow
1318	IF ( nn_snwfra == 0 ) THEN ! basic 0 or 1 snow cover
1319	WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp
1320	ELSEWHERE ; pa_s_fra = 0._wp
1321	END WHERE
1322	ELSEIF( nn_snwfra == 1 ) THEN ! snow cover depends on hsnow (met-office style)
1323	pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s )
1324	ELSEIF( nn_snwfra == 2 ) THEN ! snow cover depends on hsnow (cice style)
1325	pa_s_fra = ph_s / ( ph_s + 0.02_wp )
1326	ENDIF
1327	END SUBROUTINE ice_var_snwfra_3d
1328
1329	SUBROUTINE ice_var_snwfra_2d( ph_s, pa_s_fra )
1330	REAL(wp), DIMENSION(:,:), INTENT(in ) :: ph_s ! snow thickness
1331	REAL(wp), DIMENSION(:,:), INTENT( out) :: pa_s_fra ! ice fraction covered by snow
1332	IF ( nn_snwfra == 0 ) THEN ! basic 0 or 1 snow cover
1333	WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp
1334	ELSEWHERE ; pa_s_fra = 0._wp
1335	END WHERE
1336	ELSEIF( nn_snwfra == 1 ) THEN ! snow cover depends on hsnow (met-office style)
1337	pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s )
1338	ELSEIF( nn_snwfra == 2 ) THEN ! snow cover depends on hsnow (cice style)
1339	pa_s_fra = ph_s / ( ph_s + 0.02_wp )
1340	ENDIF
1341	END SUBROUTINE ice_var_snwfra_2d
1342
1343	SUBROUTINE ice_var_snwfra_1d( ph_s, pa_s_fra )
1344	REAL(wp), DIMENSION(:), INTENT(in ) :: ph_s ! snow thickness
1345	REAL(wp), DIMENSION(:), INTENT( out) :: pa_s_fra ! ice fraction covered by snow
1346	IF ( nn_snwfra == 0 ) THEN ! basic 0 or 1 snow cover
1347	WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp
1348	ELSEWHERE ; pa_s_fra = 0._wp
1349	END WHERE
1350	ELSEIF( nn_snwfra == 1 ) THEN ! snow cover depends on hsnow (met-office style)
1351	pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s )
1352	ELSEIF( nn_snwfra == 2 ) THEN ! snow cover depends on hsnow (cice style)
1353	pa_s_fra = ph_s / ( ph_s + 0.02_wp )
1354	ENDIF
1355	END SUBROUTINE ice_var_snwfra_1d
1356
1357	!!--------------------------------------------------------------------------
1358	!! INTERFACE ice_var_snwblow
1359	!!
1360	!! ** Purpose : Compute distribution of precip over the ice
1361	!!
1362	!! Snow accumulation in one thermodynamic time step
1363	!! snowfall is partitionned between leads and ice.
1364	!! If snow fall was uniform, a fraction (1-at_i) would fall into leads
1365	!! but because of the winds, more snow falls on leads than on sea ice
1366	!! and a greater fraction (1-at_i)^beta of the total mass of snow
1367	!! (beta < 1) falls in leads.
1368	!! In reality, beta depends on wind speed,
1369	!! and should decrease with increasing wind speed but here, it is
1370	!! considered as a constant. an average value is 0.66
1371	!!--------------------------------------------------------------------------
1372	!!gm I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE....
1373	SUBROUTINE ice_var_snwblow_2d( pin, pout )
1374	REAL(wp), DIMENSION(:,:), INTENT(in ) :: pin ! previous fraction lead ( 1. - a_i_b )
1375	REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
1376	pout = ( 1._wp - ( pin )**rn_snwblow )
1377	END SUBROUTINE ice_var_snwblow_2d
1378
1379	SUBROUTINE ice_var_snwblow_1d( pin, pout )
1380	REAL(wp), DIMENSION(:), INTENT(in ) :: pin
1381	REAL(wp), DIMENSION(:), INTENT(inout) :: pout
1382	pout = ( 1._wp - ( pin )**rn_snwblow )
1383	END SUBROUTINE ice_var_snwblow_1d
1384
1385	#else
1386	!!----------------------------------------------------------------------
1387	!! Default option Dummy module NO SI3 sea-ice model
1388	!!----------------------------------------------------------------------
1389	#endif
1390
1391	!!======================================================================
1392	END MODULE icevar

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

Download in other formats: