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

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

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

Download in other formats: