Context Navigation

icevar.F90 @ 12379

Last change on this file since 12379 was 12379, checked in by dancopsey, 4 years ago
Add meltpond lid thickness as a new prognostic.
File size: 50.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
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	REAL(wp), PARAMETER :: min_salin = 0.1_wp ! Minimum salinity in sea ice
321	REAL(wp) :: zn ! Fraction of ice depth
322	REAL(wp), PARAMETER :: msal = 0.573_wp
323	REAL(wp), PARAMETER :: nsal = 0.407_wp
324	REAL(wp), PARAMETER :: saltmax = 9.6_wp
325	! REAL(wp), PARAMETER, DIMENSION(4) :: weights = (/0.06363419, 0.25545967, 0.33155539, 0.34935076/)
326	REAL(wp), PARAMETER, DIMENSION(4) :: weights = (/0.0703805, 0.25526203, 0.32860314, 0.34575433/)
327
328	!!-------------------------------------------------------------------
329
330	!!gm Question: Remove the option 3 ? How many years since it last use ?
331
332	SELECT CASE ( nn_icesal )
333	!
334	! !---------------------------------------!
335	CASE( 1 ) ! constant salinity in time and space !
336	! !---------------------------------------!
337	sz_i(:,:,:,:) = rn_icesal
338	s_i (:,:,:) = rn_icesal
339	!
340	! !---------------------------------------------!
341	CASE( 2 ) ! time varying salinity with linear profile !
342	! !---------------------------------------------!
343	!
344	ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
345	!
346	DO jl = 1, jpl
347	DO jk = 1, nlay_i
348	sz_i(:,:,jk,jl) = s_i(:,:,jl)
349	END DO
350	END DO
351
352	! ! Slope of the linear profile
353	WHERE( h_i(:,:,:) > epsi20 ) ; z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:)
354	ELSEWHERE ; z_slope_s(:,:,:) = 0._wp
355	END WHERE
356	!
357	z1_dS = 1._wp / ( zsi1 - zsi0 )
358	DO jl = 1, jpl
359	DO jj = 1, jpj
360	DO ji = 1, jpi
361	zalpha(ji,jj,jl) = 1._wp
362	! ! force a constant profile when SSS too low (Baltic Sea)
363	IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) ) zalpha(ji,jj,jl) = 0._wp
364	END DO
365	END DO
366	END DO
367
368	!
369	! Computation of the profile
370	DO jl = 1, jpl
371	DO jk = 1, nlay_i
372	DO jj = 1, jpj
373	DO ji = 1, jpi
374
375	! Copied over from CICE
376	! zn = (real(jk,kind=dbl_kind)-0.5_wp) / real(nlay_i,kind=dbl_kind)
377	! sz_i(ji,jj,jk,jl)=(saltmax/2.0_wp)(1.0_wp-cos(pizn**(nsal/(msal+zn))))
378	! sz_i(ji,jj,jk,jl) = max(sz_i(ji,jj,jk,jl), min_salin)
379
380	! Weighting method to match CICE
381	zs0 = (s_i(ji,jj,jl)4.0_wp) weights(jk)
382	zs = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl) ! weighting the profile
383	sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
384
385
386	! zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
387	! zs = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl) ! weighting the profile
388	! sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
389
390	IF ( jl == 1 .AND. jj == 26 .AND. ji == 42 ) THEN
391	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)
392	ENDIF
393
394	END DO
395	END DO
396	END DO
397	END DO
398	!
399	DEALLOCATE( z_slope_s , zalpha )
400	!
401	! !-------------------------------------------!
402	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
403	! !-------------------------------------------! (mean = 2.30)
404	!
405	s_i(:,:,:) = 2.30_wp
406	!!gm Remark: if we keep the case 3, then compute an store one for all time-step
407	!! a array S_prof(1:nlay_i) containing the calculation and just do:
408	! DO jk = 1, nlay_i
409	! sz_i(:,:,jk,:) = S_prof(jk)
410	! END DO
411	!!gm end
412	!
413	DO jl = 1, jpl
414	DO jk = 1, nlay_i
415	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
416	sz_i(:,:,jk,jl) = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) )
417	END DO
418	END DO
419	!
420	END SELECT
421	!
422	END SUBROUTINE ice_var_salprof
423
424
425	SUBROUTINE ice_var_salprof1d
426	!!-------------------------------------------------------------------
427	!! * ROUTINE ice_var_salprof1d *
428	!!
429	!! ** Purpose : 1d computation of the sea ice salinity profile
430	!! Works with 1d vectors and is used by thermodynamic modules
431	!!-------------------------------------------------------------------
432	INTEGER :: ji, jk ! dummy loop indices
433	REAL(wp) :: zargtemp, zsal, z1_dS ! local scalars
434	REAL(wp) :: zs, zs0 ! - -
435	!
436	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_slope_s, zalpha !
437	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
438	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
439	!!-------------------------------------------------------------------
440	!
441	SELECT CASE ( nn_icesal )
442	!
443	! !---------------------------------------!
444	CASE( 1 ) ! constant salinity in time and space !
445	! !---------------------------------------!
446	sz_i_1d(1:npti,:) = rn_icesal
447	!
448	! !---------------------------------------------!
449	CASE( 2 ) ! time varying salinity with linear profile !
450	! !---------------------------------------------!
451	!
452	ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
453	!
454	! ! Slope of the linear profile
455	WHERE( h_i_1d(1:npti) > epsi20 ) ; z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
456	ELSEWHERE ; z_slope_s(1:npti) = 0._wp
457	END WHERE
458
459	z1_dS = 1._wp / ( zsi1 - zsi0 )
460	DO ji = 1, npti
461	zalpha(ji) = MAX( 0._wp , MIN( ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp ) )
462	! ! force a constant profile when SSS too low (Baltic Sea)
463	IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) ) zalpha(ji) = 0._wp
464	END DO
465	!
466	! Computation of the profile
467	DO jk = 1, nlay_i
468	DO ji = 1, npti
469	! ! linear profile with 0 surface value
470	zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
471	zs = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * s_i_1d(ji)
472	sz_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
473	END DO
474	END DO
475	!
476	DEALLOCATE( z_slope_s, zalpha )
477
478	! !-------------------------------------------!
479	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
480	! !-------------------------------------------! (mean = 2.30)
481	!
482	s_i_1d(1:npti) = 2.30_wp
483	!
484	!!gm cf remark in ice_var_salprof routine, CASE( 3 )
485	DO jk = 1, nlay_i
486	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
487	zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
488	DO ji = 1, npti
489	sz_i_1d(ji,jk) = zsal
490	END DO
491	END DO
492	!
493	END SELECT
494	!
495	END SUBROUTINE ice_var_salprof1d
496
497
498	SUBROUTINE ice_var_zapsmall
499	!!-------------------------------------------------------------------
500	!! * ROUTINE ice_var_zapsmall *
501	!!
502	!! ** Purpose : Remove too small sea ice areas and correct fluxes
503	!!-------------------------------------------------------------------
504	INTEGER :: ji, jj, jl, jk ! dummy loop indices
505	REAL(wp), DIMENSION(jpi,jpj) :: zswitch
506	!!-------------------------------------------------------------------
507	!
508	DO jl = 1, jpl !== loop over the categories ==!
509	!
510	WHERE( a_i(:,:,jl) > epsi10 ) ; h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
511	ELSEWHERE ; h_i(:,:,jl) = 0._wp
512	END WHERE
513	!
514	WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 ) ; zswitch(:,:) = 0._wp
515	ELSEWHERE ; zswitch(:,:) = 1._wp
516	END WHERE
517	!
518	!-----------------------------------------------------------------
519	! Zap ice energy and use ocean heat to melt ice
520	!-----------------------------------------------------------------
521	DO jk = 1, nlay_i
522	DO jj = 1 , jpj
523	DO ji = 1 , jpi
524	! update exchanges with ocean
525	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
526	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
527	t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
528	END DO
529	END DO
530	END DO
531	!
532	DO jk = 1, nlay_s
533	DO jj = 1 , jpj
534	DO ji = 1 , jpi
535	! update exchanges with ocean
536	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
537	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj)
538	t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
539
540	If (to_print_2d(ji,jj) == 10 .AND. jk == 1 .AND. jl == 1) THEN
541	write(numout,*)'ice_var_glo2eqv 1: t_s, e_s = ',t_s(ji,jj,jk,jl), ' ', e_s (ji,jj,jk,jl)
542	END IF
543	END DO
544	END DO
545	END DO
546	!
547	!-----------------------------------------------------------------
548	! zap ice and snow volume, add water and salt to ocean
549	!-----------------------------------------------------------------
550	DO jj = 1 , jpj
551	DO ji = 1 , jpi
552	! update exchanges with ocean
553	sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoi * r1_rdtice
554	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoi * r1_rdtice
555	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhos * r1_rdtice
556	!
557	a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
558	v_i (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
559	v_s (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
560	t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
561	oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
562	sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
563	!
564	h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
565	h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
566	!
567	a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
568	v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
569	lh_ip (ji,jj,jl) = lh_ip (ji,jj,jl) * zswitch(ji,jj)
570	!
571	END DO
572	END DO
573	!
574	END DO
575
576	! to be sure that at_i is the sum of a_i(jl)
577	at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
578	vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
579
580	! open water = 1 if at_i=0
581	WHERE( at_i(:,:) == 0._wp ) ato_i(:,:) = 1._wp
582	!
583	END SUBROUTINE ice_var_zapsmall
584
585
586	SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i )
587	!!-------------------------------------------------------------------
588	!! * ROUTINE ice_var_zapneg *
589	!!
590	!! ** Purpose : Remove negative sea ice fields and correct fluxes
591	!!-------------------------------------------------------------------
592	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
593	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area
594	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume
595	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_s ! snw volume
596	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psv_i ! salt content
597	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content
598	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration
599	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
600	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume
601	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: plh_ip ! melt pond lid thickness
602	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content
603	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content
604	!
605	INTEGER :: ji, jj, jl, jk ! dummy loop indices
606	REAL(wp) :: z1_dt
607	!!-------------------------------------------------------------------
608	!
609	z1_dt = 1._wp / pdt
610	!
611	DO jl = 1, jpl !== loop over the categories ==!
612	!
613	! make sure a_i=0 where v_i<=0
614	WHERE( pv_i(:,:,:) <= 0._wp ) pa_i(:,:,:) = 0._wp
615
616	!----------------------------------------
617	! zap ice energy and send it to the ocean
618	!----------------------------------------
619	DO jk = 1, nlay_i
620	DO jj = 1 , jpj
621	DO ji = 1 , jpi
622	IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
623	hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
624	pe_i(ji,jj,jk,jl) = 0._wp
625	ENDIF
626	END DO
627	END DO
628	END DO
629	!
630	DO jk = 1, nlay_s
631	DO jj = 1 , jpj
632	DO ji = 1 , jpi
633	IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
634	hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
635	pe_s(ji,jj,jk,jl) = 0._wp
636	ENDIF
637	END DO
638	END DO
639	END DO
640	!
641	!-----------------------------------------------------
642	! zap ice and snow volume, add water and salt to ocean
643	!-----------------------------------------------------
644	DO jj = 1 , jpj
645	DO ji = 1 , jpi
646	IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
647	wfx_res(ji,jj) = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
648	pv_i (ji,jj,jl) = 0._wp
649	ENDIF
650	IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
651	wfx_res(ji,jj) = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt
652	pv_s (ji,jj,jl) = 0._wp
653	ENDIF
654	IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
655	sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt
656	psv_i (ji,jj,jl) = 0._wp
657	ENDIF
658	END DO
659	END DO
660	!
661	END DO
662	!
663	WHERE( pato_i(:,:) < 0._wp ) pato_i(:,:) = 0._wp
664	WHERE( poa_i (:,:,:) < 0._wp ) poa_i (:,:,:) = 0._wp
665	WHERE( pa_i (:,:,:) < 0._wp ) pa_i (:,:,:) = 0._wp
666	WHERE( pa_ip (:,:,:) < 0._wp ) pa_ip (:,:,:) = 0._wp
667	WHERE( pv_ip (:,:,:) < 0._wp ) pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
668	! but it does not change conservation, so keep it this way is ok
669	WHERE( plh_ip (:,:,:) < 0._wp ) plh_ip (:,:,:) = 0._wp
670	!
671	END SUBROUTINE ice_var_zapneg
672
673
674	SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i )
675	!!-------------------------------------------------------------------
676	!! * ROUTINE ice_var_roundoff *
677	!!
678	!! ** Purpose : Remove negative sea ice values arising from roundoff errors
679	!!-------------------------------------------------------------------
680	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_i ! ice concentration
681	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_i ! ice volume
682	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_s ! snw volume
683	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: psv_i ! salt content
684	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: poa_i ! age content
685	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
686	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_ip ! melt pond volume
687	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: plh_ip ! melt pond lid thickness
688	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_s ! snw heat content
689	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_i ! ice heat content
690	!!-------------------------------------------------------------------
691	!
692	WHERE( pa_i (1:npti,:) < 0._wp .AND. pa_i (1:npti,:) > -epsi10 ) pa_i (1:npti,:) = 0._wp ! a_i must be >= 0
693	WHERE( pv_i (1:npti,:) < 0._wp .AND. pv_i (1:npti,:) > -epsi10 ) pv_i (1:npti,:) = 0._wp ! v_i must be >= 0
694	WHERE( pv_s (1:npti,:) < 0._wp .AND. pv_s (1:npti,:) > -epsi10 ) pv_s (1:npti,:) = 0._wp ! v_s must be >= 0
695	WHERE( psv_i(1:npti,:) < 0._wp .AND. psv_i(1:npti,:) > -epsi10 ) psv_i(1:npti,:) = 0._wp ! sv_i must be >= 0
696	WHERE( poa_i(1:npti,:) < 0._wp .AND. poa_i(1:npti,:) > -epsi10 ) poa_i(1:npti,:) = 0._wp ! oa_i must be >= 0
697	WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 ) pe_i (1:npti,:,:) = 0._wp ! e_i must be >= 0
698	WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 ) pe_s (1:npti,:,:) = 0._wp ! e_s must be >= 0
699	IF ( ln_pnd_H12 ) THEN
700	WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 ) pa_ip(1:npti,:) = 0._wp ! a_ip must be >= 0
701	WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 ) pv_ip(1:npti,:) = 0._wp ! v_ip must be >= 0
702	WHERE( plh_ip(1:npti,:) < 0._wp .AND. plh_ip(1:npti,:) > -epsi10 ) plh_ip(1:npti,:) = 0._wp ! lh_ip must be >= 0
703	ENDIF
704	!
705	END SUBROUTINE ice_var_roundoff
706
707
708	SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
709	!!-------------------------------------------------------------------
710	!! * ROUTINE ice_var_itd *
711	!!
712	!! ** Purpose : converting 1-cat ice to multiple ice categories
713	!!
714	!! ice thickness distribution follows a gaussian law
715	!! around the concentration of the most likely ice thickness
716	!! (similar as iceistate.F90)
717	!!
718	!! ** Method: Iterative procedure
719	!!
720	!! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
721	!!
722	!! 2) Check whether the distribution conserves area and volume, positivity and
723	!! category boundaries
724	!!
725	!! 3) If not (input ice is too thin), the last category is empty and
726	!! the number of categories is reduced (jpl-1)
727	!!
728	!! 4) Iterate until ok (SUM(itest(:) = 4)
729	!!
730	!! ** Arguments : zhti: 1-cat ice thickness
731	!! zhts: 1-cat snow depth
732	!! zati: 1-cat ice concentration
733	!!
734	!! ** Output : jpl-cat
735	!!
736	!! (Example of application: BDY forcings when input are cell averaged)
737	!!-------------------------------------------------------------------
738	INTEGER :: ji, jk, jl ! dummy loop indices
739	INTEGER :: idim, i_fill, jl0
740	REAL(wp) :: zarg, zV, zconv, zdh, zdv
741	REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables
742	REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables
743	INTEGER , DIMENSION(4) :: itest
744	!!-------------------------------------------------------------------
745	!
746	! ----------------------------------------
747	! distribution over the jpl ice categories
748	! ----------------------------------------
749	! a gaussian distribution for ice concentration is used
750	! then we check whether the distribution fullfills
751	! volume and area conservation, positivity and ice categories bounds
752	idim = SIZE( zhti , 1 )
753	zh_i(1:idim,1:jpl) = 0._wp
754	zh_s(1:idim,1:jpl) = 0._wp
755	za_i(1:idim,1:jpl) = 0._wp
756	!
757	DO ji = 1, idim
758	!
759	IF( zhti(ji) > 0._wp ) THEN
760	!
761	! find which category (jl0) the input ice thickness falls into
762	jl0 = jpl
763	DO jl = 1, jpl
764	IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
765	jl0 = jl
766	CYCLE
767	ENDIF
768	END DO
769	!
770	itest(:) = 0
771	i_fill = jpl + 1 !------------------------------------
772	DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories
773	! !------------------------------------
774	i_fill = i_fill - 1
775	!
776	zh_i(ji,1:jpl) = 0._wp
777	za_i(ji,1:jpl) = 0._wp
778	itest(:) = 0
779	!
780	IF ( i_fill == 1 ) THEN !-- case very thin ice: fill only category 1
781	zh_i(ji,1) = zhti(ji)
782	za_i (ji,1) = zati (ji)
783	ELSE !-- case ice is thicker: fill categories >1
784	! thickness
785	DO jl = 1, i_fill - 1
786	zh_i(ji,jl) = hi_mean(jl)
787	END DO
788	!
789	! concentration
790	za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
791	DO jl = 1, i_fill - 1
792	IF ( jl /= jl0 ) THEN
793	zarg = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
794	za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2)
795	ENDIF
796	END DO
797	!
798	! last category
799	za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
800	zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
801	zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )
802	!
803	! correction if concentration of upper cat is greater than lower cat
804	! (it should be a gaussian around jl0 but sometimes it is not)
805	IF ( jl0 /= jpl ) THEN
806	DO jl = jpl, jl0+1, -1
807	IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
808	zdv = zh_i(ji,jl) * za_i(ji,jl)
809	zh_i(ji,jl ) = 0._wp
810	za_i (ji,jl ) = 0._wp
811	za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
812	END IF
813	END DO
814	ENDIF
815	!
816	ENDIF
817	!
818	! Compatibility tests
819	zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) )
820	IF ( zconv < epsi06 ) itest(1) = 1 ! Test 1: area conservation
821	!
822	zconv = ABS( zhti(ji)zati(ji) - SUM( za_i(ji,1:jpl)zh_i(ji,1:jpl) ) )
823	IF ( zconv < epsi06 ) itest(2) = 1 ! Test 2: volume conservation
824	!
825	IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ?
826	!
827	itest(4) = 1
828	DO jl = 1, i_fill
829	IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations
830	END DO
831	! !----------------------------
832	END DO ! end iteration on categories
833	! !----------------------------
834	ENDIF
835	END DO
836
837	! Add Snow in each category where za_i is not 0
838	DO jl = 1, jpl
839	DO ji = 1, idim
840	IF( za_i(ji,jl) > 0._wp ) THEN
841	zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
842	! In case snow load is in excess that would lead to transformation from snow to ice
843	! Then, transfer the snow excess into the ice (different from icethd_dh)
844	zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 )
845	! recompute h_i, h_s avoiding out of bounds values
846	zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
847	zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos )
848	ENDIF
849	END DO
850	END DO
851	!
852	END SUBROUTINE ice_var_itd
853
854
855	SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )
856	!!-------------------------------------------------------------------
857	!! * ROUTINE ice_var_itd2 *
858	!!
859	!! ** Purpose : converting N-cat ice to jpl ice categories
860	!!
861	!! ice thickness distribution follows a gaussian law
862	!! around the concentration of the most likely ice thickness
863	!! (similar as iceistate.F90)
864	!!
865	!! ** Method: Iterative procedure
866	!!
867	!! 1) Fill ice cat that correspond to input thicknesses
868	!! Find the lowest(jlmin) and highest(jlmax) cat that are filled
869	!!
870	!! 2) Expand the filling to the cat jlmin-1 and jlmax+1
871	!! by removing 25% ice area from jlmin and jlmax (resp.)
872	!!
873	!! 3) Expand the filling to the empty cat between jlmin and jlmax
874	!! by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
875	!! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
876	!!
877	!! ** Arguments : zhti: N-cat ice thickness
878	!! zhts: N-cat snow depth
879	!! zati: N-cat ice concentration
880	!!
881	!! ** Output : jpl-cat
882	!!
883	!! (Example of application: BDY forcings when inputs have N-cat /= jpl)
884	!!-------------------------------------------------------------------
885	INTEGER :: ji, jl, jl1, jl2 ! dummy loop indices
886	INTEGER :: idim, icat
887	REAL(wp), PARAMETER :: ztrans = 0.25_wp
888	REAL(wp), DIMENSION(:,:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables
889	REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables
890	INTEGER , DIMENSION(:,:), ALLOCATABLE :: jlfil, jlfil2
891	INTEGER , DIMENSION(:) , ALLOCATABLE :: jlmax, jlmin
892	!!-------------------------------------------------------------------
893	!
894	idim = SIZE( zhti, 1 )
895	icat = SIZE( zhti, 2 )
896	!
897	ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays
898	ALLOCATE( jlmin(idim), jlmax(idim) )
899
900	! --- initialize output fields to 0 --- !
901	zh_i(1:idim,1:jpl) = 0._wp
902	zh_s(1:idim,1:jpl) = 0._wp
903	za_i(1:idim,1:jpl) = 0._wp
904	!
905	! --- fill the categories --- !
906	! find where cat-input = cat-output and fill cat-output fields
907	jlmax(:) = 0
908	jlmin(:) = 999
909	jlfil(:,:) = 0
910	DO jl1 = 1, jpl
911	DO jl2 = 1, icat
912	DO ji = 1, idim
913	IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
914	! fill the right category
915	zh_i(ji,jl1) = zhti(ji,jl2)
916	zh_s(ji,jl1) = zhts(ji,jl2)
917	za_i(ji,jl1) = zati(ji,jl2)
918	! record categories that are filled
919	jlmax(ji) = MAX( jlmax(ji), jl1 )
920	jlmin(ji) = MIN( jlmin(ji), jl1 )
921	jlfil(ji,jl1) = jl1
922	ENDIF
923	END DO
924	END DO
925	END DO
926	!
927	! --- fill the gaps between categories --- !
928	! transfer from categories filled at the previous step to the empty ones in between
929	DO ji = 1, idim
930	jl1 = jlmin(ji)
931	jl2 = jlmax(ji)
932	IF( jl1 > 1 ) THEN
933	! fill the lower cat (jl1-1)
934	za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
935	zh_i(ji,jl1-1) = hi_mean(jl1-1)
936	! remove from cat jl1
937	za_i(ji,jl1 ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
938	ENDIF
939	IF( jl2 < jpl ) THEN
940	! fill the upper cat (jl2+1)
941	za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
942	zh_i(ji,jl2+1) = hi_mean(jl2+1)
943	! remove from cat jl2
944	za_i(ji,jl2 ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
945	ENDIF
946	END DO
947	!
948	jlfil2(:,:) = jlfil(:,:)
949	! fill categories from low to high
950	DO jl = 2, jpl-1
951	DO ji = 1, idim
952	IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
953	! fill high
954	za_i(ji,jl) = ztrans * za_i(ji,jl-1)
955	zh_i(ji,jl) = hi_mean(jl)
956	jlfil(ji,jl) = jl
957	! remove low
958	za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
959	ENDIF
960	END DO
961	END DO
962	!
963	! fill categories from high to low
964	DO jl = jpl-1, 2, -1
965	DO ji = 1, idim
966	IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
967	! fill low
968	za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
969	zh_i(ji,jl) = hi_mean(jl)
970	jlfil2(ji,jl) = jl
971	! remove high
972	za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
973	ENDIF
974	END DO
975	END DO
976	!
977	DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays
978	DEALLOCATE( jlmin, jlmax )
979	!
980	END SUBROUTINE ice_var_itd2
981
982
983	SUBROUTINE ice_var_bv
984	!!-------------------------------------------------------------------
985	!! * ROUTINE ice_var_bv *
986	!!
987	!! ** Purpose : computes mean brine volume (%) in sea ice
988	!!
989	!! ** Method : e = - 0.054 * S (ppt) / T (C)
990	!!
991	!! References : Vancoppenolle et al., JGR, 2007
992	!!-------------------------------------------------------------------
993	INTEGER :: ji, jj, jk, jl ! dummy loop indices
994	!!-------------------------------------------------------------------
995	!
996	!!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done
997	!! instead of setting everything to zero as just below
998	bv_i (:,:,:) = 0._wp
999	DO jl = 1, jpl
1000	DO jk = 1, nlay_i
1001	WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )
1002	bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
1003	END WHERE
1004	END DO
1005	END DO
1006	WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
1007	ELSEWHERE ; bvm_i(:,:) = 0._wp
1008	END WHERE
1009	!
1010	END SUBROUTINE ice_var_bv
1011
1012
1013	SUBROUTINE ice_var_enthalpy
1014	!!-------------------------------------------------------------------
1015	!! * ROUTINE ice_var_enthalpy *
1016	!!
1017	!! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature
1018	!!
1019	!! ** Method : Formula (Bitz and Lipscomb, 1999)
1020	!!-------------------------------------------------------------------
1021	INTEGER :: ji, jk ! dummy loop indices
1022	REAL(wp) :: ztmelts ! local scalar
1023	!!-------------------------------------------------------------------
1024	!
1025	DO jk = 1, nlay_i ! Sea ice energy of melting
1026	DO ji = 1, npti
1027	ztmelts = - rTmlt * sz_i_1d(ji,jk)
1028	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
1029	! (sometimes zdf scheme produces abnormally high temperatures)
1030	e_i_1d(ji,jk) = rhoi * ( rcpi * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) ) &
1031	& + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) ) &
1032	& - rcp * ztmelts )
1033	END DO
1034	END DO
1035	DO jk = 1, nlay_s ! Snow energy of melting
1036	DO ji = 1, npti
1037	e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
1038	END DO
1039	END DO
1040	!
1041	END SUBROUTINE ice_var_enthalpy
1042
1043	FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
1044	!!---------------------------------------------------------------------
1045	!! * ROUTINE ice_var_sshdyn *
1046	!!
1047	!! ** Purpose : compute the equivalent ssh in lead when sea ice is embedded
1048	!!
1049	!! ** Method : ssh_lead = ssh + (Mice + Msnow) / rau0
1050	!!
1051	!! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira,
1052	!! Sea ice-ocean coupling using a rescaled vertical coordinate z*,
1053	!! Ocean Modelling, Volume 24, Issues 1-2, 2008
1054	!!----------------------------------------------------------------------
1055	!
1056	! input
1057	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh !: ssh [m]
1058	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass !: mass of snow and ice at current ice time step [Kg/m2]
1059	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b !: mass of snow and ice at previous ice time step [Kg/m2]
1060	!
1061	! output
1062	REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn ! equivalent ssh in lead [m]
1063	!
1064	! temporary
1065	REAL(wp) :: zintn, zintb ! time interpolation weights []
1066	REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload ! snow and ice load [m]
1067	!
1068	! compute ice load used to define the equivalent ssh in lead
1069	IF( ln_ice_embd ) THEN
1070	!
1071	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
1072	! = (1/nn_fsbc)^2 * {SUM[n] , n=0,nn_fsbc-1}
1073	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
1074	!
1075	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
1076	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
1077	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
1078	!
1079	zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0
1080	!
1081	ELSE
1082	zsnwiceload(:,:) = 0.0_wp
1083	ENDIF
1084	! compute equivalent ssh in lead
1085	ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:)
1086	!
1087	END FUNCTION ice_var_sshdyn
1088
1089
1090	#else
1091	!!----------------------------------------------------------------------
1092	!! Default option Dummy module NO SI3 sea-ice model
1093	!!----------------------------------------------------------------------
1094	#endif
1095
1096	!!======================================================================
1097	END MODULE icevar

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

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

Download in other formats: