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/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icevar.F90 @ 11362

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

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

Download in other formats: