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

Last change on this file since 10413 was 10413, checked in by clem, 5 years ago
merge dev_r9947_SI3_advection with the trunk. All sette tests passed. There is probably a conservation issue with the new advection scheme that I should solve but the routine icedyn_adv_umx.F90 is going to change anyway in the next couple of weeks. At worst, the old routine can be plugged without harm
Property svn:keywords set to `Id`
File size: 86.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 trend with the 3D advection trends using a TVD scheme
14	!! ultimate_x(_y) : compute a tracer value at velocity points using ULTIMATE scheme at various orders
15	!! macho : ???
16	!! nonosc : compute monotonic tracer fluxes by 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 lib_mpp ! MPP library
26	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
27	USE lbclnk ! lateral boundary conditions (or mpp links)
28
29	IMPLICIT NONE
30	PRIVATE
31
32	PUBLIC ice_dyn_adv_umx ! called by icedyn_adv.F90
33
34	REAL(wp) :: z1_6 = 1._wp / 6._wp ! =1/6
35	REAL(wp) :: z1_120 = 1._wp / 120._wp ! =1/120
36
37	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z1_ai, amaxu, amaxv
38
39	LOGICAL ll_dens
40
41	! advect H all the way (and get V=H*A at the end)
42	LOGICAL :: ll_thickness = .FALSE.
43
44	! look for 9 points around in nonosc limiter
45	LOGICAL :: ll_9points = .FALSE. ! false=better h?
46
47	! use HgradU in nonosc limiter
48	LOGICAL :: ll_HgradU = .TRUE. ! no effect?
49
50	! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme
51	LOGICAL :: ll_neg = .TRUE. ! keep TRUE
52
53	! limit the fluxes
54	LOGICAL :: ll_zeroup1 = .FALSE. ! false ok if Hbig otherwise needed for 2D sinon on a des valeurs de H trop fortes !!
55	LOGICAL :: ll_zeroup2 = .FALSE. ! false ok for 1D, 2D, 3D
56	LOGICAL :: ll_zeroup4 = .FALSE. ! false ok for 1D, 2D, 3D
57	LOGICAL :: ll_zeroup5 = .FALSE. ! false ok for 1D, 2D
58
59	! fluxes that are limited are uH, then (uH)*(ua/u) is used for V (only for nonosc)
60	LOGICAL :: ll_clem = .TRUE. ! simpler than rachid and works
61
62	! First advect H as H=Hdiv(u), then use H for H(n+1)=H(n)-div(uH*)
63	LOGICAL :: ll_gurvan = .FALSE. ! must be false for 1D case !!
64
65	! First guess as div(uH) (-true-) or Hdiv(u)+ugradH (-false-)
66	LOGICAL :: ll_1stguess_clem = .FALSE. ! better negative values but less good h
67
68	! advect (or not) open water. If not, retrieve it from advection of A
69	LOGICAL :: ll_ADVopw = .FALSE. !
70
71	! alternate directions for upstream
72	LOGICAL :: ll_upsxy = .TRUE.
73
74	! alternate directions for high order
75	LOGICAL :: ll_hoxy = .TRUE.
76
77	! prelimiter: use it to avoid overshoot in H
78	LOGICAL :: ll_prelimiter_zalesak = .TRUE. ! from: Zalesak(1979) eq. 14 => true is better for 1D but false is better in 3D (for h and negative values) => pb in x-y?
79	LOGICAL :: ll_prelimiter_devore = .FALSE. ! from: Devore eq. 11 (worth than zalesak)
80
81	! iterate on the limiter (only nonosc)
82	LOGICAL :: ll_limiter_it2 = .FALSE.
83
84
85	!! * Substitutions
86	# include "vectopt_loop_substitute.h90"
87	!!----------------------------------------------------------------------
88	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
89	!! $Id$
90	!! Software governed by the CeCILL licence (./LICENSE)
91	!!----------------------------------------------------------------------
92	CONTAINS
93
94	SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, &
95	& pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
96	!!----------------------------------------------------------------------
97	!! * ROUTINE ice_dyn_adv_umx *
98	!!
99	!! ** Purpose : Compute the now trend due to total advection of
100	!! tracers and add it to the general trend of tracer equations
101	!! using an "Ultimate-Macho" scheme
102	!!
103	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
104	!!----------------------------------------------------------------------
105	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
106	INTEGER , INTENT(in ) :: kt ! time step
107	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity
108	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pv_ice ! ice j-velocity
109	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area
110	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume
111	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_s ! snw volume
112	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psv_i ! salt content
113	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content
114	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration
115	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
116	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume
117	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content
118	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content
119	!
120	INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices
121	INTEGER :: icycle ! number of sub-timestep for the advection
122	REAL(wp) :: zamsk ! 1 if advection of concentration, 0 if advection of other tracers
123	REAL(wp) :: zcfl , zdt
124	REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box, zua_ho, zva_ho
125	REAL(wp), DIMENSION(jpi,jpj) :: zhvar
126	REAL(wp), DIMENSION(jpi,jpj) :: zai_b, zai_a, z1_ai_b
127	!!----------------------------------------------------------------------
128	!
129	IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
130	!
131	!
132	! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !
133	zcfl = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
134	zcfl = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
135	IF( lk_mpp ) CALL mpp_max( zcfl )
136
137	IF( zcfl > 0.5 ) THEN ; icycle = 2
138	ELSE ; icycle = 1
139	ENDIF
140
141	zdt = rdt_ice / REAL(icycle)
142
143	! --- transport --- !
144	zudy(:,:) = pu_ice(:,:) * e2u(:,:)
145	zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
146
147	! --- define velocity for advection: u*grad(H) --- !
148	DO jj = 2, jpjm1
149	DO ji = fs_2, fs_jpim1
150	IF ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN ; zcu_box(ji,jj) = 0._wp
151	ELSEIF( pu_ice(ji,jj) > 0._wp ) THEN ; zcu_box(ji,jj) = pu_ice(ji-1,jj)
152	ELSE ; zcu_box(ji,jj) = pu_ice(ji ,jj)
153	ENDIF
154
155	IF ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN ; zcv_box(ji,jj) = 0._wp
156	ELSEIF( pv_ice(ji,jj) > 0._wp ) THEN ; zcv_box(ji,jj) = pv_ice(ji,jj-1)
157	ELSE ; zcv_box(ji,jj) = pv_ice(ji,jj )
158	ENDIF
159	END DO
160	END DO
161
162	IF(.NOT. ALLOCATED(z1_ai)) ALLOCATE(z1_ai(jpi,jpj))
163	IF( ll_zeroup2 ) THEN
164	IF(.NOT. ALLOCATED(amaxu)) ALLOCATE(amaxu (jpi,jpj))
165	IF(.NOT. ALLOCATED(amaxv)) ALLOCATE(amaxv (jpi,jpj))
166	ENDIF
167	!---------------!
168	!== advection ==!
169	!---------------!
170	DO jt = 1, icycle
171
172	IF( ll_ADVopw ) THEN
173	ll_dens=.FALSE.
174	zamsk = 1._wp
175	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) ) ! Open water area
176	ELSE
177	zai_b(:,:) = SUM( pa_i(:,:,:), dim=3 )
178	ENDIF
179
180	DO jl = 1, jpl
181	!
182	WHERE( pa_i(:,:,jl) >= epsi20 ) ; z1_ai_b(:,:) = 1._wp / pa_i(:,:,jl)
183	ELSEWHERE ; z1_ai_b(:,:) = 0.
184	END WHERE
185	!
186	IF( ll_zeroup2 ) THEN
187	DO jj = 2, jpjm1
188	DO ji = fs_2, fs_jpim1
189	amaxu(ji,jj)=MAX( pa_i(ji,jj,jl), pa_i(ji,jj-1,jl), pa_i(ji,jj+1,jl), &
190	& pa_i(ji+1,jj,jl), pa_i(ji+1,jj-1,jl), pa_i(ji+1,jj+1,jl) )
191	amaxv(ji,jj)=MAX( pa_i(ji,jj,jl), pa_i(ji-1,jj,jl), pa_i(ji+1,jj,jl), &
192	& pa_i(ji,jj+1,jl), pa_i(ji-1,jj+1,jl), pa_i(ji+1,jj+1,jl) )
193	END DO
194	END DO
195	CALL lbc_lnk_multi(amaxu, 'T', 1., amaxv, 'T', 1.)
196	ENDIF
197	!
198	zamsk = 1._wp
199	ll_dens=.TRUE.
200	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl), pa_i(:,:,jl), zua_ho, zva_ho ) ! Ice area
201	ll_dens=.FALSE.
202
203	WHERE( pa_i(:,:,jl) >= epsi20 ) ; z1_ai(:,:) = 1._wp / pa_i(:,:,jl)
204	ELSEWHERE ; z1_ai(:,:) = 0.
205	END WHERE
206
207	IF( ll_thickness ) THEN
208	zua_ho(:,:) = zudy(:,:)
209	zva_ho(:,:) = zvdx(:,:)
210	ENDIF
211
212	zamsk = 0._wp ; zhvar(:,:) = pv_i (:,:,jl) * z1_ai_b(:,:)
213	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), pv_i (:,:,jl) ) ! Ice volume
214	IF( ll_thickness ) pv_i(:,:,jl) = zhvar(:,:) * pa_i(:,:,jl)
215
216	zamsk = 0._wp ; zhvar(:,:) = pv_s (:,:,jl) * z1_ai_b(:,:)
217	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), pv_s (:,:,jl) ) ! Snw volume
218	IF( ll_thickness ) pv_s(:,:,jl) = zhvar(:,:) * pa_i(:,:,jl)
219
220	zamsk = 0._wp ; zhvar(:,:) = psv_i(:,:,jl) * z1_ai_b(:,:)
221	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), psv_i(:,:,jl) ) ! Salt content
222
223	zamsk = 0._wp ; zhvar(:,:) = poa_i(:,:,jl) * z1_ai_b(:,:)
224	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), poa_i(:,:,jl) ) ! Age content
225
226	DO jk = 1, nlay_i
227	zamsk = 0._wp ; zhvar(:,:) = pe_i(:,:,jk,jl) * z1_ai_b(:,:)
228	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar(:,:), pe_i(:,:,jk,jl) ) ! Ice heat content
229	END DO
230
231	DO jk = 1, nlay_s
232	zamsk = 0._wp ; zhvar(:,:) = pe_s(:,:,jk,jl) * z1_ai_b(:,:)
233	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar(:,:), pe_s(:,:,jk,jl) ) ! Snw heat content
234	END DO
235	!
236	IF ( ln_pnd_H12 ) THEN ! melt ponds (must be the last ones to be advected because of z1_ai_b...)
237	!
238	WHERE( pa_ip(:,:,jl) >= epsi20 ) ; z1_ai_b(:,:) = 1._wp / pa_ip(:,:,jl)
239	ELSEWHERE ; z1_ai_b(:,:) = 0.
240	END WHERE
241	!
242	zamsk = 1._wp
243	ll_dens=.TRUE.
244	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl), pa_ip(:,:,jl), zua_ho, zva_ho ) ! mp fractio!n
245	ll_dens=.FALSE.
246
247	WHERE( pa_ip(:,:,jl) >= epsi20 ) ; z1_ai(:,:) = 1._wp / pa_ip(:,:,jl)
248	ELSEWHERE ; z1_ai(:,:) = 0.
249	END WHERE
250
251	zamsk = 0._wp ; zhvar(:,:) = pv_ip(:,:,jl) * z1_ai_b(:,:)
252	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), pv_ip(:,:,jl) ) ! mp volume
253	ENDIF
254	!
255	!
256	END DO
257
258	IF( .NOT. ll_ADVopw ) THEN
259	zai_a(:,:) = SUM( pa_i(:,:,:), dim=3 )
260	DO jj = 2, jpjm1
261	DO ji = fs_2, fs_jpim1
262	pato_i(ji,jj) = pato_i(ji,jj) - ( zai_a(ji,jj) - zai_b(ji,jj) ) &
263	& - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
264	END DO
265	END DO
266	CALL lbc_lnk( pato_i(:,:), 'T', 1. )
267	ENDIF
268
269	END DO
270	!
271	END SUBROUTINE ice_dyn_adv_umx
272
273
274	SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho )
275	!!----------------------------------------------------------------------
276	!! * ROUTINE adv_umx *
277	!!
278	!! ** Purpose : Compute the now trend due to total advection of
279	!! tracers and add it to the general trend of tracer equations
280	!!
281	!! ** Method : TVD scheme, i.e. 2nd order centered scheme with
282	!! corrected flux (monotonic correction)
283	!! note: - this advection scheme needs a leap-frog time scheme
284	!!
285	!! ** Action : - pt the after advective tracer
286	!!----------------------------------------------------------------------
287	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
288	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
289	INTEGER , INTENT(in ) :: jt ! number of sub-iteration
290	INTEGER , INTENT(in ) :: kt ! number of iteration
291	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
292	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pu , pv ! 2 ice velocity components => u*e2
293	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc , pvc ! 2 ice velocity components => ue2 or ua*e2u
294	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pubox, pvbox ! upstream velocity
295	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pt ! tracer field
296	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: ptc ! tracer content field
297	REAL(wp), DIMENSION(jpi,jpj), INTENT( out), OPTIONAL :: pua_ho, pva_ho ! high order u*a fluxes
298	!
299	INTEGER :: ji, jj ! dummy loop indices
300	REAL(wp) :: ztra ! local scalar
301	INTEGER :: kn_limiter = 1 ! 1=nonosc ; 2=superbee ; 3=h3(rachid)
302	REAL(wp), DIMENSION(jpi,jpj) :: zfu_ho , zfv_ho , zt_u, zt_v, zpt
303	REAL(wp), DIMENSION(jpi,jpj) :: zfu_ups, zfv_ups, zt_ups ! only for nonosc
304	!!----------------------------------------------------------------------
305	!
306	! upstream (_ups) advection with initial mass fluxes
307	! ---------------------------------------------------
308	IF( ll_clem ) zfu_ups=0.; zfv_ups=0.
309
310	IF( ll_gurvan .AND. pamsk==0. ) THEN
311	DO jj = 2, jpjm1
312	DO ji = fs_2, fs_jpim1
313	pt(ji,jj) = ( pt (ji,jj) + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) &
314	& + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
315	END DO
316	END DO
317	CALL lbc_lnk( pt, 'T', 1. )
318	ENDIF
319
320
321	IF( .NOT. ll_upsxy ) THEN
322
323	! fluxes
324	DO jj = 1, jpjm1
325	DO ji = 1, fs_jpim1
326	IF( ll_clem ) THEN
327	zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj)
328	zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1)
329	ELSE
330	zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj)
331	zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1)
332	ENDIF
333	END DO
334	END DO
335
336	ELSE
337	! 1 if advection of A
338	! z1_ai already defined IF advection of other tracers
339	IF( pamsk == 1. ) z1_ai(:,:) = 1._wp
340	!
341	IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==!
342	! flux in x-direction
343	DO jj = 1, jpjm1
344	DO ji = 1, fs_jpim1
345	IF( ll_clem ) THEN
346	zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj)
347	ELSE
348	zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj)
349	ENDIF
350	END DO
351	END DO
352
353	! first guess of tracer content from u-flux
354	DO jj = 2, jpjm1
355	DO ji = fs_2, fs_jpim1 ! vector opt.
356	IF( ll_clem ) THEN
357	IF( ll_gurvan ) THEN
358	zpt(ji,jj) = ( pt(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
359	ELSE
360	zpt(ji,jj) = ( pt(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) &
361	& + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) &
362	& * tmask(ji,jj,1)
363	ENDIF
364	ELSE
365	zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) &
366	& * tmask(ji,jj,1)
367	ENDIF
368	!! IF( ji==26 .AND. jj==86) THEN
369	!! WRITE(numout,) '***********************'
370	!! WRITE(numout,*) 'zpt upstream',zpt(ji,jj)
371	!! ENDIF
372	END DO
373	END DO
374	CALL lbc_lnk( zpt, 'T', 1. )
375	!
376	! flux in y-direction
377	DO jj = 1, jpjm1
378	DO ji = 1, fs_jpim1
379	IF( ll_clem ) THEN
380	zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1)
381	ELSE
382	zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * zpt(ji,jj+1)
383	ENDIF
384	END DO
385	END DO
386
387	!
388	ELSE !== even ice time step: adv_y then adv_x ==!
389	! flux in y-direction
390	DO jj = 1, jpjm1
391	DO ji = 1, fs_jpim1
392	IF( ll_clem ) THEN
393	zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1)
394	ELSE
395	zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1)
396	ENDIF
397	END DO
398	END DO
399
400	! first guess of tracer content from v-flux
401	DO jj = 2, jpjm1
402	DO ji = fs_2, fs_jpim1 ! vector opt.
403	IF( ll_clem ) THEN
404	IF( ll_gurvan ) THEN
405	zpt(ji,jj) = ( pt(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
406	ELSE
407	zpt(ji,jj) = ( pt(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) &
408	& + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) &
409	& * tmask(ji,jj,1)
410	ENDIF
411	ELSE
412	zpt(ji,jj) = ( ptc(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) &
413	& * tmask(ji,jj,1)
414	ENDIF
415	!! IF( ji==26 .AND. jj==86) THEN
416	!! WRITE(numout,) '***********************'
417	!! WRITE(numout,*) 'zpt upstream',zpt(ji,jj)
418	!! ENDIF
419	END DO
420	END DO
421	CALL lbc_lnk( zpt, 'T', 1. )
422	!
423	! flux in x-direction
424	DO jj = 1, jpjm1
425	DO ji = 1, fs_jpim1
426	IF( ll_clem ) THEN
427	zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj)
428	ELSE
429	zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * zpt(ji+1,jj)
430	ENDIF
431	END DO
432	END DO
433	!
434	ENDIF
435
436	ENDIF
437
438	IF( ll_clem .AND. kn_limiter /= 1 ) &
439	& CALL ctl_stop( 'STOP', 'icedyn_adv_umx: ll_clem incompatible with limiters other than nonosc' )
440
441	IF( ll_zeroup2 ) THEN
442	DO jj = 1, jpjm1
443	DO ji = 1, fs_jpim1 ! vector opt.
444	IF( amaxu(ji,jj) == 0._wp ) zfu_ups(ji,jj) = 0._wp
445	IF( amaxv(ji,jj) == 0._wp ) zfv_ups(ji,jj) = 0._wp
446	END DO
447	END DO
448	ENDIF
449
450	! guess after content field with upstream scheme
451	DO jj = 2, jpjm1
452	DO ji = fs_2, fs_jpim1
453	ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj ) &
454	& + zfv_ups(ji,jj) - zfv_ups(ji ,jj-1) ) * r1_e1e2t(ji,jj)
455	IF( ll_clem ) THEN
456	IF( ll_gurvan ) THEN
457	zt_ups(ji,jj) = ( pt (ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
458	ELSE
459	zt_ups(ji,jj) = ( pt (ji,jj) + pdt * ztra + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
460	& + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
461	ENDIF
462	ELSE
463	zt_ups(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
464	ENDIF
465	!! IF( ji==26 .AND. jj==86) THEN
466	!! WRITE(numout,) '*************************'
467	!! WRITE(numout,*) 'zt upstream',zt_ups(ji,jj)
468	!! ENDIF
469	END DO
470	END DO
471	CALL lbc_lnk( zt_ups, 'T', 1. )
472
473	! High order (_ho) fluxes
474	! -----------------------
475	SELECT CASE( kn_umx )
476	!
477	CASE ( 20 ) !== centered second order ==!
478	!
479	CALL cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho, &
480	& zt_ups, zfu_ups, zfv_ups )
481	!
482	CASE ( 1:5 ) !== 1st to 5th order ULTIMATE-MACHO scheme ==!
483	!
484	CALL macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, zt_u, zt_v, zfu_ho, zfv_ho, &
485	& zt_ups, zfu_ups, zfv_ups )
486	!
487	END SELECT
488
489	IF( ll_thickness ) THEN
490	! final trend with corrected fluxes
491	! ------------------------------------
492	DO jj = 2, jpjm1
493	DO ji = fs_2, fs_jpim1
494	IF( ll_gurvan ) THEN
495	ztra = - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj)
496	ELSE
497	ztra = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) &
498	& + ( pt(ji,jj) * ( pu(ji,jj) - pu(ji-1,jj) ) * (1.-pamsk) ) &
499	& + ( pt(ji,jj) * ( pv(ji,jj) - pv(ji,jj-1) ) * (1.-pamsk) ) ) * r1_e1e2t(ji,jj)
500	ENDIF
501	pt(ji,jj) = ( pt(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
502
503	IF( pt(ji,jj) < -epsi20 ) THEN
504	WRITE(numout,*) 'T<0 ',pt(ji,jj)
505	ENDIF
506
507	IF( pt(ji,jj) < 0._wp .AND. pt(ji,jj) >= -epsi20 ) pt(ji,jj) = 0._wp
508
509	!! IF( ji==26 .AND. jj==86) THEN
510	!! WRITE(numout,*) 'zt high order',pt(ji,jj)
511	!! ENDIF
512	END DO
513	END DO
514	CALL lbc_lnk( pt, 'T', 1. )
515	ENDIF
516
517	! Rachid trick
518	! ------------
519	IF( ll_clem ) THEN
520	IF( pamsk == 0. ) THEN
521	DO jj = 1, jpjm1
522	DO ji = 1, fs_jpim1
523	IF( ABS( puc(ji,jj) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN
524	zfu_ho (ji,jj) = zfu_ho (ji,jj) * puc(ji,jj) / pu(ji,jj)
525	zfu_ups(ji,jj) = zfu_ups(ji,jj) * puc(ji,jj) / pu(ji,jj)
526	ELSE
527	zfu_ho (ji,jj) = 0._wp
528	zfu_ups(ji,jj) = 0._wp
529	ENDIF
530	!
531	IF( ABS( pvc(ji,jj) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN
532	zfv_ho (ji,jj) = zfv_ho (ji,jj) * pvc(ji,jj) / pv(ji,jj)
533	zfv_ups(ji,jj) = zfv_ups(ji,jj) * pvc(ji,jj) / pv(ji,jj)
534	ELSE
535	zfv_ho (ji,jj) = 0._wp
536	zfv_ups(ji,jj) = 0._wp
537	ENDIF
538	ENDDO
539	ENDDO
540	ENDIF
541	ENDIF
542
543	IF( ll_zeroup5 ) THEN
544	DO jj = 2, jpjm1
545	DO ji = 2, fs_jpim1 ! vector opt.
546	zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) &
547	& - ( zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
548	IF( zpt(ji,jj) < 0. ) THEN
549	zfu_ho(ji,jj) = zfu_ups(ji,jj)
550	zfu_ho(ji-1,jj) = zfu_ups(ji-1,jj)
551	zfv_ho(ji,jj) = zfv_ups(ji,jj)
552	zfv_ho(ji,jj-1) = zfv_ups(ji,jj-1)
553	ENDIF
554	END DO
555	END DO
556	CALL lbc_lnk_multi( zfu_ho, 'U', -1., zfv_ho, 'V', -1. )
557	ENDIF
558
559	! output upstream trend of concentration and high order fluxes
560	! ------------------------------------------------------------
561	IF( ll_dens ) THEN
562	! high order u*a
563	DO jj = 1, jpjm1
564	DO ji = 1, fs_jpim1
565	pua_ho (ji,jj) = zfu_ho (ji,jj)
566	pva_ho (ji,jj) = zfv_ho (ji,jj)
567	END DO
568	END DO
569	!!CALL lbc_lnk( pua_ho, 'U', -1. ) ! clem: not needed I think
570	!!CALL lbc_lnk( pva_ho, 'V', -1. )
571	ENDIF
572
573
574	IF( .NOT.ll_thickness ) THEN
575	! final trend with corrected fluxes
576	! ------------------------------------
577	DO jj = 2, jpjm1
578	DO ji = fs_2, fs_jpim1
579	ztra = - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) & ! Div(uaH) or Div(ua)
580	& * r1_e1e2t(ji,jj) * pdt
581
582	!!IF( ptc(ji,jj)+ztra < 0._wp ) THEN
583	!! ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) + zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) & ! Div(uaH) or Div(ua)
584	!! & * r1_e1e2t(ji,jj) * pdt
585	!!ENDIF
586	!!IF( ptc(ji,jj)+ztra < 0._wp ) THEN
587	!! WRITE(numout,*) 'Tc<0 ',ptc(ji,jj)+ztra
588	!! ztra = 0._wp
589	!!ENDIF
590
591	ptc(ji,jj) = ( ptc(ji,jj) + ztra ) * tmask(ji,jj,1)
592
593	!! IF( ji==26 .AND. jj==86) THEN
594	!! WRITE(numout,*) 'ztc high order',ptc(ji,jj)
595	!! ENDIF
596
597	END DO
598	END DO
599	CALL lbc_lnk( ptc, 'T', 1. )
600	ENDIF
601
602	!
603	END SUBROUTINE adv_umx
604
605	SUBROUTINE cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, &
606	& pt_ups, pfu_ups, pfv_ups )
607	!!---------------------------------------------------------------------
608	!! * ROUTINE macho *
609	!!
610	!! ** Purpose : compute
611	!!
612	!! ** Method : ... ???
613	!! TIM = transient interpolation Modeling
614	!!
615	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
616	!!----------------------------------------------------------------------
617	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
618	INTEGER , INTENT(in ) :: kn_limiter ! limiter
619	INTEGER , INTENT(in ) :: jt ! number of sub-iteration
620	INTEGER , INTENT(in ) :: kt ! number of iteration
621	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
622	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields
623	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pu, pv ! 2 ice velocity components
624	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc, pvc ! 2 ice velocity * A components
625	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ptc ! tracer content at before time step
626	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes
627	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt_ups ! upstream guess of tracer content
628	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pfu_ups, pfv_ups ! upstream fluxes
629	!
630	INTEGER :: ji, jj ! dummy loop indices
631	LOGICAL :: ll_xy = .TRUE.
632	REAL(wp), DIMENSION(jpi,jpj) :: zzt
633	!!----------------------------------------------------------------------
634	!
635	IF( .NOT.ll_xy ) THEN !-- no alternate directions --!
636	!
637	DO jj = 1, jpjm1
638	DO ji = 1, fs_jpim1
639	pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) )
640	pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) )
641	END DO
642	END DO
643	IF ( kn_limiter == 1 ) THEN
644	IF( ll_clem ) THEN
645	CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
646	ELSE
647	CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
648	ENDIF
649	ELSEIF( kn_limiter == 2 ) THEN
650	CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
651	CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
652	ELSEIF( kn_limiter == 3 ) THEN
653	CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
654	CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
655	ENDIF
656	!
657	ELSE !-- alternate directions --!
658	!
659	IF( pamsk == 1. ) z1_ai(:,:) = 1._wp
660	!
661	IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==!
662	!
663	! flux in x-direction
664	DO jj = 1, jpjm1
665	DO ji = 1, fs_jpim1
666	pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) )
667	END DO
668	END DO
669	IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
670	IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
671
672	! first guess of tracer content from u-flux
673	DO jj = 2, jpjm1
674	DO ji = fs_2, fs_jpim1 ! vector opt.
675	IF( ll_clem ) THEN
676	zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) * z1_ai(ji,jj) ) * tmask(ji,jj,1)
677	ELSE
678	zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) * z1_ai(ji,jj)
679	ENDIF
680	END DO
681	END DO
682	CALL lbc_lnk( zzt, 'T', 1. )
683
684	! flux in y-direction
685	DO jj = 1, jpjm1
686	DO ji = 1, fs_jpim1
687	pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( zzt(ji,jj) + zzt(ji,jj+1) )
688	END DO
689	END DO
690	IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
691	IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
692
693	ELSE !== even ice time step: adv_y then adv_x ==!
694	!
695	! flux in y-direction
696	DO jj = 1, jpjm1
697	DO ji = 1, fs_jpim1
698	pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) )
699	END DO
700	END DO
701	IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
702	IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
703	!
704	! first guess of tracer content from v-flux
705	DO jj = 2, jpjm1
706	DO ji = fs_2, fs_jpim1 ! vector opt.
707	IF( ll_clem ) THEN
708	zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) * z1_ai(ji,jj) ) * tmask(ji,jj,1)
709	ELSE
710	zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) * z1_ai(ji,jj)
711	ENDIF
712	END DO
713	END DO
714	CALL lbc_lnk( zzt, 'T', 1. )
715	!
716	! flux in x-direction
717	DO jj = 1, jpjm1
718	DO ji = 1, fs_jpim1
719	pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( zzt(ji,jj) + zzt(ji+1,jj) )
720	END DO
721	END DO
722	IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
723	IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
724
725	ENDIF
726	IF( ll_clem ) THEN
727	IF( kn_limiter == 1 ) CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
728	ELSE
729	IF( kn_limiter == 1 ) CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
730	ENDIF
731
732	ENDIF
733
734	END SUBROUTINE cen2
735
736
737	SUBROUTINE macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, pt_u, pt_v, pfu_ho, pfv_ho, &
738	& pt_ups, pfu_ups, pfv_ups )
739	!!---------------------------------------------------------------------
740	!! * ROUTINE macho *
741	!!
742	!! ** Purpose : compute
743	!!
744	!! ** Method : ... ???
745	!! TIM = transient interpolation Modeling
746	!!
747	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
748	!!----------------------------------------------------------------------
749	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
750	INTEGER , INTENT(in ) :: kn_limiter ! limiter
751	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
752	INTEGER , INTENT(in ) :: jt ! number of sub-iteration
753	INTEGER , INTENT(in ) :: kt ! number of iteration
754	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
755	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields
756	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pu, pv ! 2 ice velocity components
757	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc, pvc ! 2 ice velocity * A components
758	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pubox, pvbox ! upstream velocity
759	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ptc ! tracer content at before time step
760	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_u, pt_v ! tracer at u- and v-points
761	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes
762	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt_ups ! upstream guess of tracer content
763	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pfu_ups, pfv_ups ! upstream fluxes
764	!
765	INTEGER :: ji, jj ! dummy loop indices
766	REAL(wp) :: ztra
767	REAL(wp), DIMENSION(jpi,jpj) :: zzt, zzfu_ho, zzfv_ho
768	!!----------------------------------------------------------------------
769	!
770	IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==!
771	!
772	! !-- ultimate interpolation of pt at u-point --!
773	CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
774	! !-- limiter in x --!
775	IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
776	IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
777	! !-- advective form update in zzt --!
778
779	IF( ll_1stguess_clem ) THEN
780
781	! first guess of tracer content from u-flux
782	DO jj = 2, jpjm1
783	DO ji = fs_2, fs_jpim1 ! vector opt.
784	IF( ll_clem ) THEN
785	IF( ll_gurvan ) THEN
786	zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
787	ELSE
788	zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) &
789	& + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
790	ENDIF
791	ELSE
792	zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
793	ENDIF
794	END DO
795	END DO
796	CALL lbc_lnk( zzt, 'T', 1. )
797
798	ELSE
799
800	DO jj = 2, jpjm1
801	DO ji = fs_2, fs_jpim1 ! vector opt.
802	IF( ll_gurvan ) THEN
803	zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj) &
804	& - pt (ji,jj) * pdt * ( pu (ji,jj) - pu (ji-1,jj) ) * r1_e1e2t(ji,jj)
805	ELSE
806	zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj) &
807	& - pt (ji,jj) * pdt * ( pu (ji,jj) - pu (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk
808	ENDIF
809	zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
810	END DO
811	END DO
812	CALL lbc_lnk( zzt, 'T', 1. )
813	ENDIF
814	!
815	! !-- ultimate interpolation of pt at v-point --!
816	IF( ll_hoxy ) THEN
817	CALL ultimate_y( kn_umx, pdt, zzt, pv, pvc, pt_v, pfv_ho )
818	ELSE
819	CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
820	ENDIF
821	! !-- limiter in y --!
822	IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
823	IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
824	!
825	!
826	ELSE !== even ice time step: adv_y then adv_x ==!
827	!
828	! !-- ultimate interpolation of pt at v-point --!
829	CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
830	! !-- limiter in y --!
831	IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
832	IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
833	! !-- advective form update in zzt --!
834	IF( ll_1stguess_clem ) THEN
835
836	! first guess of tracer content from v-flux
837	DO jj = 2, jpjm1
838	DO ji = fs_2, fs_jpim1 ! vector opt.
839	IF( ll_clem ) THEN
840	IF( ll_gurvan ) THEN
841	zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
842	ELSE
843	zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) &
844	& + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
845	ENDIF
846	ELSE
847	zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) &
848	& * tmask(ji,jj,1)
849	ENDIF
850	END DO
851	END DO
852	CALL lbc_lnk( zzt, 'T', 1. )
853
854	ELSE
855
856	DO jj = 2, jpjm1
857	DO ji = fs_2, fs_jpim1
858	IF( ll_gurvan ) THEN
859	zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj) &
860	& - pt (ji,jj) * pdt * ( pv (ji,jj) - pv (ji,jj-1) ) * r1_e1e2t(ji,jj)
861	ELSE
862	zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj) &
863	& - pt (ji,jj) * pdt * ( pv (ji,jj) - pv (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk
864	ENDIF
865	zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
866	END DO
867	END DO
868	CALL lbc_lnk( zzt, 'T', 1. )
869	ENDIF
870	!
871	! !-- ultimate interpolation of pt at u-point --!
872	IF( ll_hoxy ) THEN
873	CALL ultimate_x( kn_umx, pdt, zzt, pu, puc, pt_u, pfu_ho )
874	ELSE
875	CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
876	ENDIF
877	! !-- limiter in x --!
878	IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
879	IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
880	!
881	!
882	ENDIF
883
884
885	IF( kn_limiter == 1 ) THEN
886	IF( .NOT. ll_limiter_it2 ) THEN
887	IF( ll_clem ) THEN
888	CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
889	ELSE
890	CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
891	ENDIF
892	ELSE
893	zzfu_ho(:,:) = pfu_ho(:,:)
894	zzfv_ho(:,:) = pfv_ho(:,:)
895	! 1st iteration of nonosc (limit the flux with the upstream solution)
896	IF( ll_clem ) THEN
897	CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho )
898	ELSE
899	CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho )
900	ENDIF
901	! guess after content field with high order
902	DO jj = 2, jpjm1
903	DO ji = fs_2, fs_jpim1
904	ztra = - ( zzfu_ho(ji,jj) - zzfu_ho(ji-1,jj) + zzfv_ho(ji,jj) - zzfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj)
905	zzt(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
906	END DO
907	END DO
908	CALL lbc_lnk( zzt, 'T', 1. )
909	! 2nd iteration of nonosc (limit the flux with the limited high order solution)
910	IF( ll_clem ) THEN
911	CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho )
912	ELSE
913	CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho )
914	ENDIF
915	ENDIF
916	ENDIF
917	!
918	END SUBROUTINE macho
919
920
921	SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
922	!!---------------------------------------------------------------------
923	!! * ROUTINE ultimate_x *
924	!!
925	!! ** Purpose : compute
926	!!
927	!! ** Method : ... ???
928	!! TIM = transient interpolation Modeling
929	!!
930	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
931	!!----------------------------------------------------------------------
932	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
933	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
934	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pu ! ice i-velocity component
935	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc ! ice i-velocity * A component
936	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields
937	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_u ! tracer at u-point
938	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pfu_ho ! high order flux
939	!
940	INTEGER :: ji, jj ! dummy loop indices
941	REAL(wp) :: zcu, zdx2, zdx4 ! - -
942	REAL(wp), DIMENSION(jpi,jpj) :: ztu1, ztu2, ztu3, ztu4
943	!!----------------------------------------------------------------------
944	!
945	! !-- Laplacian in i-direction --!
946	DO jj = 2, jpjm1 ! First derivative (gradient)
947	DO ji = 1, fs_jpim1
948	ztu1(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
949	END DO
950	! ! Second derivative (Laplacian)
951	DO ji = fs_2, fs_jpim1
952	ztu2(ji,jj) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj)
953	END DO
954	END DO
955	CALL lbc_lnk( ztu2, 'T', 1. )
956	!
957	! !-- BiLaplacian in i-direction --!
958	DO jj = 2, jpjm1 ! Third derivative
959	DO ji = 1, fs_jpim1
960	ztu3(ji,jj) = ( ztu2(ji+1,jj) - ztu2(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
961	END DO
962	! ! Fourth derivative
963	DO ji = fs_2, fs_jpim1
964	ztu4(ji,jj) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj)
965	END DO
966	END DO
967	CALL lbc_lnk( ztu4, 'T', 1. )
968	!
969	!
970	SELECT CASE (kn_umx )
971	!
972	CASE( 1 ) !== 1st order central TIM ==! (Eq. 21)
973	!
974	DO jj = 1, jpjm1
975	DO ji = 1, fs_jpim1 ! vector opt.
976	pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj) + pt(ji,jj) &
977	& - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) )
978	END DO
979	END DO
980	!
981	CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23)
982	!
983	DO jj = 1, jpjm1
984	DO ji = 1, fs_jpim1 ! vector opt.
985	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
986	pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj) + pt(ji,jj) &
987	& - zcu * ( pt(ji+1,jj) - pt(ji,jj) ) )
988	END DO
989	END DO
990	!
991	CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24)
992	!
993	DO jj = 1, jpjm1
994	DO ji = 1, fs_jpim1 ! vector opt.
995	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
996	zdx2 = e1u(ji,jj) * e1u(ji,jj)
997	!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
998	pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj) + pt (ji,jj) &
999	& - zcu * ( pt (ji+1,jj) - pt (ji,jj) ) ) &
1000	& + z1_6 * zdx2 * ( zcuzcu - 1._wp ) ( ztu2(ji+1,jj) + ztu2(ji,jj) &
1001	& - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) ) )
1002	END DO
1003	END DO
1004	!
1005	CASE( 4 ) !== 4th order central TIM ==! (Eq. 27)
1006	!
1007	DO jj = 1, jpjm1
1008	DO ji = 1, fs_jpim1 ! vector opt.
1009	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
1010	zdx2 = e1u(ji,jj) * e1u(ji,jj)
1011	!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
1012	pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj) + pt (ji,jj) &
1013	& - zcu * ( pt (ji+1,jj) - pt (ji,jj) ) ) &
1014	& + z1_6 * zdx2 * ( zcuzcu - 1._wp ) ( ztu2(ji+1,jj) + ztu2(ji,jj) &
1015	& - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) ) )
1016	END DO
1017	END DO
1018	!
1019	CASE( 5 ) !== 5th order central TIM ==! (Eq. 29)
1020	!
1021	DO jj = 1, jpjm1
1022	DO ji = 1, fs_jpim1 ! vector opt.
1023	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
1024	zdx2 = e1u(ji,jj) * e1u(ji,jj)
1025	!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
1026	zdx4 = zdx2 * zdx2
1027	pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj) + pt (ji,jj) &
1028	& - zcu * ( pt (ji+1,jj) - pt (ji,jj) ) ) &
1029	& + z1_6 * zdx2 * ( zcuzcu - 1._wp ) ( ztu2(ji+1,jj) + ztu2(ji,jj) &
1030	& - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) ) &
1031	& + z1_120 * zdx4 * ( zcuzcu - 1._wp ) ( zcuzcu - 4._wp ) ( ztu4(ji+1,jj) + ztu4(ji,jj) &
1032	& - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj) - ztu4(ji,jj) ) ) )
1033	END DO
1034	END DO
1035	!
1036	END SELECT
1037	! !-- High order flux in i-direction --!
1038	IF( ll_neg ) THEN
1039	DO jj = 1, jpjm1
1040	DO ji = 1, fs_jpim1
1041	IF( pt_u(ji,jj) < 0._wp ) THEN
1042	pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj) + pt(ji,jj) &
1043	& - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) )
1044	ENDIF
1045	END DO
1046	END DO
1047	ENDIF
1048
1049	DO jj = 1, jpjm1
1050	DO ji = 1, fs_jpim1 ! vector opt.
1051	IF( ll_clem ) THEN
1052	pfu_ho(ji,jj) = pu(ji,jj) * pt_u(ji,jj)
1053	ELSE
1054	pfu_ho(ji,jj) = puc(ji,jj) * pt_u(ji,jj)
1055	ENDIF
1056	END DO
1057	END DO
1058	!
1059	END SUBROUTINE ultimate_x
1060
1061
1062	SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
1063	!!---------------------------------------------------------------------
1064	!! * ROUTINE ultimate_y *
1065	!!
1066	!! ** Purpose : compute
1067	!!
1068	!! ** Method : ... ???
1069	!! TIM = transient interpolation Modeling
1070	!!
1071	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
1072	!!----------------------------------------------------------------------
1073	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
1074	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
1075	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pv ! ice j-velocity component
1076	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pvc ! ice j-velocity*A component
1077	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields
1078	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_v ! tracer at v-point
1079	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pfv_ho ! high order flux
1080	!
1081	INTEGER :: ji, jj ! dummy loop indices
1082	REAL(wp) :: zcv, zdy2, zdy4 ! - -
1083	REAL(wp), DIMENSION(jpi,jpj) :: ztv1, ztv2, ztv3, ztv4
1084	!!----------------------------------------------------------------------
1085	!
1086	! !-- Laplacian in j-direction --!
1087	DO jj = 1, jpjm1 ! First derivative (gradient)
1088	DO ji = fs_2, fs_jpim1
1089	ztv1(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
1090	END DO
1091	END DO
1092	DO jj = 2, jpjm1 ! Second derivative (Laplacian)
1093	DO ji = fs_2, fs_jpim1
1094	ztv2(ji,jj) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj)
1095	END DO
1096	END DO
1097	CALL lbc_lnk( ztv2, 'T', 1. )
1098	!
1099	! !-- BiLaplacian in j-direction --!
1100	DO jj = 1, jpjm1 ! First derivative
1101	DO ji = fs_2, fs_jpim1
1102	ztv3(ji,jj) = ( ztv2(ji,jj+1) - ztv2(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
1103	END DO
1104	END DO
1105	DO jj = 2, jpjm1 ! Second derivative
1106	DO ji = fs_2, fs_jpim1
1107	ztv4(ji,jj) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj)
1108	END DO
1109	END DO
1110	CALL lbc_lnk( ztv4, 'T', 1. )
1111	!
1112	!
1113	SELECT CASE (kn_umx )
1114	!
1115	CASE( 1 ) !== 1st order central TIM ==! (Eq. 21)
1116	DO jj = 1, jpjm1
1117	DO ji = 1, fs_jpim1
1118	pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1) + pt(ji,jj) ) &
1119	& - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) )
1120	END DO
1121	END DO
1122	!
1123	CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23)
1124	DO jj = 1, jpjm1
1125	DO ji = 1, fs_jpim1
1126	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1127	pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1) + pt(ji,jj) ) &
1128	& - zcv * ( pt(ji,jj+1) - pt(ji,jj) ) )
1129	END DO
1130	END DO
1131	CALL lbc_lnk( pt_v, 'V', 1. )
1132	!
1133	CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24)
1134	DO jj = 1, jpjm1
1135	DO ji = 1, fs_jpim1
1136	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1137	zdy2 = e2v(ji,jj) * e2v(ji,jj)
1138	!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
1139	pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1) + pt (ji,jj) &
1140	& - zcv * ( pt (ji,jj+1) - pt (ji,jj) ) ) &
1141	& + z1_6 * zdy2 * ( zcvzcv - 1._wp ) ( ztv2(ji,jj+1) + ztv2(ji,jj) &
1142	& - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) )
1143	END DO
1144	END DO
1145	!
1146	CASE( 4 ) !== 4th order central TIM ==! (Eq. 27)
1147	DO jj = 1, jpjm1
1148	DO ji = 1, fs_jpim1
1149	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1150	zdy2 = e2v(ji,jj) * e2v(ji,jj)
1151	!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
1152	pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1) + pt (ji,jj) &
1153	& - zcv * ( pt (ji,jj+1) - pt (ji,jj) ) ) &
1154	& + z1_6 * zdy2 * ( zcvzcv - 1._wp ) ( ztv2(ji,jj+1) + ztv2(ji,jj) &
1155	& - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) )
1156	END DO
1157	END DO
1158	!
1159	CASE( 5 ) !== 5th order central TIM ==! (Eq. 29)
1160	DO jj = 1, jpjm1
1161	DO ji = 1, fs_jpim1
1162	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1163	zdy2 = e2v(ji,jj) * e2v(ji,jj)
1164	!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
1165	zdy4 = zdy2 * zdy2
1166	pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1) + pt (ji,jj) &
1167	& - zcv * ( pt (ji,jj+1) - pt (ji,jj) ) ) &
1168	& + z1_6 * zdy2 * ( zcvzcv - 1._wp ) ( ztv2(ji,jj+1) + ztv2(ji,jj) &
1169	& - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) &
1170	& + z1_120 * zdy4 * ( zcvzcv - 1._wp ) ( zcvzcv - 4._wp ) ( ztv4(ji,jj+1) + ztv4(ji,jj) &
1171	& - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1) - ztv4(ji,jj) ) ) )
1172	END DO
1173	END DO
1174	!
1175	END SELECT
1176	! !-- High order flux in j-direction --!
1177	IF( ll_neg ) THEN
1178	DO jj = 1, jpjm1
1179	DO ji = 1, fs_jpim1
1180	IF( pt_v(ji,jj) < 0._wp ) THEN
1181	pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1) + pt(ji,jj) ) &
1182	& - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) )
1183	ENDIF
1184	END DO
1185	END DO
1186	ENDIF
1187
1188	DO jj = 1, jpjm1
1189	DO ji = 1, fs_jpim1 ! vector opt.
1190	IF( ll_clem ) THEN
1191	pfv_ho(ji,jj) = pv(ji,jj) * pt_v(ji,jj)
1192	ELSE
1193	pfv_ho(ji,jj) = pvc(ji,jj) * pt_v(ji,jj)
1194	ENDIF
1195	END DO
1196	END DO
1197	!
1198	END SUBROUTINE ultimate_y
1199
1200
1201	SUBROUTINE nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho )
1202	!!---------------------------------------------------------------------
1203	!! * ROUTINE nonosc *
1204	!!
1205	!! ** Purpose : compute monotonic tracer fluxes from the upstream
1206	!! scheme and the before field by a nonoscillatory algorithm
1207	!!
1208	!! ** Method : ... ???
1209	!! warning : pt and pt_low must be masked, but the boundaries
1210	!! conditions on the fluxes are not necessary zalezak (1979)
1211	!! drange (1995) multi-dimensional forward-in-time and upstream-
1212	!! in-space based differencing for fluid
1213	!!----------------------------------------------------------------------
1214	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
1215	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
1216	REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pu ! ice i-velocity => u*e2
1217	REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: puc ! ice i-velocity A => ue2*a
1218	REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pv ! ice j-velocity => v*e1
1219	REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pvc ! ice j-velocity A => ve1*a
1220	REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: ptc, pt, pt_low ! before field & upstream guess of after field
1221	REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) :: pfv_low, pfu_low ! upstream flux
1222	REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) :: pfv_ho, pfu_ho ! monotonic flux
1223	!
1224	INTEGER :: ji, jj ! dummy loop indices
1225	REAL(wp) :: zpos, zneg, zbig, zsml, z1_dt, zpos2, zneg2 ! local scalars
1226	REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign, zcoef ! - -
1227	REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zti_low, ztj_low, zzt
1228	!!----------------------------------------------------------------------
1229	zbig = 1.e+40_wp
1230	zsml = epsi20
1231
1232	IF( ll_zeroup2 ) THEN
1233	DO jj = 1, jpjm1
1234	DO ji = 1, fs_jpim1 ! vector opt.
1235	IF( amaxu(ji,jj) == 0._wp ) pfu_ho(ji,jj) = 0._wp
1236	IF( amaxv(ji,jj) == 0._wp ) pfv_ho(ji,jj) = 0._wp
1237	END DO
1238	END DO
1239	ENDIF
1240
1241	IF( ll_zeroup4 ) THEN
1242	DO jj = 1, jpjm1
1243	DO ji = 1, fs_jpim1 ! vector opt.
1244	IF( pfu_low(ji,jj) == 0._wp ) pfu_ho(ji,jj) = 0._wp
1245	IF( pfv_low(ji,jj) == 0._wp ) pfv_ho(ji,jj) = 0._wp
1246	END DO
1247	END DO
1248	ENDIF
1249
1250
1251	IF( ll_zeroup1 ) THEN
1252	DO jj = 2, jpjm1
1253	DO ji = fs_2, fs_jpim1
1254	IF( ll_gurvan ) THEN
1255	zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) &
1256	& - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
1257	ELSE
1258	zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) &
1259	& - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) &
1260	& + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1261	& + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
1262	ENDIF
1263	IF( zzt(ji,jj) < 0._wp ) THEN
1264	pfu_ho(ji,jj) = pfu_low(ji,jj)
1265	pfv_ho(ji,jj) = pfv_low(ji,jj)
1266	WRITE(numout,) ' 1 negative high order zzt *',ji,jj,zzt(ji,jj)
1267	ENDIF
1268	!! IF( ji==26 .AND. jj==86) THEN
1269	!! WRITE(numout,*) 'zzt high order',zzt(ji,jj)
1270	!! WRITE(numout,) 'pfu_ho',(pfu_ho(ji,jj)) r1_e1e2t(ji,jj) * pdt
1271	!! WRITE(numout,) 'pfv_ho',(pfv_ho(ji,jj)) r1_e1e2t(ji,jj) * pdt
1272	!! WRITE(numout,) 'pfu_hom1',(pfu_ho(ji-1,jj)) r1_e1e2t(ji,jj) * pdt
1273	!! WRITE(numout,) 'pfv_hom1',(pfv_ho(ji,jj-1)) r1_e1e2t(ji,jj) * pdt
1274	!! ENDIF
1275	IF( ll_gurvan ) THEN
1276	zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) &
1277	& - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
1278	ELSE
1279	zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) &
1280	& - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) &
1281	& + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1282	& + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
1283	ENDIF
1284	IF( zzt(ji,jj) < 0._wp ) THEN
1285	pfu_ho(ji-1,jj) = pfu_low(ji-1,jj)
1286	pfv_ho(ji,jj-1) = pfv_low(ji,jj-1)
1287	WRITE(numout,) ' 2 negative high order zzt *',ji,jj,zzt(ji,jj)
1288	ENDIF
1289	IF( ll_gurvan ) THEN
1290	zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) &
1291	& - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
1292	ELSE
1293	zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) &
1294	& - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) &
1295	& + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1296	& + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
1297	ENDIF
1298	IF( zzt(ji,jj) < 0._wp ) THEN
1299	WRITE(numout,) ' 3 negative high order zzt *',ji,jj,zzt(ji,jj)
1300	ENDIF
1301	END DO
1302	END DO
1303	CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. )
1304	ENDIF
1305
1306
1307	! antidiffusive flux : high order minus low order
1308	! --------------------------------------------------
1309	DO jj = 1, jpjm1
1310	DO ji = 1, fs_jpim1 ! vector opt.
1311	pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_low(ji,jj)
1312	pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_low(ji,jj)
1313	END DO
1314	END DO
1315
1316	! extreme case where pfu_ho has to be zero
1317	! ----------------------------------------
1318	! pfu_ho
1319	! * --->
1320	! \| \| * \| \|
1321	! \| \| \| * \|
1322	! \| \| \| \| *
1323	! t_low : i-1 i i+1 i+2
1324	IF( ll_prelimiter_zalesak ) THEN
1325
1326	DO jj = 2, jpjm1
1327	DO ji = fs_2, fs_jpim1
1328	zti_low(ji,jj)= pt_low(ji+1,jj )
1329	ztj_low(ji,jj)= pt_low(ji ,jj+1)
1330	END DO
1331	END DO
1332	CALL lbc_lnk_multi( zti_low, 'T', 1., ztj_low, 'T', 1. )
1333
1334
1335	!! this does not work
1336	!! DO jj = 2, jpjm1
1337	!! DO ji = fs_2, fs_jpim1
1338	!! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji+1,jj ) - pt_low (ji ,jj) ) .AND. &
1339	!! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji ,jj+1) - pt_low (ji ,jj) ) &
1340	!! & ) THEN
1341	!! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_low(ji+1,jj ) - zti_low(ji ,jj) ) .AND. &
1342	!! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_low(ji,jj+1 ) - ztj_low(ji ,jj) ) &
1343	!! & ) THEN
1344	!! pfu_ho(ji,jj) = 0. ; pfv_ho(ji,jj) = 0.
1345	!! ENDIF
1346	!! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji ,jj) - pt_low (ji-1,jj ) ) .AND. &
1347	!! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji ,jj) - pt_low (ji ,jj-1) ) &
1348	!! & ) THEN
1349	!! pfu_ho(ji,jj) = 0. ; pfv_ho(ji,jj) = 0.
1350	!! ENDIF
1351	!! ENDIF
1352	!! END DO
1353	!! END DO
1354
1355	DO jj = 2, jpjm1
1356	DO ji = fs_2, fs_jpim1
1357	IF ( pfu_ho(ji,jj) * ( pt_low(ji+1,jj) - pt_low(ji,jj) ) <= 0. .AND. &
1358	& pfv_ho(ji,jj) * ( pt_low(ji,jj+1) - pt_low(ji,jj) ) <= 0. ) THEN
1359	!
1360	IF( pfu_ho(ji,jj) * ( zti_low(ji+1,jj) - zti_low(ji,jj) ) <= 0 .AND. &
1361	& pfv_ho(ji,jj) * ( ztj_low(ji,jj+1) - ztj_low(ji,jj) ) <= 0) pfu_ho(ji,jj)=0. ; pfv_ho(ji,jj)=0.
1362	!
1363	IF( pfu_ho(ji,jj) * ( pt_low(ji ,jj) - pt_low(ji-1,jj) ) <= 0 .AND. &
1364	& pfv_ho(ji,jj) * ( pt_low(ji ,jj) - pt_low(ji,jj-1) ) <= 0) pfu_ho(ji,jj)=0. ; pfv_ho(ji,jj)=0.
1365	!
1366	ENDIF
1367	END DO
1368	END DO
1369	CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond.
1370
1371	ELSEIF( ll_prelimiter_devore ) THEN
1372	DO jj = 2, jpjm1
1373	DO ji = fs_2, fs_jpim1
1374	zti_low(ji,jj)= pt_low(ji+1,jj )
1375	ztj_low(ji,jj)= pt_low(ji ,jj+1)
1376	END DO
1377	END DO
1378	CALL lbc_lnk_multi( zti_low, 'T', 1., ztj_low, 'T', 1. )
1379
1380	z1_dt = 1._wp / pdt
1381	DO jj = 2, jpjm1
1382	DO ji = fs_2, fs_jpim1
1383	zsign = SIGN( 1., pt_low(ji+1,jj) - pt_low(ji,jj) )
1384	pfu_ho(ji,jj) = zsign * MAX( 0. , MIN( ABS(pfu_ho(ji,jj)) , &
1385	& zsign * ( pt_low (ji ,jj) - pt_low (ji-1,jj) ) * e1e2t(ji ,jj) * z1_dt , &
1386	& zsign * ( zti_low(ji+1,jj) - zti_low(ji ,jj) ) * e1e2t(ji+1,jj) * z1_dt ) )
1387
1388	zsign = SIGN( 1., pt_low(ji,jj+1) - pt_low(ji,jj) )
1389	pfv_ho(ji,jj) = zsign * MAX( 0. , MIN( ABS(pfv_ho(ji,jj)) , &
1390	& zsign * ( pt_low (ji,jj ) - pt_low (ji,jj-1) ) * e1e2t(ji,jj ) * z1_dt , &
1391	& zsign * ( ztj_low(ji,jj+1) - ztj_low(ji,jj ) ) * e1e2t(ji,jj+1) * z1_dt ) )
1392	END DO
1393	END DO
1394	CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond.
1395
1396	ENDIF
1397
1398
1399	! Search local extrema
1400	! --------------------
1401	! max/min of pt & pt_low with large negative/positive value (-/+zbig) outside ice cover
1402	DO jj = 1, jpj
1403	DO ji = 1, jpi
1404	IF ( pt(ji,jj) <= 0._wp .AND. pt_low(ji,jj) <= 0._wp ) THEN
1405	zbup(ji,jj) = -zbig
1406	zbdo(ji,jj) = zbig
1407	ELSEIF( pt(ji,jj) <= 0._wp .AND. pt_low(ji,jj) > 0._wp ) THEN
1408	zbup(ji,jj) = pt_low(ji,jj)
1409	zbdo(ji,jj) = pt_low(ji,jj)
1410	ELSEIF( pt(ji,jj) > 0._wp .AND. pt_low(ji,jj) <= 0._wp ) THEN
1411	zbup(ji,jj) = pt(ji,jj)
1412	zbdo(ji,jj) = pt(ji,jj)
1413	ELSE
1414	zbup(ji,jj) = MAX( pt(ji,jj) , pt_low(ji,jj) )
1415	zbdo(ji,jj) = MIN( pt(ji,jj) , pt_low(ji,jj) )
1416	ENDIF
1417	END DO
1418	END DO
1419
1420
1421	z1_dt = 1._wp / pdt
1422	DO jj = 2, jpjm1
1423	DO ji = fs_2, fs_jpim1 ! vector opt.
1424	!
1425	IF( .NOT. ll_9points ) THEN
1426	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
1427	zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj ), zbdo(ji+1,jj ), zbdo(ji ,jj-1), zbdo(ji ,jj+1) )
1428	!
1429	ELSE
1430	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
1431	& zbup(ji-1,jj-1), zbup(ji+1,jj+1), zbup(ji+1,jj-1), zbup(ji-1,jj+1) )
1432	zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj ), zbdo(ji+1,jj ), zbdo(ji ,jj-1), zbdo(ji ,jj+1), &
1433	& zbdo(ji-1,jj-1), zbdo(ji+1,jj+1), zbdo(ji+1,jj-1), zbdo(ji-1,jj+1) )
1434	ENDIF
1435	!
1436	zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji ,jj) ) & ! positive/negative part of the flux
1437	& + MAX( 0., pfv_ho(ji,jj-1) ) - MIN( 0., pfv_ho(ji,jj ) )
1438	zneg = MAX( 0., pfu_ho(ji ,jj) ) - MIN( 0., pfu_ho(ji-1,jj) ) &
1439	& + MAX( 0., pfv_ho(ji,jj ) ) - MIN( 0., pfv_ho(ji,jj-1) )
1440	!
1441	IF( ll_HgradU .AND. .NOT.ll_gurvan ) THEN
1442	zneg2 = ( pt(ji,jj) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) ) * ( 1. - pamsk )
1443	zpos2 = ( - pt(ji,jj) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) - pt(ji,jj) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) ) * ( 1. - pamsk )
1444	ELSE
1445	zneg2 = 0. ; zpos2 = 0.
1446	ENDIF
1447	!
1448	! ! up & down beta terms
1449	! zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / ( zpos + zsml ) * e1e2t(ji,jj) * z1_dt
1450	! zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / ( zneg + zsml ) * e1e2t(ji,jj) * z1_dt
1451
1452	IF( (zpos+zpos2) > 0. ) THEN ; zbetup(ji,jj) = MAX( 0._wp, zup - pt_low(ji,jj) ) / (zpos+zpos2) * e1e2t(ji,jj) * z1_dt
1453	ELSE ; zbetup(ji,jj) = 0. ! zbig
1454	ENDIF
1455	!
1456	IF( (zneg+zneg2) > 0. ) THEN ; zbetdo(ji,jj) = MAX( 0._wp, pt_low(ji,jj) - zdo ) / (zneg+zneg2) * e1e2t(ji,jj) * z1_dt
1457	ELSE ; zbetdo(ji,jj) = 0. ! zbig
1458	ENDIF
1459	!
1460	! if all the points are outside ice cover
1461	IF( zup == -zbig ) zbetup(ji,jj) = 0. ! zbig
1462	IF( zdo == zbig ) zbetdo(ji,jj) = 0. ! zbig
1463	!
1464
1465	!! IF( ji==26 .AND. jj==86) THEN
1466	! WRITE(numout,*) '-----------------'
1467	! WRITE(numout,*) 'zpos',zpos,zpos2
1468	! WRITE(numout,*) 'zneg',zneg,zneg2
1469	! WRITE(numout,*) 'puc/pu',ABS(puc(ji,jj))/MAX(epsi20, ABS(pu(ji,jj)))
1470	! WRITE(numout,*) 'pvc/pv',ABS(pvc(ji,jj))/MAX(epsi20, ABS(pv(ji,jj)))
1471	! WRITE(numout,*) 'pucm1/pu',ABS(puc(ji-1,jj))/MAX(epsi20, ABS(pu(ji-1,jj)))
1472	! WRITE(numout,*) 'pvcm1/pv',ABS(pvc(ji,jj-1))/MAX(epsi20, ABS(pv(ji,jj-1)))
1473	! WRITE(numout,) 'pfu_ho',(pfu_ho(ji,jj)+pfu_low(ji,jj)) r1_e1e2t(ji,jj) * pdt
1474	! WRITE(numout,) 'pfv_ho',(pfv_ho(ji,jj)+pfv_low(ji,jj)) r1_e1e2t(ji,jj) * pdt
1475	! WRITE(numout,) 'pfu_hom1',(pfu_ho(ji-1,jj)+pfu_low(ji-1,jj)) r1_e1e2t(ji,jj) * pdt
1476	! WRITE(numout,) 'pfv_hom1',(pfv_ho(ji,jj-1)+pfv_low(ji,jj-1)) r1_e1e2t(ji,jj) * pdt
1477	! WRITE(numout,) 'pfu_low',pfu_low(ji,jj) r1_e1e2t(ji,jj) * pdt
1478	! WRITE(numout,) 'pfv_low',pfv_low(ji,jj) r1_e1e2t(ji,jj) * pdt
1479	! WRITE(numout,) 'pfu_lowm1',pfu_low(ji-1,jj) r1_e1e2t(ji,jj) * pdt
1480	! WRITE(numout,) 'pfv_lowm1',pfv_low(ji,jj-1) r1_e1e2t(ji,jj) * pdt
1481	!
1482	! WRITE(numout,*) 'pt',pt(ji,jj)
1483	! WRITE(numout,*) 'ptim1',pt(ji-1,jj)
1484	! WRITE(numout,*) 'ptjm1',pt(ji,jj-1)
1485	! WRITE(numout,*) 'pt_low',pt_low(ji,jj)
1486	! WRITE(numout,*) 'zbetup',zbetup(ji,jj)
1487	! WRITE(numout,*) 'zbetdo',zbetdo(ji,jj)
1488	! WRITE(numout,*) 'zup',zup
1489	! WRITE(numout,*) 'zdo',zdo
1490	! ENDIF
1491	!
1492	END DO
1493	END DO
1494	CALL lbc_lnk_multi( zbetup, 'T', 1., zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)
1495
1496
1497	! monotonic flux in the y direction
1498	! ---------------------------------
1499	DO jj = 1, jpjm1
1500	DO ji = 1, fs_jpim1 ! vector opt.
1501	zau = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji+1,jj) )
1502	zbu = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji+1,jj) )
1503	zcu = 0.5 + SIGN( 0.5 , pfu_ho(ji,jj) )
1504	!
1505	zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu )
1506
1507	pfu_ho(ji,jj) = pfu_ho(ji,jj) * zcoef + pfu_low(ji,jj)
1508
1509	!! IF( ji==26 .AND. jj==86) THEN
1510	!! WRITE(numout,*) 'coefU',zcoef
1511	!! WRITE(numout,) 'pfu_ho',(pfu_ho(ji,jj)) r1_e1e2t(ji,jj) * pdt
1512	!! WRITE(numout,) 'pfu_hom1',(pfu_ho(ji-1,jj)) r1_e1e2t(ji,jj) * pdt
1513	!! ENDIF
1514
1515	END DO
1516	END DO
1517
1518	DO jj = 1, jpjm1
1519	DO ji = 1, fs_jpim1 ! vector opt.
1520	zav = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji,jj+1) )
1521	zbv = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji,jj+1) )
1522	zcv = 0.5 + SIGN( 0.5 , pfv_ho(ji,jj) )
1523	!
1524	zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv )
1525
1526	pfv_ho(ji,jj) = pfv_ho(ji,jj) * zcoef + pfv_low(ji,jj)
1527
1528	!! IF( ji==26 .AND. jj==86) THEN
1529	!! WRITE(numout,*) 'coefV',zcoef
1530	!! WRITE(numout,) 'pfv_ho',(pfv_ho(ji,jj)) r1_e1e2t(ji,jj) * pdt
1531	!! WRITE(numout,) 'pfv_hom1',(pfv_ho(ji,jj-1)) r1_e1e2t(ji,jj) * pdt
1532	!! ENDIF
1533	END DO
1534	END DO
1535
1536	! clem test
1537	DO jj = 2, jpjm1
1538	DO ji = 2, fs_jpim1 ! vector opt.
1539	IF( ll_gurvan ) THEN
1540	zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) &
1541	& - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
1542	ELSE
1543	zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) &
1544	& - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) &
1545	& + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1546	& + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
1547	ENDIF
1548	IF( zzt(ji,jj) < -epsi20 ) THEN
1549	WRITE(numout,*) 'T<0 nonosc',zzt(ji,jj)
1550	ENDIF
1551	END DO
1552	END DO
1553
1554	!
1555	!
1556	END SUBROUTINE nonosc_2d
1557
1558	SUBROUTINE limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
1559	!!---------------------------------------------------------------------
1560	!! * ROUTINE limiter_x *
1561	!!
1562	!! ** Purpose : compute flux limiter
1563	!!----------------------------------------------------------------------
1564	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
1565	REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pu ! ice i-velocity => u*e2
1566	REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: puc ! ice i-velocity A => ue2*a
1567	REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pt ! ice tracer
1568	REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) :: pfu_ho ! high order flux
1569	REAL(wp), DIMENSION (jpi,jpj), INTENT(in ), OPTIONAL :: pfu_ups ! upstream flux
1570	!
1571	REAL(wp) :: Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr
1572	INTEGER :: ji, jj ! dummy loop indices
1573	REAL(wp), DIMENSION (jpi,jpj) :: zslpx ! tracer slopes
1574	!!----------------------------------------------------------------------
1575	!
1576	DO jj = 2, jpjm1
1577	DO ji = fs_2, fs_jpim1 ! vector opt.
1578	zslpx(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * umask(ji,jj,1)
1579	END DO
1580	END DO
1581	CALL lbc_lnk( zslpx, 'U', -1.) ! lateral boundary cond.
1582
1583	DO jj = 2, jpjm1
1584	DO ji = fs_2, fs_jpim1 ! vector opt.
1585	uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj)
1586
1587	Rjm = zslpx(ji-1,jj)
1588	Rj = zslpx(ji ,jj)
1589	Rjp = zslpx(ji+1,jj)
1590
1591	IF( PRESENT(pfu_ups) ) THEN
1592
1593	IF( pu(ji,jj) > 0. ) THEN ; Rr = Rjm
1594	ELSE ; Rr = Rjp
1595	ENDIF
1596
1597	zh3 = pfu_ho(ji,jj) - pfu_ups(ji,jj)
1598	IF( Rj > 0. ) THEN
1599	zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(puc(ji,jj)), &
1600	& MIN( 2. * Rr * 0.5 * ABS(puc(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) )
1601	ELSE
1602	zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(puc(ji,jj)), &
1603	& MIN(-2. * Rr * 0.5 * ABS(puc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) )
1604	ENDIF
1605	pfu_ho(ji,jj) = pfu_ups(ji,jj) + zlimiter
1606
1607	ELSE
1608	IF( Rj /= 0. ) THEN
1609	IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj
1610	ELSE ; Cr = Rjp / Rj
1611	ENDIF
1612	ELSE
1613	Cr = 0.
1614	!IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm * 1.e20
1615	!ELSE ; Cr = Rjp * 1.e20
1616	!ENDIF
1617	ENDIF
1618
1619	! -- superbee --
1620	zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
1621	! -- van albada 2 --
1622	!!zpsi = 2.Cr / (CrCr+1.)
1623
1624	! -- sweby (with beta=1) --
1625	!!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
1626	! -- van Leer --
1627	!!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
1628	! -- ospre --
1629	!!zpsi = 1.5 * ( CrCr + Cr ) / ( CrCr + Cr + 1. )
1630	! -- koren --
1631	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( (1.+2Cr)/3., 2. ) ) )
1632	! -- charm --
1633	!IF( Cr > 0. ) THEN ; zpsi = Cr * (3.Cr + 1.) / ( (Cr + 1.) (Cr + 1.) )
1634	!ELSE ; zpsi = 0.
1635	!ENDIF
1636	! -- van albada 1 --
1637	!!zpsi = (CrCr + Cr) / (CrCr +1)
1638	! -- smart --
1639	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, 4. ) ) )
1640	! -- umist --
1641	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
1642
1643	! high order flux corrected by the limiter
1644	pfu_ho(ji,jj) = pfu_ho(ji,jj) - ABS( puc(ji,jj) ) * ( (1.-zpsi) + uCFLzpsi ) Rj * 0.5
1645
1646	ENDIF
1647	END DO
1648	END DO
1649	CALL lbc_lnk( pfu_ho, 'U', -1.) ! lateral boundary cond.
1650	!
1651	END SUBROUTINE limiter_x
1652
1653	SUBROUTINE limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
1654	!!---------------------------------------------------------------------
1655	!! * ROUTINE limiter_y *
1656	!!
1657	!! ** Purpose : compute flux limiter
1658	!!----------------------------------------------------------------------
1659	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
1660	REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pv ! ice i-velocity => u*e2
1661	REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pvc ! ice i-velocity A => ue2*a
1662	REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pt ! ice tracer
1663	REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) :: pfv_ho ! high order flux
1664	REAL(wp), DIMENSION (jpi,jpj), INTENT(in ), OPTIONAL :: pfv_ups ! upstream flux
1665	!
1666	REAL(wp) :: Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr
1667	INTEGER :: ji, jj ! dummy loop indices
1668	REAL(wp), DIMENSION (jpi,jpj) :: zslpy ! tracer slopes
1669	!!----------------------------------------------------------------------
1670	!
1671	DO jj = 2, jpjm1
1672	DO ji = fs_2, fs_jpim1 ! vector opt.
1673	zslpy(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * vmask(ji,jj,1)
1674	END DO
1675	END DO
1676	CALL lbc_lnk( zslpy, 'V', -1.) ! lateral boundary cond.
1677
1678	DO jj = 2, jpjm1
1679	DO ji = fs_2, fs_jpim1 ! vector opt.
1680	vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj)
1681
1682	Rjm = zslpy(ji,jj-1)
1683	Rj = zslpy(ji,jj )
1684	Rjp = zslpy(ji,jj+1)
1685
1686	IF( PRESENT(pfv_ups) ) THEN
1687
1688	IF( pv(ji,jj) > 0. ) THEN ; Rr = Rjm
1689	ELSE ; Rr = Rjp
1690	ENDIF
1691
1692	zh3 = pfv_ho(ji,jj) - pfv_ups(ji,jj)
1693	IF( Rj > 0. ) THEN
1694	zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pvc(ji,jj)), &
1695	& MIN( 2. * Rr * 0.5 * ABS(pvc(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) )
1696	ELSE
1697	zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pvc(ji,jj)), &
1698	& MIN(-2. * Rr * 0.5 * ABS(pvc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) )
1699	ENDIF
1700	pfv_ho(ji,jj) = pfv_ups(ji,jj) + zlimiter
1701
1702	ELSE
1703
1704	IF( Rj /= 0. ) THEN
1705	IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj
1706	ELSE ; Cr = Rjp / Rj
1707	ENDIF
1708	ELSE
1709	Cr = 0.
1710	!IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm * 1.e20
1711	!ELSE ; Cr = Rjp * 1.e20
1712	!ENDIF
1713	ENDIF
1714
1715	! -- superbee --
1716	zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
1717	! -- van albada 2 --
1718	!!zpsi = 2.Cr / (CrCr+1.)
1719
1720	! -- sweby (with beta=1) --
1721	!!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
1722	! -- van Leer --
1723	!!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
1724	! -- ospre --
1725	!!zpsi = 1.5 * ( CrCr + Cr ) / ( CrCr + Cr + 1. )
1726	! -- koren --
1727	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( (1.+2Cr)/3., 2. ) ) )
1728	! -- charm --
1729	!IF( Cr > 0. ) THEN ; zpsi = Cr * (3.Cr + 1.) / ( (Cr + 1.) (Cr + 1.) )
1730	!ELSE ; zpsi = 0.
1731	!ENDIF
1732	! -- van albada 1 --
1733	!!zpsi = (CrCr + Cr) / (CrCr +1)
1734	! -- smart --
1735	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, 4. ) ) )
1736	! -- umist --
1737	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
1738
1739	! high order flux corrected by the limiter
1740	pfv_ho(ji,jj) = pfv_ho(ji,jj) - ABS( pvc(ji,jj) ) * ( (1.-zpsi) + vCFLzpsi ) Rj * 0.5
1741
1742	ENDIF
1743	END DO
1744	END DO
1745	CALL lbc_lnk( pfv_ho, 'V', -1.) ! lateral boundary cond.
1746	!
1747	END SUBROUTINE limiter_y
1748
1749	#else
1750	!!----------------------------------------------------------------------
1751	!! Default option Dummy module NO SI3 sea-ice model
1752	!!----------------------------------------------------------------------
1753	#endif
1754
1755	!!======================================================================
1756	END MODULE icedyn_adv_umx

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

Download in other formats: