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/releases/release-4.0/src/ICE – NEMO

Context Navigation

source: NEMO/releases/release-4.0/src/ICE/icevar.F90 @ 10929

Last change on this file since 10929 was 10929, checked in by clem, 5 years ago
solve a problem of conservation when ice advection is called twice (because of a cfl that is too large)
Property svn:keywords set to `Id`
File size: 45.8 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	!----------------------------------------
566	! zap ice energy and send it to the ocean
567	!----------------------------------------
568	DO jk = 1, nlay_i
569	DO jj = 1 , jpj
570	DO ji = 1 , jpi
571	IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
572	hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
573	pe_i(ji,jj,jk,jl) = 0._wp
574	ENDIF
575	END DO
576	END DO
577	END DO
578	!
579	DO jk = 1, nlay_s
580	DO jj = 1 , jpj
581	DO ji = 1 , jpi
582	IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
583	hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
584	pe_s(ji,jj,jk,jl) = 0._wp
585	ENDIF
586	END DO
587	END DO
588	END DO
589	!
590	!-----------------------------------------------------
591	! zap ice and snow volume, add water and salt to ocean
592	!-----------------------------------------------------
593	DO jj = 1 , jpj
594	DO ji = 1 , jpi
595	IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
596	wfx_res(ji,jj) = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
597	pv_i (ji,jj,jl) = 0._wp
598	ENDIF
599	IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
600	wfx_res(ji,jj) = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt
601	pv_s (ji,jj,jl) = 0._wp
602	ENDIF
603	IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN
604	sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt
605	psv_i (ji,jj,jl) = 0._wp
606	ENDIF
607	END DO
608	END DO
609	!
610	END DO
611	!
612	WHERE( pato_i(:,:) < 0._wp ) pato_i(:,:) = 0._wp
613	WHERE( poa_i (:,:,:) < 0._wp ) poa_i (:,:,:) = 0._wp
614	WHERE( pa_i (:,:,:) < 0._wp ) pa_i (:,:,:) = 0._wp
615	WHERE( pa_ip (:,:,:) < 0._wp ) pa_ip (:,:,:) = 0._wp
616	WHERE( pv_ip (:,:,:) < 0._wp ) pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
617	! but it does not change conservation, so keep it this way is ok
618	!
619	END SUBROUTINE ice_var_zapneg
620
621
622	SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
623	!!-------------------------------------------------------------------
624	!! * ROUTINE ice_var_itd *
625	!!
626	!! ** Purpose : converting 1-cat ice to multiple ice categories
627	!!
628	!! ice thickness distribution follows a gaussian law
629	!! around the concentration of the most likely ice thickness
630	!! (similar as iceistate.F90)
631	!!
632	!! ** Method: Iterative procedure
633	!!
634	!! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
635	!!
636	!! 2) Check whether the distribution conserves area and volume, positivity and
637	!! category boundaries
638	!!
639	!! 3) If not (input ice is too thin), the last category is empty and
640	!! the number of categories is reduced (jpl-1)
641	!!
642	!! 4) Iterate until ok (SUM(itest(:) = 4)
643	!!
644	!! ** Arguments : zhti: 1-cat ice thickness
645	!! zhts: 1-cat snow depth
646	!! zati: 1-cat ice concentration
647	!!
648	!! ** Output : jpl-cat
649	!!
650	!! (Example of application: BDY forcings when input are cell averaged)
651	!!-------------------------------------------------------------------
652	INTEGER :: ji, jk, jl ! dummy loop indices
653	INTEGER :: idim, i_fill, jl0
654	REAL(wp) :: zarg, zV, zconv, zdh, zdv
655	REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables
656	REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables
657	INTEGER , DIMENSION(4) :: itest
658	!!-------------------------------------------------------------------
659	!
660	! ----------------------------------------
661	! distribution over the jpl ice categories
662	! ----------------------------------------
663	! a gaussian distribution for ice concentration is used
664	! then we check whether the distribution fullfills
665	! volume and area conservation, positivity and ice categories bounds
666	idim = SIZE( zhti , 1 )
667	zh_i(1:idim,1:jpl) = 0._wp
668	zh_s(1:idim,1:jpl) = 0._wp
669	za_i(1:idim,1:jpl) = 0._wp
670	!
671	DO ji = 1, idim
672	!
673	IF( zhti(ji) > 0._wp ) THEN
674	!
675	! find which category (jl0) the input ice thickness falls into
676	jl0 = jpl
677	DO jl = 1, jpl
678	IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
679	jl0 = jl
680	CYCLE
681	ENDIF
682	END DO
683	!
684	itest(:) = 0
685	i_fill = jpl + 1 !------------------------------------
686	DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories
687	! !------------------------------------
688	i_fill = i_fill - 1
689	!
690	zh_i(ji,1:jpl) = 0._wp
691	za_i(ji,1:jpl) = 0._wp
692	itest(:) = 0
693	!
694	IF ( i_fill == 1 ) THEN !-- case very thin ice: fill only category 1
695	zh_i(ji,1) = zhti(ji)
696	za_i (ji,1) = zati (ji)
697	ELSE !-- case ice is thicker: fill categories >1
698	! thickness
699	DO jl = 1, i_fill - 1
700	zh_i(ji,jl) = hi_mean(jl)
701	END DO
702	!
703	! concentration
704	za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
705	DO jl = 1, i_fill - 1
706	IF ( jl /= jl0 ) THEN
707	zarg = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
708	za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2)
709	ENDIF
710	END DO
711	!
712	! last category
713	za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
714	zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
715	zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )
716	!
717	! correction if concentration of upper cat is greater than lower cat
718	! (it should be a gaussian around jl0 but sometimes it is not)
719	IF ( jl0 /= jpl ) THEN
720	DO jl = jpl, jl0+1, -1
721	IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
722	zdv = zh_i(ji,jl) * za_i(ji,jl)
723	zh_i(ji,jl ) = 0._wp
724	za_i (ji,jl ) = 0._wp
725	za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
726	END IF
727	END DO
728	ENDIF
729	!
730	ENDIF
731	!
732	! Compatibility tests
733	zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) )
734	IF ( zconv < epsi06 ) itest(1) = 1 ! Test 1: area conservation
735	!
736	zconv = ABS( zhti(ji)zati(ji) - SUM( za_i(ji,1:jpl)zh_i(ji,1:jpl) ) )
737	IF ( zconv < epsi06 ) itest(2) = 1 ! Test 2: volume conservation
738	!
739	IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ?
740	!
741	itest(4) = 1
742	DO jl = 1, i_fill
743	IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations
744	END DO
745	! !----------------------------
746	END DO ! end iteration on categories
747	! !----------------------------
748	ENDIF
749	END DO
750
751	! Add Snow in each category where za_i is not 0
752	DO jl = 1, jpl
753	DO ji = 1, idim
754	IF( za_i(ji,jl) > 0._wp ) THEN
755	zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
756	! In case snow load is in excess that would lead to transformation from snow to ice
757	! Then, transfer the snow excess into the ice (different from icethd_dh)
758	zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 )
759	! recompute h_i, h_s avoiding out of bounds values
760	zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
761	zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos )
762	ENDIF
763	END DO
764	END DO
765	!
766	END SUBROUTINE ice_var_itd
767
768
769	SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )
770	!!-------------------------------------------------------------------
771	!! * ROUTINE ice_var_itd2 *
772	!!
773	!! ** Purpose : converting N-cat ice to jpl ice categories
774	!!
775	!! ice thickness distribution follows a gaussian law
776	!! around the concentration of the most likely ice thickness
777	!! (similar as iceistate.F90)
778	!!
779	!! ** Method: Iterative procedure
780	!!
781	!! 1) Fill ice cat that correspond to input thicknesses
782	!! Find the lowest(jlmin) and highest(jlmax) cat that are filled
783	!!
784	!! 2) Expand the filling to the cat jlmin-1 and jlmax+1
785	!! by removing 25% ice area from jlmin and jlmax (resp.)
786	!!
787	!! 3) Expand the filling to the empty cat between jlmin and jlmax
788	!! by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
789	!! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
790	!!
791	!! ** Arguments : zhti: N-cat ice thickness
792	!! zhts: N-cat snow depth
793	!! zati: N-cat ice concentration
794	!!
795	!! ** Output : jpl-cat
796	!!
797	!! (Example of application: BDY forcings when inputs have N-cat /= jpl)
798	!!-------------------------------------------------------------------
799	INTEGER :: ji, jl, jl1, jl2 ! dummy loop indices
800	INTEGER :: idim, icat
801	REAL(wp), PARAMETER :: ztrans = 0.25_wp
802	REAL(wp), DIMENSION(:,:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables
803	REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables
804	INTEGER , DIMENSION(:,:), ALLOCATABLE :: jlfil, jlfil2
805	INTEGER , DIMENSION(:) , ALLOCATABLE :: jlmax, jlmin
806	!!-------------------------------------------------------------------
807	!
808	idim = SIZE( zhti, 1 )
809	icat = SIZE( zhti, 2 )
810	!
811	ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays
812	ALLOCATE( jlmin(idim), jlmax(idim) )
813
814	! --- initialize output fields to 0 --- !
815	zh_i(1:idim,1:jpl) = 0._wp
816	zh_s(1:idim,1:jpl) = 0._wp
817	za_i(1:idim,1:jpl) = 0._wp
818	!
819	! --- fill the categories --- !
820	! find where cat-input = cat-output and fill cat-output fields
821	jlmax(:) = 0
822	jlmin(:) = 999
823	jlfil(:,:) = 0
824	DO jl1 = 1, jpl
825	DO jl2 = 1, icat
826	DO ji = 1, idim
827	IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
828	! fill the right category
829	zh_i(ji,jl1) = zhti(ji,jl2)
830	zh_s(ji,jl1) = zhts(ji,jl2)
831	za_i(ji,jl1) = zati(ji,jl2)
832	! record categories that are filled
833	jlmax(ji) = MAX( jlmax(ji), jl1 )
834	jlmin(ji) = MIN( jlmin(ji), jl1 )
835	jlfil(ji,jl1) = jl1
836	ENDIF
837	END DO
838	END DO
839	END DO
840	!
841	! --- fill the gaps between categories --- !
842	! transfer from categories filled at the previous step to the empty ones in between
843	DO ji = 1, idim
844	jl1 = jlmin(ji)
845	jl2 = jlmax(ji)
846	IF( jl1 > 1 ) THEN
847	! fill the lower cat (jl1-1)
848	za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
849	zh_i(ji,jl1-1) = hi_mean(jl1-1)
850	! remove from cat jl1
851	za_i(ji,jl1 ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
852	ENDIF
853	IF( jl2 < jpl ) THEN
854	! fill the upper cat (jl2+1)
855	za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
856	zh_i(ji,jl2+1) = hi_mean(jl2+1)
857	! remove from cat jl2
858	za_i(ji,jl2 ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
859	ENDIF
860	END DO
861	!
862	jlfil2(:,:) = jlfil(:,:)
863	! fill categories from low to high
864	DO jl = 2, jpl-1
865	DO ji = 1, idim
866	IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
867	! fill high
868	za_i(ji,jl) = ztrans * za_i(ji,jl-1)
869	zh_i(ji,jl) = hi_mean(jl)
870	jlfil(ji,jl) = jl
871	! remove low
872	za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
873	ENDIF
874	END DO
875	END DO
876	!
877	! fill categories from high to low
878	DO jl = jpl-1, 2, -1
879	DO ji = 1, idim
880	IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
881	! fill low
882	za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
883	zh_i(ji,jl) = hi_mean(jl)
884	jlfil2(ji,jl) = jl
885	! remove high
886	za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
887	ENDIF
888	END DO
889	END DO
890	!
891	DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays
892	DEALLOCATE( jlmin, jlmax )
893	!
894	END SUBROUTINE ice_var_itd2
895
896
897	SUBROUTINE ice_var_bv
898	!!-------------------------------------------------------------------
899	!! * ROUTINE ice_var_bv *
900	!!
901	!! ** Purpose : computes mean brine volume (%) in sea ice
902	!!
903	!! ** Method : e = - 0.054 * S (ppt) / T (C)
904	!!
905	!! References : Vancoppenolle et al., JGR, 2007
906	!!-------------------------------------------------------------------
907	INTEGER :: ji, jj, jk, jl ! dummy loop indices
908	!!-------------------------------------------------------------------
909	!
910	!!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done
911	!! instead of setting everything to zero as just below
912	bv_i (:,:,:) = 0._wp
913	DO jl = 1, jpl
914	DO jk = 1, nlay_i
915	WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )
916	bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
917	END WHERE
918	END DO
919	END DO
920	WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
921	ELSEWHERE ; bvm_i(:,:) = 0._wp
922	END WHERE
923	!
924	END SUBROUTINE ice_var_bv
925
926
927	SUBROUTINE ice_var_enthalpy
928	!!-------------------------------------------------------------------
929	!! * ROUTINE ice_var_enthalpy *
930	!!
931	!! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature
932	!!
933	!! ** Method : Formula (Bitz and Lipscomb, 1999)
934	!!-------------------------------------------------------------------
935	INTEGER :: ji, jk ! dummy loop indices
936	REAL(wp) :: ztmelts ! local scalar
937	!!-------------------------------------------------------------------
938	!
939	DO jk = 1, nlay_i ! Sea ice energy of melting
940	DO ji = 1, npti
941	ztmelts = - rTmlt * sz_i_1d(ji,jk)
942	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
943	! (sometimes zdf scheme produces abnormally high temperatures)
944	e_i_1d(ji,jk) = rhoi * ( rcpi * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) ) &
945	& + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) ) &
946	& - rcp * ztmelts )
947	END DO
948	END DO
949	DO jk = 1, nlay_s ! Snow energy of melting
950	DO ji = 1, npti
951	e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
952	END DO
953	END DO
954	!
955	END SUBROUTINE ice_var_enthalpy
956
957	FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
958	!!---------------------------------------------------------------------
959	!! * ROUTINE ice_var_sshdyn *
960	!!
961	!! ** Purpose : compute the equivalent ssh in lead when sea ice is embedded
962	!!
963	!! ** Method : ssh_lead = ssh + (Mice + Msnow) / rau0
964	!!
965	!! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira,
966	!! Sea ice-ocean coupling using a rescaled vertical coordinate z*,
967	!! Ocean Modelling, Volume 24, Issues 1-2, 2008
968	!!----------------------------------------------------------------------
969	!
970	! input
971	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh !: ssh [m]
972	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass !: mass of snow and ice at current ice time step [Kg/m2]
973	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b !: mass of snow and ice at previous ice time step [Kg/m2]
974	!
975	! output
976	REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn ! equivalent ssh in lead [m]
977	!
978	! temporary
979	REAL(wp) :: zintn, zintb ! time interpolation weights []
980	REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload ! snow and ice load [m]
981	!
982	! compute ice load used to define the equivalent ssh in lead
983	IF( ln_ice_embd ) THEN
984	!
985	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
986	! = (1/nn_fsbc)^2 * {SUM[n] , n=0,nn_fsbc-1}
987	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
988	!
989	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
990	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
991	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
992	!
993	zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0
994	!
995	ELSE
996	zsnwiceload(:,:) = 0.0_wp
997	ENDIF
998	! compute equivalent ssh in lead
999	ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:)
1000	!
1001	END FUNCTION ice_var_sshdyn
1002
1003
1004	#else
1005	!!----------------------------------------------------------------------
1006	!! Default option Dummy module NO SI3 sea-ice model
1007	!!----------------------------------------------------------------------
1008	#endif
1009
1010	!!======================================================================
1011	END MODULE icevar

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

Download in other formats: