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

Context Navigation

source: NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icevar.F90 @ 10411

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

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

Download in other formats: