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/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ICE/icevar.F90 @ 13181

Last change on this file since 13181 was 13181, checked in by orioltp, 4 years ago
Several bugfixes added by Sam Hatfield (ECMWF).
Property svn:keywords set to `Id`
File size: 56.3 KB

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

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

Download in other formats: