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 @ 12371

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

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

Download in other formats: