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/branches/2019/dev_r11943_MERGE_2019/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn_adv_umx.F90 @ 12252

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

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

Download in other formats: