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.

icedyn_adv_umx.F90 in NEMO/trunk/src/ICE – NEMO

Context Navigation

source: NEMO/trunk/src/ICE/icedyn_adv_umx.F90 @ 10911

Last change on this file since 10911 was 10911, checked in by clem, 5 years ago
Major change: the advection scheme UMx has been revisited to clean all the unphysical values which occured. Minor change: the ref config SPITZ12 has been slightly modified to test more parameterizations that are available in the code (melt ponds and landfast)
Property svn:keywords set to `Id`
File size: 83.3 KB

Line
1	MODULE icedyn_adv_umx
2	!!==============================================================================
3	!! * MODULE icedyn_adv_umx *
4	!! sea-ice : advection using the ULTIMATE-MACHO scheme
5	!!==============================================================================
6	!! History : 3.6 ! 2014-11 (C. Rousset, G. Madec) 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	!! ice_dyn_adv_umx : update the tracer fields
14	!! ultimate_x(_y) : compute a tracer value at velocity points using ULTIMATE scheme at various orders
15	!! macho : compute the fluxes
16	!! nonosc_ice : limit the fluxes using a non-oscillatory algorithm
17	!!----------------------------------------------------------------------
18	USE phycst ! physical constant
19	USE dom_oce ! ocean domain
20	USE sbc_oce , ONLY : nn_fsbc ! update frequency of surface boundary condition
21	USE ice ! sea-ice variables
22	USE icevar ! sea-ice: operations
23	!
24	USE in_out_manager ! I/O manager
25	USE iom ! I/O manager library
26	USE lib_mpp ! MPP library
27	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
28	USE lbclnk ! lateral boundary conditions (or mpp links)
29
30	IMPLICIT NONE
31	PRIVATE
32
33	PUBLIC ice_dyn_adv_umx ! called by icedyn_adv.F90
34
35	REAL(wp) :: z1_6 = 1._wp / 6._wp ! =1/6
36	REAL(wp) :: z1_120 = 1._wp / 120._wp ! =1/120
37
38	! limiter: 1=nonosc_ice, 2=superbee, 3=h3(rachid)
39	INTEGER :: kn_limiter = 1
40
41	! if there is an outward velocity in a grid cell where there is no ice initially (typically at the ice edge),
42	! interpolated T at u/v points can be non-zero while it should
43	! (because of the high order of the advection scheme). Thus set it to 0 in this case
44	LOGICAL :: ll_icedge = .TRUE.
45
46	! if T interpolated at u/v points is negative or v_i < 1.e-6
47	! then interpolate T at u/v points using the upstream scheme
48	LOGICAL :: ll_neg = .TRUE.
49
50	! alternate directions for upstream
51	LOGICAL :: ll_upsxy = .TRUE.
52
53	! alternate directions for high order
54	LOGICAL :: ll_hoxy = .TRUE.
55
56	! prelimiter: use it to avoid overshoot in H
57	LOGICAL :: ll_prelimiter_zalesak = .FALSE. ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D
58
59	! advection for S and T: dVS/dt = -div( uA * uHS / u ) => ll_advS = F
60	! or dVS/dt = -div( uV * uS / u ) => ll_advS = T
61	LOGICAL :: ll_advS = .FALSE.
62
63	!! * Substitutions
64	# include "vectopt_loop_substitute.h90"
65	!!----------------------------------------------------------------------
66	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
67	!! $Id$
68	!! Software governed by the CeCILL licence (./LICENSE)
69	!!----------------------------------------------------------------------
70	CONTAINS
71
72	SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip, &
73	& pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
74	!!----------------------------------------------------------------------
75	!! * ROUTINE ice_dyn_adv_umx *
76	!!
77	!! ** Purpose : Compute the now trend due to total advection of
78	!! tracers and add it to the general trend of tracer equations
79	!! using an "Ultimate-Macho" scheme
80	!!
81	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
82	!!----------------------------------------------------------------------
83	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
84	INTEGER , INTENT(in ) :: kt ! time step
85	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity
86	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pv_ice ! ice j-velocity
87	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_i ! ice thickness
88	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_s ! snw thickness
89	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_ip ! ice pond thickness
90	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area
91	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume
92	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_s ! snw volume
93	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psv_i ! salt content
94	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content
95	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration
96	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
97	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume
98	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content
99	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content
100	!
101	INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices
102	INTEGER :: icycle ! number of sub-timestep for the advection
103	REAL(wp) :: zamsk ! 1 if advection of concentration, 0 if advection of other tracers
104	REAL(wp) :: zdt
105	REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication
106	REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box
107	REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2
108	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zua_ho, zva_ho, zua_ups, zva_ups
109	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zuv_ho, zvv_ho, zuv_ups, zvv_ups
110	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_ai, z1_aip, z1_vi, z1_vs
111	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhvar
112	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max
113	!!----------------------------------------------------------------------
114	!
115	IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
116	!
117	! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- !
118	DO jl = 1, jpl
119	DO jj = 2, jpjm1
120	DO ji = fs_2, fs_jpim1
121	zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), &
122	& ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), &
123	& ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), &
124	& ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) )
125	zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), &
126	& ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), &
127	& ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), &
128	& ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) )
129	zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), &
130	& ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), &
131	& ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), &
132	& ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) )
133	END DO
134	END DO
135	END DO
136	CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. )
137	!
138	!
139	! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- !
140	! Note: the advection split is applied at the next time-step in order to avoid blocking global comm.
141	! this should not affect too much the stability
142	zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
143	zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
144
145	! non-blocking global communication send zcflnow and receive zcflprv
146	CALL mpp_delay_max( 'icedyn_adv_umx', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 )
147
148	IF( zcflprv(1) > .5 ) THEN ; icycle = 2
149	ELSE ; icycle = 1
150	ENDIF
151	zdt = rdt_ice / REAL(icycle)
152
153	! --- transport --- !
154	zudy(:,:) = pu_ice(:,:) * e2u(:,:)
155	zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
156
157	! --- define velocity for advection: u*grad(H) --- !
158	DO jj = 2, jpjm1
159	DO ji = fs_2, fs_jpim1
160	IF ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN ; zcu_box(ji,jj) = 0._wp
161	ELSEIF( pu_ice(ji,jj) > 0._wp ) THEN ; zcu_box(ji,jj) = pu_ice(ji-1,jj)
162	ELSE ; zcu_box(ji,jj) = pu_ice(ji ,jj)
163	ENDIF
164
165	IF ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN ; zcv_box(ji,jj) = 0._wp
166	ELSEIF( pv_ice(ji,jj) > 0._wp ) THEN ; zcv_box(ji,jj) = pv_ice(ji,jj-1)
167	ELSE ; zcv_box(ji,jj) = pv_ice(ji,jj )
168	ENDIF
169	END DO
170	END DO
171
172	!---------------!
173	!== advection ==!
174	!---------------!
175	DO jt = 1, icycle
176
177	! record at_i before advection (for open water)
178	zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
179
180	! inverse of A and Ap
181	WHERE( pa_i(:,:,:) >= epsi20 ) ; z1_ai(:,:,:) = 1._wp / pa_i(:,:,:)
182	ELSEWHERE ; z1_ai(:,:,:) = 0.
183	END WHERE
184	WHERE( pa_ip(:,:,:) >= epsi20 ) ; z1_aip(:,:,:) = 1._wp / pa_ip(:,:,:)
185	ELSEWHERE ; z1_aip(:,:,:) = 0.
186	END WHERE
187	!
188	! set u*A=u for advection of A only
189	DO jl = 1, jpl
190	zua_ho(:,:,jl) = zudy(:,:)
191	zva_ho(:,:,jl) = zvdx(:,:)
192	END DO
193
194	!== Ice area ==!
195	zamsk = 1._wp
196	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
197	& pa_i, pa_i, zua_ups, zva_ups, zua_ho , zva_ho )
198	zamsk = 0._wp
199
200	IF( .NOT. ll_advS ) THEN !-- advection form: uA * uHS / u
201	!== Ice volume ==!
202	zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
203	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
204	& zhvar, pv_i, zua_ups, zva_ups )
205	!
206	!== Snw volume ==!
207	zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
208	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
209	& zhvar, pv_s, zua_ups, zva_ups )
210	!
211	!== Salt content ==!
212	zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:)
213	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
214	& zhvar, psv_i, zua_ups, zva_ups )
215	!
216	!== Ice heat content ==!
217	DO jk = 1, nlay_i
218	zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:)
219	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, &
220	& zhvar, pe_i(:,:,jk,:), zua_ups, zva_ups )
221	END DO
222	!
223	!== Snw heat content ==!
224	DO jk = 1, nlay_s
225	zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:)
226	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, &
227	& zhvar, pe_s(:,:,jk,:), zua_ups, zva_ups )
228	END DO
229	!
230	ELSE !-- advection form: uV * uS / u
231	! inverse of Vi
232	WHERE( pv_i(:,:,:) >= epsi20 ) ; z1_vi(:,:,:) = 1._wp / pv_i(:,:,:)
233	ELSEWHERE ; z1_vi(:,:,:) = 0.
234	END WHERE
235	! inverse of Vs
236	WHERE( pv_s(:,:,:) >= epsi20 ) ; z1_vs(:,:,:) = 1._wp / pv_s(:,:,:)
237	ELSEWHERE ; z1_vs(:,:,:) = 0.
238	END WHERE
239	!
240	! It is important to first calculate the ice fields and then the snow fields (because we use the same arrays)
241	!
242	!== Ice volume ==!
243	zuv_ups = zua_ups
244	zvv_ups = zva_ups
245	zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
246	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
247	& zhvar, pv_i, zuv_ups, zvv_ups, zuv_ho , zvv_ho )
248	!
249	!== Salt content ==!
250	zhvar(:,:,:) = psv_i(:,:,:) * z1_vi(:,:,:)
251	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zuv_ho , zvv_ho , zcu_box, zcv_box, &
252	& zhvar, psv_i, zuv_ups, zvv_ups )
253	!
254	!== Ice heat content ==!
255	DO jk = 1, nlay_i
256	zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_vi(:,:,:)
257	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, &
258	& zhvar, pe_i(:,:,jk,:), zuv_ups, zvv_ups )
259	END DO
260	!
261	!== Snow volume ==!
262	zuv_ups = zua_ups
263	zvv_ups = zva_ups
264	zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
265	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
266	& zhvar, pv_s, zuv_ups, zvv_ups, zuv_ho , zvv_ho )
267	!
268	!== Snw heat content ==!
269	DO jk = 1, nlay_s
270	zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_vs(:,:,:)
271	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, &
272	& zhvar, pe_s(:,:,jk,:), zuv_ups, zvv_ups )
273	END DO
274	!
275	ENDIF
276	!
277	!== Ice age ==!
278	IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN
279	! clem: in theory we should use the formulation below to advect the ice age, but the code is unable to deal with
280	! fields that do not depend on volume (here oa_i depends on concentration). It creates abnormal ages that
281	! spread into the domain. Therefore we cheat and consider that ice age should be advected as ice concentration
282	!!zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:)
283	!!CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, poa_i )
284	! set u*A=u for advection of ice age
285	DO jl = 1, jpl
286	zua_ho(:,:,jl) = zudy(:,:)
287	zva_ho(:,:,jl) = zvdx(:,:)
288	END DO
289	zamsk = 1._wp
290	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho, zva_ho, zcu_box, zcv_box, &
291	& poa_i, poa_i )
292	zamsk = 0._wp
293	ENDIF
294	!
295	!== melt ponds ==!
296	IF ( ln_pnd_H12 ) THEN
297	! set u*A=u for advection of Ap only
298	DO jl = 1, jpl
299	zua_ho(:,:,jl) = zudy(:,:)
300	zva_ho(:,:,jl) = zvdx(:,:)
301	END DO
302	! fraction
303	zamsk = 1._wp
304	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
305	& pa_ip, pa_ip, zua_ups, zva_ups, zua_ho , zva_ho )
306	zamsk = 0._wp
307	! volume
308	zhvar(:,:,:) = pv_ip(:,:,:) * z1_aip(:,:,:)
309	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
310	& zhvar, pv_ip, zua_ups, zva_ups )
311	ENDIF
312	!
313	!== Open water area ==!
314	zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
315	DO jj = 2, jpjm1
316	DO ji = fs_2, fs_jpim1
317	pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &
318	& - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
319	END DO
320	END DO
321	CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T', 1. )
322	!
323	!
324	! --- Ensure non-negative fields and in-bound thicknesses --- !
325	!
326	! Remove negative values (conservation is ensured)
327	! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
328	CALL ice_var_zapneg( pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
329	!
330	! Make sure ice thickness is not too big
331	! (because ice thickness can be too large where ice concentration is very small)
332	CALL Hbig( zhi_max, zhs_max, zhip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
333
334	END DO
335	!
336	END SUBROUTINE ice_dyn_adv_umx
337
338
339	SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, &
340	& pt, ptc, pua_ups, pva_ups, pua_ho, pva_ho )
341	!!----------------------------------------------------------------------
342	!! * ROUTINE adv_umx *
343	!!
344	!! ** Purpose : Compute the now trend due to total advection of
345	!! tracers and add it to the general trend of tracer equations
346	!!
347	!! ** Method : - calculate upstream fluxes and upstream solution for tracers V/A(=H) etc
348	!! - calculate tracer H at u and v points (Ultimate)
349	!! - calculate the high order fluxes using alterning directions (Macho)
350	!! - apply a limiter on the fluxes (nonosc_ice)
351	!! - convert this tracer flux to a "volume" flux (uH -> uV)
352	!! - apply a limiter a second time on the volumes fluxes (nonosc_ice)
353	!! - calculate the high order solution for V
354	!!
355	!! ** Action : solve 3 equations => a) dA/dt = -div(uA)
356	!! b) dV/dt = -div(uV) using dH/dt = -u.grad(H)
357	!! c) dVS/dt = -div(uVS) using either dHS/dt = -u.grad(HS) or dS/dt = -u.grad(S)
358	!!
359	!! in eq. b), - fluxes uH are evaluated (with UMx) and limited with nonosc_ice. This step is necessary to get a good H.
360	!! - then we convert this flux to a "volume" flux this way => uH * uA / u
361	!! where uA is the flux from eq. a)
362	!! this "volume" flux is also limited with nonosc_ice (otherwise overshoots can occur)
363	!! - at last we estimate dV/dt = -div(uH * uA / u)
364	!!
365	!! in eq. c), one can solve the equation for S (ln_advS=T), then dVS/dt = -div(uV * uS / u)
366	!! or for HS (ln_advS=F), then dVS/dt = -div(uA * uHS / u)
367	!!
368	!! ** Note : - this method can lead to tiny negative V (-1.e-20) => set it to 0 while conserving mass etc.
369	!! - At the ice edge, Ultimate scheme can lead to:
370	!! 1) negative interpolated tracers at u-v points
371	!! 2) non-zero interpolated tracers at u-v points eventhough there is no ice and velocity is outward
372	!! Solution for 1): apply an upstream scheme when it occurs. A better solution would be to degrade the order of
373	!! the scheme automatically by applying a mask of the ice cover inside Ultimate (not done).
374	!! Solution for 2): we set it to 0 in this case
375	!! - Eventhough 1D tests give very good results (typically the one from Schar & Smolarkiewiecz), the 2D is less good.
376	!! Large values of H can appear for very small ice concentration, and when it does it messes the things up since we
377	!! work on H (and not V). It is partly related to the multi-category approach
378	!! Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice
379	!! concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter
380	!! since sv_i and e_i are still good.
381	!!----------------------------------------------------------------------
382	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
383	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
384	INTEGER , INTENT(in ) :: jt ! number of sub-iteration
385	INTEGER , INTENT(in ) :: kt ! number of iteration
386	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
387	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu , pv ! 2 ice velocity components => u*e2
388	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: puc , pvc ! 2 ice velocity components => ue2 or ua*e2u
389	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pubox, pvbox ! upstream velocity
390	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pt ! tracer field
391	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: ptc ! tracer content field
392	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout), OPTIONAL :: pua_ups, pva_ups ! upstream u*a fluxes
393	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out), OPTIONAL :: pua_ho, pva_ho ! high order u*a fluxes
394	!
395	INTEGER :: ji, jj, jl ! dummy loop indices
396	REAL(wp) :: ztra ! local scalar
397	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zfu_ho , zfv_ho , zpt
398	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zfu_ups, zfv_ups, zt_ups
399	!!----------------------------------------------------------------------
400	!
401	! Upstream (_ups) fluxes
402	! -----------------------
403	CALL upstream( pamsk, jt, kt, pdt, pt, pu, pv, zt_ups, zfu_ups, zfv_ups )
404
405	! High order (_ho) fluxes
406	! -----------------------
407	SELECT CASE( kn_umx )
408	!
409	CASE ( 20 ) !== centered second order ==!
410	!
411	CALL cen2( pamsk, jt, kt, pdt, pt, pu, pv, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho )
412	!
413	CASE ( 1:5 ) !== 1st to 5th order ULTIMATE-MACHO scheme ==!
414	!
415	CALL macho( pamsk, kn_umx, jt, kt, pdt, pt, pu, pv, pubox, pvbox, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho )
416	!
417	END SELECT
418	!
419	! --ho --ho
420	! new fluxes = uH u*a / u
421	! ----------------------------
422	IF( pamsk == 0._wp ) THEN
423	DO jl = 1, jpl
424	DO jj = 1, jpjm1
425	DO ji = 1, fs_jpim1
426	IF( ABS( pu(ji,jj) ) > epsi10 ) THEN
427	zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc (ji,jj,jl) / pu(ji,jj)
428	zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj)
429	ELSE
430	zfu_ho (ji,jj,jl) = 0._wp
431	zfu_ups(ji,jj,jl) = 0._wp
432	ENDIF
433	!
434	IF( ABS( pv(ji,jj) ) > epsi10 ) THEN
435	zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc (ji,jj,jl) / pv(ji,jj)
436	zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj)
437	ELSE
438	zfv_ho (ji,jj,jl) = 0._wp
439	zfv_ups(ji,jj,jl) = 0._wp
440	ENDIF
441	END DO
442	END DO
443	END DO
444
445	! the new "volume" fluxes must also be "flux corrected"
446	! thus we calculate the upstream solution and apply a limiter again
447	DO jl = 1, jpl
448	DO jj = 2, jpjm1
449	DO ji = fs_2, fs_jpim1
450	ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) )
451	!
452	zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)
453	END DO
454	END DO
455	END DO
456	CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. )
457	!
458	IF ( kn_limiter == 1 ) THEN
459	CALL nonosc_ice( 1._wp, pdt, pu, pv, ptc, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho )
460	ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN
461	CALL limiter_x( pdt, pu, ptc, zfu_ups, zfu_ho )
462	CALL limiter_y( pdt, pv, ptc, zfv_ups, zfv_ho )
463	ENDIF
464	!
465	ENDIF
466	! --ho --ups
467	! in case of advection of A: output ua and ua
468	! -----------------------------------------------
469	IF( PRESENT( pua_ho ) ) THEN
470	DO jl = 1, jpl
471	DO jj = 1, jpjm1
472	DO ji = 1, fs_jpim1
473	pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl)
474	pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl)
475	END DO
476	END DO
477	END DO
478	ENDIF
479	!
480	! final trend with corrected fluxes
481	! ---------------------------------
482	DO jl = 1, jpl
483	DO jj = 2, jpjm1
484	DO ji = fs_2, fs_jpim1
485	ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) )
486	!
487	ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)
488	END DO
489	END DO
490	END DO
491	CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T', 1. )
492	!
493	END SUBROUTINE adv_umx
494
495
496	SUBROUTINE upstream( pamsk, jt, kt, pdt, pt, pu, pv, pt_ups, pfu_ups, pfv_ups )
497	!!---------------------------------------------------------------------
498	!! * ROUTINE upstream *
499	!!
500	!! ** Purpose : compute the upstream fluxes and upstream guess of tracer
501	!!----------------------------------------------------------------------
502	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
503	INTEGER , INTENT(in ) :: jt ! number of sub-iteration
504	INTEGER , INTENT(in ) :: kt ! number of iteration
505	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
506	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields
507	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu, pv ! 2 ice velocity components
508	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_ups ! upstream guess of tracer
509	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ups, pfv_ups ! upstream fluxes
510	!
511	INTEGER :: ji, jj, jl ! dummy loop indices
512	REAL(wp) :: ztra ! local scalar
513	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zpt
514	!!----------------------------------------------------------------------
515
516	IF( .NOT. ll_upsxy ) THEN ! no alternate directions !
517	!
518	DO jl = 1, jpl
519	DO jj = 1, jpjm1
520	DO ji = 1, fs_jpim1
521	pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl)
522	pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl)
523	END DO
524	END DO
525	END DO
526	!
527	ELSE ! alternate directions !
528	!
529	IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==!
530	!
531	DO jl = 1, jpl !-- flux in x-direction
532	DO jj = 1, jpjm1
533	DO ji = 1, fs_jpim1
534	pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl)
535	END DO
536	END DO
537	END DO
538	!
539	DO jl = 1, jpl !-- first guess of tracer from u-flux
540	DO jj = 2, jpjm1
541	DO ji = fs_2, fs_jpim1
542	ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) ) &
543	& + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk)
544	!
545	zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
546	END DO
547	END DO
548	END DO
549	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
550	!
551	DO jl = 1, jpl !-- flux in y-direction
552	DO jj = 1, jpjm1
553	DO ji = 1, fs_jpim1
554	pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl)
555	END DO
556	END DO
557	END DO
558	!
559	ELSE !== even ice time step: adv_y then adv_x ==!
560	!
561	DO jl = 1, jpl !-- flux in y-direction
562	DO jj = 1, jpjm1
563	DO ji = 1, fs_jpim1
564	pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl)
565	END DO
566	END DO
567	END DO
568	!
569	DO jl = 1, jpl !-- first guess of tracer from v-flux
570	DO jj = 2, jpjm1
571	DO ji = fs_2, fs_jpim1
572	ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) ) &
573	& + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk)
574	!
575	zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
576	END DO
577	END DO
578	END DO
579	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
580	!
581	DO jl = 1, jpl !-- flux in x-direction
582	DO jj = 1, jpjm1
583	DO ji = 1, fs_jpim1
584	pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl)
585	END DO
586	END DO
587	END DO
588	!
589	ENDIF
590
591	ENDIF
592	!
593	DO jl = 1, jpl !-- after tracer with upstream scheme
594	DO jj = 2, jpjm1
595	DO ji = fs_2, fs_jpim1
596	ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj ,jl) &
597	& + pfv_ups(ji,jj,jl) - pfv_ups(ji ,jj-1,jl) ) &
598	& + ( pu (ji,jj ) - pu (ji-1,jj ) &
599	& + pv (ji,jj ) - pv (ji ,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk)
600	!
601	pt_ups(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
602	END DO
603	END DO
604	END DO
605	CALL lbc_lnk( 'icedyn_adv_umx', pt_ups, 'T', 1. )
606
607	END SUBROUTINE upstream
608
609
610	SUBROUTINE cen2( pamsk, jt, kt, pdt, pt, pu, pv, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
611	!!---------------------------------------------------------------------
612	!! * ROUTINE cen2 *
613	!!
614	!! ** Purpose : compute the high order fluxes using a centered
615	!! second order scheme
616	!!----------------------------------------------------------------------
617	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
618	INTEGER , INTENT(in ) :: jt ! number of sub-iteration
619	INTEGER , INTENT(in ) :: kt ! number of iteration
620	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
621	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields
622	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu, pv ! 2 ice velocity components
623	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt_ups ! upstream guess of tracer
624	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pfu_ups, pfv_ups ! upstream fluxes
625	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes
626	!
627	INTEGER :: ji, jj, jl ! dummy loop indices
628	REAL(wp) :: ztra ! local scalar
629	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zpt
630	!!----------------------------------------------------------------------
631	!
632	IF( .NOT.ll_hoxy ) THEN ! no alternate directions !
633	!
634	DO jl = 1, jpl
635	DO jj = 1, jpjm1
636	DO ji = 1, fs_jpim1
637	pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj ,jl) )
638	pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji ,jj+1,jl) )
639	END DO
640	END DO
641	END DO
642	!
643	IF ( kn_limiter == 1 ) THEN
644	CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
645	ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN
646	CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
647	CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
648	ENDIF
649	!
650	ELSE ! alternate directions !
651	!
652	IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==!
653	!
654	DO jl = 1, jpl !-- flux in x-direction
655	DO jj = 1, jpjm1
656	DO ji = 1, fs_jpim1
657	pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
658	END DO
659	END DO
660	END DO
661	IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
662
663	DO jl = 1, jpl !-- first guess of tracer from u-flux
664	DO jj = 2, jpjm1
665	DO ji = fs_2, fs_jpim1
666	ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) &
667	& + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk)
668	!
669	zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
670	END DO
671	END DO
672	END DO
673	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
674
675	DO jl = 1, jpl !-- flux in y-direction
676	DO jj = 1, jpjm1
677	DO ji = 1, fs_jpim1
678	pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) )
679	END DO
680	END DO
681	END DO
682	IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
683
684	ELSE !== even ice time step: adv_y then adv_x ==!
685	!
686	DO jl = 1, jpl !-- flux in y-direction
687	DO jj = 1, jpjm1
688	DO ji = 1, fs_jpim1
689	pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
690	END DO
691	END DO
692	END DO
693	IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
694	!
695	DO jl = 1, jpl !-- first guess of tracer from v-flux
696	DO jj = 2, jpjm1
697	DO ji = fs_2, fs_jpim1
698	ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) &
699	& + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk)
700	!
701	zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
702	END DO
703	END DO
704	END DO
705	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
706	!
707	DO jl = 1, jpl !-- flux in x-direction
708	DO jj = 1, jpjm1
709	DO ji = 1, fs_jpim1
710	pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) )
711	END DO
712	END DO
713	END DO
714	IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
715
716	ENDIF
717	IF( kn_limiter == 1 ) CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
718
719	ENDIF
720
721	END SUBROUTINE cen2
722
723
724	SUBROUTINE macho( pamsk, kn_umx, jt, kt, pdt, pt, pu, pv, pubox, pvbox, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
725	!!---------------------------------------------------------------------
726	!! * ROUTINE macho *
727	!!
728	!! ** Purpose : compute the high order fluxes using Ultimate-Macho scheme
729	!!
730	!! ** Method : ...
731	!!
732	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
733	!!----------------------------------------------------------------------
734	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
735	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
736	INTEGER , INTENT(in ) :: jt ! number of sub-iteration
737	INTEGER , INTENT(in ) :: kt ! number of iteration
738	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
739	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields
740	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu, pv ! 2 ice velocity components
741	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pubox, pvbox ! upstream velocity
742	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt_ups ! upstream guess of tracer
743	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pfu_ups, pfv_ups ! upstream fluxes
744	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes
745	!
746	INTEGER :: ji, jj, jl ! dummy loop indices
747	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zt_u, zt_v, zpt
748	!!----------------------------------------------------------------------
749	!
750	IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==!
751	!
752	! !-- ultimate interpolation of pt at u-point --!
753	CALL ultimate_x( pamsk, kn_umx, pdt, pt, pu, zt_u, pfu_ho )
754	! !-- limiter in x --!
755	IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
756	! !-- advective form update in zpt --!
757	DO jl = 1, jpl
758	DO jj = 2, jpjm1
759	DO ji = fs_2, fs_jpim1
760	zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pubox(ji,jj ) * ( zt_u(ji,jj,jl) - zt_u(ji-1,jj,jl) ) * r1_e1t (ji,jj) &
761	& + pt (ji,jj,jl) * ( pu (ji,jj ) - pu (ji-1,jj ) ) * r1_e1e2t(ji,jj) &
762	& * pamsk &
763	& ) * pdt ) * tmask(ji,jj,1)
764	!!clem test
765	!!zpt(ji,jj,jl) = MAX( 0._wp, zpt(ji,jj,jl) )
766	!!clem test
767	END DO
768	END DO
769	END DO
770	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
771	!
772	! !-- ultimate interpolation of pt at v-point --!
773	IF( ll_hoxy ) THEN
774	CALL ultimate_y( pamsk, kn_umx, pdt, zpt, pv, zt_v, pfv_ho )
775	ELSE
776	CALL ultimate_y( pamsk, kn_umx, pdt, pt , pv, zt_v, pfv_ho )
777	ENDIF
778	! !-- limiter in y --!
779	IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
780	!
781	!
782	ELSE !== even ice time step: adv_y then adv_x ==!
783	!
784	! !-- ultimate interpolation of pt at v-point --!
785	CALL ultimate_y( pamsk, kn_umx, pdt, pt, pv, zt_v, pfv_ho )
786	! !-- limiter in y --!
787	IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
788	! !-- advective form update in zpt --!
789	DO jl = 1, jpl
790	DO jj = 2, jpjm1
791	DO ji = fs_2, fs_jpim1
792	zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pvbox(ji,jj ) * ( zt_v(ji,jj,jl) - zt_v(ji,jj-1,jl) ) * r1_e2t (ji,jj) &
793	& + pt (ji,jj,jl) * ( pv (ji,jj ) - pv (ji,jj-1 ) ) * r1_e1e2t(ji,jj) &
794	& * pamsk &
795	& ) * pdt ) * tmask(ji,jj,1)
796	!!clem test
797	!!zpt(ji,jj,jl) = MAX( 0._wp, zpt(ji,jj,jl) )
798	!!clem test
799	END DO
800	END DO
801	END DO
802	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
803	!
804	! !-- ultimate interpolation of pt at u-point --!
805	IF( ll_hoxy ) THEN
806	CALL ultimate_x( pamsk, kn_umx, pdt, zpt, pu, zt_u, pfu_ho )
807	ELSE
808	CALL ultimate_x( pamsk, kn_umx, pdt, pt , pu, zt_u, pfu_ho )
809	ENDIF
810	! !-- limiter in x --!
811	IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
812	!
813	ENDIF
814
815	IF( kn_limiter == 1 ) CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
816	!
817	END SUBROUTINE macho
818
819
820	SUBROUTINE ultimate_x( pamsk, kn_umx, pdt, pt, pu, pt_u, pfu_ho )
821	!!---------------------------------------------------------------------
822	!! * ROUTINE ultimate_x *
823	!!
824	!! ** Purpose : compute tracer at u-points
825	!!
826	!! ** Method : ...
827	!!
828	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
829	!!----------------------------------------------------------------------
830	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
831	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
832	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
833	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu ! ice i-velocity component
834	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields
835	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_u ! tracer at u-point
836	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho ! high order flux
837	!
838	INTEGER :: ji, jj, jl ! dummy loop indices
839	REAL(wp) :: zcu, zdx2, zdx4, zvi_cen2 ! - -
840	REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztu1, ztu2, ztu3, ztu4
841	!!----------------------------------------------------------------------
842	!
843	! !-- Laplacian in i-direction --!
844	DO jl = 1, jpl
845	DO jj = 2, jpjm1 ! First derivative (gradient)
846	DO ji = 1, fs_jpim1
847	ztu1(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
848	END DO
849	! ! Second derivative (Laplacian)
850	DO ji = fs_2, fs_jpim1
851	ztu2(ji,jj,jl) = ( ztu1(ji,jj,jl) - ztu1(ji-1,jj,jl) ) * r1_e1t(ji,jj)
852	END DO
853	END DO
854	END DO
855	CALL lbc_lnk( 'icedyn_adv_umx', ztu2, 'T', 1. )
856	!
857	! !-- BiLaplacian in i-direction --!
858	DO jl = 1, jpl
859	DO jj = 2, jpjm1 ! Third derivative
860	DO ji = 1, fs_jpim1
861	ztu3(ji,jj,jl) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
862	END DO
863	! ! Fourth derivative
864	DO ji = fs_2, fs_jpim1
865	ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj)
866	END DO
867	END DO
868	END DO
869	CALL lbc_lnk( 'icedyn_adv_umx', ztu4, 'T', 1. )
870	!
871	!
872	SELECT CASE (kn_umx )
873	!
874	CASE( 1 ) !== 1st order central TIM ==! (Eq. 21)
875	!
876	DO jl = 1, jpl
877	DO jj = 1, jpjm1
878	DO ji = 1, fs_jpim1 ! vector opt.
879	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) &
880	& - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
881	END DO
882	END DO
883	END DO
884	!
885	CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23)
886	!
887	DO jl = 1, jpl
888	DO jj = 1, jpjm1
889	DO ji = 1, fs_jpim1 ! vector opt.
890	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
891	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) &
892	& - zcu * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
893	END DO
894	END DO
895	END DO
896	!
897	CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24)
898	!
899	DO jl = 1, jpl
900	DO jj = 1, jpjm1
901	DO ji = 1, fs_jpim1 ! vector opt.
902	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
903	zdx2 = e1u(ji,jj) * e1u(ji,jj)
904	!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
905	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
906	& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
907	& + z1_6 * zdx2 * ( zcuzcu - 1._wp ) ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
908	& - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
909	END DO
910	END DO
911	END DO
912	!
913	CASE( 4 ) !== 4th order central TIM ==! (Eq. 27)
914	!
915	DO jl = 1, jpl
916	DO jj = 1, jpjm1
917	DO ji = 1, fs_jpim1 ! vector opt.
918	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
919	zdx2 = e1u(ji,jj) * e1u(ji,jj)
920	!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
921	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
922	& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
923	& + z1_6 * zdx2 * ( zcuzcu - 1._wp ) ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
924	& - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
925	END DO
926	END DO
927	END DO
928	!
929	CASE( 5 ) !== 5th order central TIM ==! (Eq. 29)
930	!
931	DO jl = 1, jpl
932	DO jj = 1, jpjm1
933	DO ji = 1, fs_jpim1 ! vector opt.
934	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
935	zdx2 = e1u(ji,jj) * e1u(ji,jj)
936	!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
937	zdx4 = zdx2 * zdx2
938	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
939	& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
940	& + z1_6 * zdx2 * ( zcuzcu - 1._wp ) ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
941	& - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) &
942	& + z1_120 * zdx4 * ( zcuzcu - 1._wp ) ( zcuzcu - 4._wp ) ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl) &
943	& - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) )
944	END DO
945	END DO
946	END DO
947	!
948	END SELECT
949	!
950	! if there is an outward velocity in a grid cell where there is no ice initially (typically at the ice edge),
951	! interpolated T at u/v points can be non-zero while it should
952	! (because of the high order of the advection scheme). Thus set it to 0 in this case
953	IF( ll_icedge ) THEN
954	DO jl = 1, jpl
955	DO jj = 1, jpjm1
956	DO ji = 1, fs_jpim1
957	IF( pt(ji,jj,jl) <= 0._wp .AND. pu(ji,jj) >= 0._wp ) THEN
958	pt_u(ji,jj,jl) = 0._wp
959	ENDIF
960	END DO
961	END DO
962	END DO
963	ENDIF
964	!
965	! if pt at u-point is negative then use the upstream value
966	! this should not be necessary if a proper sea-ice mask is set in Ultimate
967	! to degrade the order of the scheme when necessary (for ex. at the ice edge)
968	IF( ll_neg ) THEN
969	DO jl = 1, jpl
970	DO jj = 1, jpjm1
971	DO ji = 1, fs_jpim1
972	zvi_cen2 = 0.5_wp * ( v_i(ji+1,jj,jl) + v_i(ji,jj,jl) )
973	IF( pt_u(ji,jj,jl) < 0._wp .OR. ( zvi_cen2 < epsi06 .AND. pamsk == 0._wp ) ) THEN
974	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) &
975	& - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
976	ENDIF
977	END DO
978	END DO
979	END DO
980	ENDIF
981	! !-- High order flux in i-direction --!
982	DO jl = 1, jpl
983	DO jj = 1, jpjm1
984	DO ji = 1, fs_jpim1 ! vector opt.
985	pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl)
986	END DO
987	END DO
988	END DO
989	!
990	END SUBROUTINE ultimate_x
991
992
993	SUBROUTINE ultimate_y( pamsk, kn_umx, pdt, pt, pv, pt_v, pfv_ho )
994	!!---------------------------------------------------------------------
995	!! * ROUTINE ultimate_y *
996	!!
997	!! ** Purpose : compute tracer at v-points
998	!!
999	!! ** Method : ...
1000	!!
1001	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
1002	!!----------------------------------------------------------------------
1003	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
1004	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
1005	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
1006	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pv ! ice j-velocity component
1007	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields
1008	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_v ! tracer at v-point
1009	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfv_ho ! high order flux
1010	!
1011	INTEGER :: ji, jj, jl ! dummy loop indices
1012	REAL(wp) :: zcv, zdy2, zdy4, zvi_cen2 ! - -
1013	REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztv1, ztv2, ztv3, ztv4
1014	!!----------------------------------------------------------------------
1015	!
1016	! !-- Laplacian in j-direction --!
1017	DO jl = 1, jpl
1018	DO jj = 1, jpjm1 ! First derivative (gradient)
1019	DO ji = fs_2, fs_jpim1
1020	ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
1021	END DO
1022	END DO
1023	DO jj = 2, jpjm1 ! Second derivative (Laplacian)
1024	DO ji = fs_2, fs_jpim1
1025	ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj)
1026	END DO
1027	END DO
1028	END DO
1029	CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1. )
1030	!
1031	! !-- BiLaplacian in j-direction --!
1032	DO jl = 1, jpl
1033	DO jj = 1, jpjm1 ! First derivative
1034	DO ji = fs_2, fs_jpim1
1035	ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
1036	END DO
1037	END DO
1038	DO jj = 2, jpjm1 ! Second derivative
1039	DO ji = fs_2, fs_jpim1
1040	ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj)
1041	END DO
1042	END DO
1043	END DO
1044	CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1. )
1045	!
1046	!
1047	SELECT CASE (kn_umx )
1048	!
1049	CASE( 1 ) !== 1st order central TIM ==! (Eq. 21)
1050	DO jl = 1, jpl
1051	DO jj = 1, jpjm1
1052	DO ji = 1, fs_jpim1
1053	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) &
1054	& - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
1055	END DO
1056	END DO
1057	END DO
1058	!
1059	CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23)
1060	DO jl = 1, jpl
1061	DO jj = 1, jpjm1
1062	DO ji = 1, fs_jpim1
1063	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1064	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) &
1065	& - zcv * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
1066	END DO
1067	END DO
1068	END DO
1069	!
1070	CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24)
1071	DO jl = 1, jpl
1072	DO jj = 1, jpjm1
1073	DO ji = 1, fs_jpim1
1074	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1075	zdy2 = e2v(ji,jj) * e2v(ji,jj)
1076	!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
1077	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
1078	& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
1079	& + z1_6 * zdy2 * ( zcvzcv - 1._wp ) ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
1080	& - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
1081	END DO
1082	END DO
1083	END DO
1084	!
1085	CASE( 4 ) !== 4th order central TIM ==! (Eq. 27)
1086	DO jl = 1, jpl
1087	DO jj = 1, jpjm1
1088	DO ji = 1, fs_jpim1
1089	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1090	zdy2 = e2v(ji,jj) * e2v(ji,jj)
1091	!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
1092	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
1093	& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
1094	& + z1_6 * zdy2 * ( zcvzcv - 1._wp ) ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
1095	& - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
1096	END DO
1097	END DO
1098	END DO
1099	!
1100	CASE( 5 ) !== 5th order central TIM ==! (Eq. 29)
1101	DO jl = 1, jpl
1102	DO jj = 1, jpjm1
1103	DO ji = 1, fs_jpim1
1104	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1105	zdy2 = e2v(ji,jj) * e2v(ji,jj)
1106	!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
1107	zdy4 = zdy2 * zdy2
1108	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
1109	& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
1110	& + z1_6 * zdy2 * ( zcvzcv - 1._wp ) ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
1111	& - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) &
1112	& + z1_120 * zdy4 * ( zcvzcv - 1._wp ) ( zcvzcv - 4._wp ) ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl) &
1113	& - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) )
1114	END DO
1115	END DO
1116	END DO
1117	!
1118	END SELECT
1119	!
1120	! if there is an outward velocity in a grid cell where there is no ice initially (typically at the ice edge),
1121	! interpolated T at u/v points can be non-zero while it should
1122	! (because of the high order of the advection scheme). Thus set it to 0 in this case
1123	IF( ll_icedge ) THEN
1124	DO jl = 1, jpl
1125	DO jj = 1, jpjm1
1126	DO ji = 1, fs_jpim1
1127	IF( pt(ji,jj,jl) <= 0._wp .AND. pv(ji,jj) >= 0._wp ) THEN
1128	pt_v(ji,jj,jl) = 0._wp
1129	ENDIF
1130	END DO
1131	END DO
1132	END DO
1133	ENDIF
1134	!
1135	! if pt at v-point is negative then use the upstream value
1136	! this should not be necessary if a proper sea-ice mask is set in Ultimate
1137	! to degrade the order of the scheme when necessary (for ex. at the ice edge)
1138	IF( ll_neg ) THEN
1139	DO jl = 1, jpl
1140	DO jj = 1, jpjm1
1141	DO ji = 1, fs_jpim1
1142	zvi_cen2 = 0.5_wp * ( v_i(ji,jj+1,jl) + v_i(ji,jj,jl) )
1143	IF( pt_v(ji,jj,jl) < 0._wp .OR. ( zvi_cen2 < epsi06 .AND. pamsk == 0._wp ) ) THEN
1144	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) &
1145	& - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
1146	ENDIF
1147	END DO
1148	END DO
1149	END DO
1150	ENDIF
1151	! !-- High order flux in j-direction --!
1152	DO jl = 1, jpl
1153	DO jj = 1, jpjm1
1154	DO ji = 1, fs_jpim1 ! vector opt.
1155	pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl)
1156	END DO
1157	END DO
1158	END DO
1159	!
1160	END SUBROUTINE ultimate_y
1161
1162
1163	SUBROUTINE nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
1164	!!---------------------------------------------------------------------
1165	!! * ROUTINE nonosc_ice *
1166	!!
1167	!! ** Purpose : compute monotonic tracer fluxes from the upstream
1168	!! scheme and the before field by a non-oscillatory algorithm
1169	!!
1170	!! ** Method : ...
1171	!!----------------------------------------------------------------------
1172	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
1173	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
1174	REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pu ! ice i-velocity => u*e2
1175	REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pv ! ice j-velocity => v*e1
1176	REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pt, pt_ups ! before field & upstream guess of after field
1177	REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pfv_ups, pfu_ups ! upstream flux
1178	REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfv_ho, pfu_ho ! monotonic flux
1179	!
1180	INTEGER :: ji, jj, jl ! dummy loop indices
1181	REAL(wp) :: zpos, zneg, zbig, zup, zdo, z1_dt ! local scalars
1182	REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zcoef, zzt ! - -
1183	REAL(wp), DIMENSION(jpi,jpj ) :: zbup, zbdo
1184	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zbetup, zbetdo, zti_ups, ztj_ups
1185	!!----------------------------------------------------------------------
1186	zbig = 1.e+40_wp
1187
1188	! antidiffusive flux : high order minus low order
1189	! --------------------------------------------------
1190	DO jl = 1, jpl
1191	DO jj = 1, jpjm1
1192	DO ji = 1, fs_jpim1 ! vector opt.
1193	pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)
1194	pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)
1195	END DO
1196	END DO
1197	END DO
1198
1199	! extreme case where pfu_ho has to be zero
1200	! ----------------------------------------
1201	! pfu_ho
1202	! * --->
1203	! \| \| * \| \|
1204	! \| \| \| * \|
1205	! \| \| \| \| *
1206	! t_ups : i-1 i i+1 i+2
1207	IF( ll_prelimiter_zalesak ) THEN
1208
1209	DO jl = 1, jpl
1210	DO jj = 2, jpjm1
1211	DO ji = fs_2, fs_jpim1
1212	zti_ups(ji,jj,jl)= pt_ups(ji+1,jj ,jl)
1213	ztj_ups(ji,jj,jl)= pt_ups(ji ,jj+1,jl)
1214	END DO
1215	END DO
1216	END DO
1217	CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_ups, 'T', 1., ztj_ups, 'T', 1. )
1218
1219	DO jl = 1, jpl
1220	DO jj = 2, jpjm1
1221	DO ji = fs_2, fs_jpim1
1222	IF ( pfu_ho(ji,jj,jl) * ( pt_ups(ji+1,jj ,jl) - pt_ups(ji,jj,jl) ) <= 0._wp .AND. &
1223	& pfv_ho(ji,jj,jl) * ( pt_ups(ji ,jj+1,jl) - pt_ups(ji,jj,jl) ) <= 0._wp ) THEN
1224	!
1225	IF( pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj ,jl) - zti_ups(ji,jj,jl) ) <= 0._wp .AND. &
1226	& pfv_ho(ji,jj,jl) * ( ztj_ups(ji ,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0._wp ) THEN
1227	pfu_ho(ji,jj,jl)=0._wp
1228	pfv_ho(ji,jj,jl)=0._wp
1229	ENDIF
1230	!
1231	IF( pfu_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji-1,jj ,jl) ) <= 0._wp .AND. &
1232	& pfv_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji ,jj-1,jl) ) <= 0._wp ) THEN
1233	pfu_ho(ji,jj,jl)=0._wp
1234	pfv_ho(ji,jj,jl)=0._wp
1235	ENDIF
1236	!
1237	ENDIF
1238	END DO
1239	END DO
1240	END DO
1241	CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond.
1242
1243	ENDIF
1244
1245	! Search local extrema
1246	! --------------------
1247	! max/min of pt & pt_ups with large negative/positive value (-/+zbig) outside ice cover
1248	z1_dt = 1._wp / pdt
1249	DO jl = 1, jpl
1250
1251	DO jj = 1, jpj
1252	DO ji = 1, jpi
1253	IF ( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN
1254	zbup(ji,jj) = -zbig
1255	zbdo(ji,jj) = zbig
1256	ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) > 0._wp ) THEN
1257	zbup(ji,jj) = pt_ups(ji,jj,jl)
1258	zbdo(ji,jj) = pt_ups(ji,jj,jl)
1259	ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN
1260	zbup(ji,jj) = pt(ji,jj,jl)
1261	zbdo(ji,jj) = pt(ji,jj,jl)
1262	ELSE
1263	zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_ups(ji,jj,jl) )
1264	zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_ups(ji,jj,jl) )
1265	ENDIF
1266	END DO
1267	END DO
1268
1269	DO jj = 2, jpjm1
1270	DO ji = fs_2, fs_jpim1 ! vector opt.
1271	!
1272	zup = MAX( zbup(ji,jj), zbup(ji-1,jj), zbup(ji+1,jj), zbup(ji,jj-1), zbup(ji,jj+1) ) ! search max/min in neighbourhood
1273	zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj), zbdo(ji+1,jj), zbdo(ji,jj-1), zbdo(ji,jj+1) )
1274	!
1275	zpos = MAX( 0._wp, pfu_ho(ji-1,jj ,jl) ) - MIN( 0._wp, pfu_ho(ji ,jj ,jl) ) & ! positive/negative part of the flux
1276	& + MAX( 0._wp, pfv_ho(ji ,jj-1,jl) ) - MIN( 0._wp, pfv_ho(ji ,jj ,jl) )
1277	zneg = MAX( 0._wp, pfu_ho(ji ,jj ,jl) ) - MIN( 0._wp, pfu_ho(ji-1,jj ,jl) ) &
1278	& + MAX( 0._wp, pfv_ho(ji ,jj ,jl) ) - MIN( 0._wp, pfv_ho(ji ,jj-1,jl) )
1279	!
1280	zpos = zpos - (pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) &
1281	& ) * ( 1. - pamsk )
1282	zneg = zneg + (pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) &
1283	& ) * ( 1. - pamsk )
1284	!
1285	! ! up & down beta terms
1286	! clem: zbetup and zbetdo must be 0 for zpos>1.e-10 & zneg>1.e-10 (do not put 0 instead of 1.e-10 !!!)
1287	IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt
1288	ELSE ; zbetup(ji,jj,jl) = 0._wp ! zbig
1289	ENDIF
1290	!
1291	IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
1292	ELSE ; zbetdo(ji,jj,jl) = 0._wp ! zbig
1293	ENDIF
1294	!
1295	! if all the points are outside ice cover
1296	IF( zup == -zbig ) zbetup(ji,jj,jl) = 0._wp ! zbig
1297	IF( zdo == zbig ) zbetdo(ji,jj,jl) = 0._wp ! zbig
1298	!
1299	END DO
1300	END DO
1301	END DO
1302	CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)
1303
1304
1305	! monotonic flux in the y direction
1306	! ---------------------------------
1307	DO jl = 1, jpl
1308	DO jj = 1, jpjm1
1309	DO ji = 1, fs_jpim1 ! vector opt.
1310	zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) )
1311	zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) )
1312	zcu = 0.5_wp + SIGN( 0.5_wp , pfu_ho(ji,jj,jl) )
1313	!
1314	zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu )
1315	!
1316	pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_ups(ji,jj,jl)
1317	!
1318	END DO
1319	END DO
1320
1321	DO jj = 1, jpjm1
1322	DO ji = 1, fs_jpim1 ! vector opt.
1323	zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) )
1324	zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) )
1325	zcv = 0.5_wp + SIGN( 0.5_wp , pfv_ho(ji,jj,jl) )
1326	!
1327	zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv )
1328	!
1329	pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_ups(ji,jj,jl)
1330	!
1331	END DO
1332	END DO
1333
1334	END DO
1335	!
1336	END SUBROUTINE nonosc_ice
1337
1338
1339	SUBROUTINE limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
1340	!!---------------------------------------------------------------------
1341	!! * ROUTINE limiter_x *
1342	!!
1343	!! ** Purpose : compute flux limiter
1344	!!----------------------------------------------------------------------
1345	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
1346	REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pu ! ice i-velocity => u*e2
1347	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt ! ice tracer
1348	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pfu_ups ! upstream flux
1349	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pfu_ho ! high order flux
1350	!
1351	REAL(wp) :: Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr
1352	INTEGER :: ji, jj, jl ! dummy loop indices
1353	REAL(wp), DIMENSION (jpi,jpj,jpl) :: zslpx ! tracer slopes
1354	!!----------------------------------------------------------------------
1355	!
1356	DO jl = 1, jpl
1357	DO jj = 2, jpjm1
1358	DO ji = fs_2, fs_jpim1 ! vector opt.
1359	zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1)
1360	END DO
1361	END DO
1362	END DO
1363	CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.) ! lateral boundary cond.
1364
1365	DO jl = 1, jpl
1366	DO jj = 2, jpjm1
1367	DO ji = fs_2, fs_jpim1 ! vector opt.
1368	uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj)
1369
1370	Rjm = zslpx(ji-1,jj,jl)
1371	Rj = zslpx(ji ,jj,jl)
1372	Rjp = zslpx(ji+1,jj,jl)
1373
1374	IF( kn_limiter == 3 ) THEN
1375
1376	IF( pu(ji,jj) > 0. ) THEN ; Rr = Rjm
1377	ELSE ; Rr = Rjp
1378	ENDIF
1379
1380	zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)
1381	IF( Rj > 0. ) THEN
1382	zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pu(ji,jj)), &
1383	& MIN( 2. * Rr * 0.5 * ABS(pu(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
1384	ELSE
1385	zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)), &
1386	& MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
1387	ENDIF
1388	pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter
1389
1390	ELSEIF( kn_limiter == 2 ) THEN
1391	IF( Rj /= 0. ) THEN
1392	IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj
1393	ELSE ; Cr = Rjp / Rj
1394	ENDIF
1395	ELSE
1396	Cr = 0.
1397	ENDIF
1398
1399	! -- superbee --
1400	zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
1401	! -- van albada 2 --
1402	!!zpsi = 2.Cr / (CrCr+1.)
1403	! -- sweby (with beta=1) --
1404	!!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
1405	! -- van Leer --
1406	!!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
1407	! -- ospre --
1408	!!zpsi = 1.5 * ( CrCr + Cr ) / ( CrCr + Cr + 1. )
1409	! -- koren --
1410	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( (1.+2Cr)/3., 2. ) ) )
1411	! -- charm --
1412	!IF( Cr > 0. ) THEN ; zpsi = Cr * (3.Cr + 1.) / ( (Cr + 1.) (Cr + 1.) )
1413	!ELSE ; zpsi = 0.
1414	!ENDIF
1415	! -- van albada 1 --
1416	!!zpsi = (CrCr + Cr) / (CrCr +1)
1417	! -- smart --
1418	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, 4. ) ) )
1419	! -- umist --
1420	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
1421
1422	! high order flux corrected by the limiter
1423	pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFLzpsi ) Rj * 0.5
1424
1425	ENDIF
1426	END DO
1427	END DO
1428	END DO
1429	CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.) ! lateral boundary cond.
1430	!
1431	END SUBROUTINE limiter_x
1432
1433
1434	SUBROUTINE limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
1435	!!---------------------------------------------------------------------
1436	!! * ROUTINE limiter_y *
1437	!!
1438	!! ** Purpose : compute flux limiter
1439	!!----------------------------------------------------------------------
1440	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
1441	REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pv ! ice i-velocity => u*e2
1442	REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pt ! ice tracer
1443	REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pfv_ups ! upstream flux
1444	REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfv_ho ! high order flux
1445	!
1446	REAL(wp) :: Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr
1447	INTEGER :: ji, jj, jl ! dummy loop indices
1448	REAL(wp), DIMENSION (jpi,jpj,jpl) :: zslpy ! tracer slopes
1449	!!----------------------------------------------------------------------
1450	!
1451	DO jl = 1, jpl
1452	DO jj = 2, jpjm1
1453	DO ji = fs_2, fs_jpim1 ! vector opt.
1454	zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1)
1455	END DO
1456	END DO
1457	END DO
1458	CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.) ! lateral boundary cond.
1459
1460	DO jl = 1, jpl
1461	DO jj = 2, jpjm1
1462	DO ji = fs_2, fs_jpim1 ! vector opt.
1463	vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj)
1464
1465	Rjm = zslpy(ji,jj-1,jl)
1466	Rj = zslpy(ji,jj ,jl)
1467	Rjp = zslpy(ji,jj+1,jl)
1468
1469	IF( kn_limiter == 3 ) THEN
1470
1471	IF( pv(ji,jj) > 0. ) THEN ; Rr = Rjm
1472	ELSE ; Rr = Rjp
1473	ENDIF
1474
1475	zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)
1476	IF( Rj > 0. ) THEN
1477	zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pv(ji,jj)), &
1478	& MIN( 2. * Rr * 0.5 * ABS(pv(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
1479	ELSE
1480	zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)), &
1481	& MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
1482	ENDIF
1483	pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter
1484
1485	ELSEIF( kn_limiter == 2 ) THEN
1486
1487	IF( Rj /= 0. ) THEN
1488	IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj
1489	ELSE ; Cr = Rjp / Rj
1490	ENDIF
1491	ELSE
1492	Cr = 0.
1493	ENDIF
1494
1495	! -- superbee --
1496	zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
1497	! -- van albada 2 --
1498	!!zpsi = 2.Cr / (CrCr+1.)
1499	! -- sweby (with beta=1) --
1500	!!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
1501	! -- van Leer --
1502	!!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
1503	! -- ospre --
1504	!!zpsi = 1.5 * ( CrCr + Cr ) / ( CrCr + Cr + 1. )
1505	! -- koren --
1506	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( (1.+2Cr)/3., 2. ) ) )
1507	! -- charm --
1508	!IF( Cr > 0. ) THEN ; zpsi = Cr * (3.Cr + 1.) / ( (Cr + 1.) (Cr + 1.) )
1509	!ELSE ; zpsi = 0.
1510	!ENDIF
1511	! -- van albada 1 --
1512	!!zpsi = (CrCr + Cr) / (CrCr +1)
1513	! -- smart --
1514	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, 4. ) ) )
1515	! -- umist --
1516	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
1517
1518	! high order flux corrected by the limiter
1519	pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFLzpsi ) Rj * 0.5
1520
1521	ENDIF
1522	END DO
1523	END DO
1524	END DO
1525	CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.) ! lateral boundary cond.
1526	!
1527	END SUBROUTINE limiter_y
1528
1529
1530	SUBROUTINE Hbig( phi_max, phs_max, phip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
1531	!!-------------------------------------------------------------------
1532	!! * ROUTINE Hbig *
1533	!!
1534	!! ** Purpose : Thickness correction in case advection scheme creates
1535	!! abnormally tick ice or snow
1536	!!
1537	!! ** Method : 1- check whether ice thickness is larger than the surrounding 9-points
1538	!! (before advection) and reduce it by adapting ice concentration
1539	!! 2- check whether snow thickness is larger than the surrounding 9-points
1540	!! (before advection) and reduce it by sending the excess in the ocean
1541	!! 3- check whether snow load deplets the snow-ice interface below sea level$
1542	!! and reduce it by sending the excess in the ocean
1543	!! 4- correct pond fraction to avoid a_ip > a_i
1544	!!
1545	!! ** input : Max thickness of the surrounding 9-points
1546	!!-------------------------------------------------------------------
1547	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max ! max ice thick from surrounding 9-pts
1548	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip
1549	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s
1550	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i
1551	!
1552	INTEGER :: ji, jj, jk, jl ! dummy loop indices
1553	REAL(wp) :: zhip, zhi, zhs, zvs_excess, zfra
1554	REAL(wp), DIMENSION(jpi,jpj) :: zswitch
1555	!!-------------------------------------------------------------------
1556	!
1557	!
1558	DO jl = 1, jpl
1559
1560	DO jj = 1, jpj
1561	DO ji = 1, jpi
1562	IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
1563	!
1564	! ! -- check h_ip -- !
1565	! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
1566	IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
1567	zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) )
1568	IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN
1569	pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl)
1570	ENDIF
1571	ENDIF
1572	!
1573	! ! -- check h_i -- !
1574	! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
1575	zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl)
1576	IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
1577	pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) ) !-- bound h_i to hi_max (99 m)
1578	ENDIF
1579	!
1580	! ! -- check h_s -- !
1581	! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
1582	zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl)
1583	IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
1584	zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
1585	!
1586	wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice
1587	hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
1588	!
1589	pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
1590	pv_s(ji,jj,jl) = pa_i(ji,jj,jl) * phs_max(ji,jj,jl)
1591	ENDIF
1592	!
1593	! ! -- check snow load -- !
1594	! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean
1595	! this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin)
1596	! this imposed mini can artificially make the snow very thick (if concentration decreases drastically)
1597	zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
1598	IF( zvs_excess > 0._wp ) THEN
1599	zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 )
1600	wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice
1601	hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
1602	!
1603	pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
1604	pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess
1605	ENDIF
1606
1607	ENDIF
1608	END DO
1609	END DO
1610	END DO
1611	! !-- correct pond fraction to avoid a_ip > a_i
1612	WHERE( pa_ip(:,:,:) > pa_i(:,:,:) ) pa_ip(:,:,:) = pa_i(:,:,:)
1613	!
1614	!
1615	END SUBROUTINE Hbig
1616
1617	#else
1618	!!----------------------------------------------------------------------
1619	!! Default option Dummy module NO SI3 sea-ice model
1620	!!----------------------------------------------------------------------
1621	#endif
1622
1623	!!======================================================================
1624	END MODULE icedyn_adv_umx

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

Download in other formats: