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/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icevar.F90 @ 12369

Last change on this file since 12369 was 12369, checked in by dancopsey, 4 years ago
Add print statements
File size: 49.4 KB

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

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

Download in other formats: