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_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icevar.F90 @ 11273

Last change on this file since 11273 was 11223, checked in by smasson, 5 years ago
dev_r10984_HPC-13 : cleaning of rewriting of bdydta, see #2285
Property svn:keywords set to `Id`
File size: 51.7 KB

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

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

Download in other formats: