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/trunk/src/ICE – NEMO

Context Navigation

source: NEMO/trunk/src/ICE/icevar.F90 @ 10945

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

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

Download in other formats: