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

Last change on this file since 10519 was 10519, checked in by clem, 5 years ago
solve a compilation issue with agrif. It looks like agrif does not like 2 subroutines to have the same name eventhough they are not in the same module and they are not declared as public. I would like to say wtf but lets remain polite and say whats going on? Does agrif consider all the subroutines as public?
Property svn:keywords set to `Id`
File size: 67.9 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_ice : 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	! limiter: 1=nonosc_ice, 2=superbee, 3=h3(rachid)
38	INTEGER :: kn_limiter = 1
39
40	! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme
41	! clem: if set to true, the 2D test case "diagonal advection" does not work (I do not understand why)
42	LOGICAL :: ll_neg = .FALSE.
43
44	! alternate directions for upstream
45	LOGICAL :: ll_upsxy = .TRUE.
46
47	! alternate directions for high order
48	LOGICAL :: ll_hoxy = .TRUE.
49
50	! prelimiter: use it to avoid overshoot in H
51	! clem: if set to true, the 2D test case "diagnoal advection" does not work (I do not understand why)
52	LOGICAL :: ll_prelimiter_zalesak = .FALSE. ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D
53
54
55	!! * Substitutions
56	# include "vectopt_loop_substitute.h90"
57	!!----------------------------------------------------------------------
58	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
59	!! $Id$
60	!! Software governed by the CeCILL licence (./LICENSE)
61	!!----------------------------------------------------------------------
62	CONTAINS
63
64	SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, &
65	& pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
66	!!----------------------------------------------------------------------
67	!! * ROUTINE ice_dyn_adv_umx *
68	!!
69	!! ** Purpose : Compute the now trend due to total advection of
70	!! tracers and add it to the general trend of tracer equations
71	!! using an "Ultimate-Macho" scheme
72	!!
73	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
74	!!----------------------------------------------------------------------
75	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
76	INTEGER , INTENT(in ) :: kt ! time step
77	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity
78	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pv_ice ! ice j-velocity
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 fraction
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
94	REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! send zcflnow and receive zcflprv
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) :: zua_ho, zva_ho
98	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_ai, z1_aip
99	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhvar
100	!!----------------------------------------------------------------------
101	!
102	IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
103	!
104	! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !
105	! When needed, the advection split is applied at the next time-step in order to avoid blocking global comm.
106	! ...this should not affect too much the stability... Was ok on the tests we did...
107	zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
108	zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
109
110	! non-blocking global communication send zcflnow and receive zcflprv
111	CALL mpp_delay_max( 'icedyn_adv_umx', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 )
112
113	IF( zcflprv(1) > .5 ) THEN ; icycle = 2
114	ELSE ; icycle = 1
115	ENDIF
116
117	zdt = rdt_ice / REAL(icycle)
118
119	! --- transport --- !
120	zudy(:,:) = pu_ice(:,:) * e2u(:,:)
121	zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
122
123	! --- define velocity for advection: u*grad(H) --- !
124	DO jj = 2, jpjm1
125	DO ji = fs_2, fs_jpim1
126	IF ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN ; zcu_box(ji,jj) = 0._wp
127	ELSEIF( pu_ice(ji,jj) > 0._wp ) THEN ; zcu_box(ji,jj) = pu_ice(ji-1,jj)
128	ELSE ; zcu_box(ji,jj) = pu_ice(ji ,jj)
129	ENDIF
130
131	IF ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN ; zcv_box(ji,jj) = 0._wp
132	ELSEIF( pv_ice(ji,jj) > 0._wp ) THEN ; zcv_box(ji,jj) = pv_ice(ji,jj-1)
133	ELSE ; zcv_box(ji,jj) = pv_ice(ji,jj )
134	ENDIF
135	END DO
136	END DO
137
138	!---------------!
139	!== advection ==!
140	!---------------!
141	DO jt = 1, icycle
142
143	! record at_i before advection (for open water)
144	zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
145
146	! inverse of A and Ap
147	WHERE( pa_i(:,:,:) >= epsi20 ) ; z1_ai(:,:,:) = 1._wp / pa_i(:,:,:)
148	ELSEWHERE ; z1_ai(:,:,:) = 0.
149	END WHERE
150	WHERE( pa_ip(:,:,:) >= epsi20 ) ; z1_aip(:,:,:) = 1._wp / pa_ip(:,:,:)
151	ELSEWHERE ; z1_aip(:,:,:) = 0.
152	END WHERE
153	!
154	! set u*a=u for advection of A only
155	DO jl = 1, jpl
156	zua_ho(:,:,jl) = zudy(:,:)
157	zva_ho(:,:,jl) = zvdx(:,:)
158	END DO
159
160	zamsk = 1._wp
161	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_i, pa_i, zua_ho, zva_ho ) ! Ice area
162	zamsk = 0._wp
163	!
164	zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
165	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_i ) ! Ice volume
166	!
167	zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
168	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_s ) ! Snw volume
169	!
170	zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:)
171	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, psv_i ) ! Salt content
172	!
173	zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:)
174	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, poa_i ) ! Age content
175	!
176	DO jk = 1, nlay_i
177	zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:)
178	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_i(:,:,jk,:) ) ! Ice heat content
179	END DO
180	!
181	DO jk = 1, nlay_s
182	zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:)
183	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_s(:,:,jk,:) ) ! Snw heat content
184	END DO
185	!
186	IF ( ln_pnd_H12 ) THEN
187	! set u*a=u for advection of Ap only
188	DO jl = 1, jpl
189	zua_ho(:,:,jl) = zudy(:,:)
190	zva_ho(:,:,jl) = zvdx(:,:)
191	END DO
192
193	zamsk = 1._wp
194	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_ip, pa_ip, zua_ho, zva_ho ) ! mp fraction
195	zamsk = 0._wp
196	!
197	zhvar(:,:,:) = pv_ip(:,:,:) * z1_aip(:,:,:)
198	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_ip ) ! mp volume
199	ENDIF
200	!
201	zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
202	DO jj = 2, jpjm1
203	DO ji = fs_2, fs_jpim1
204	pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & ! Open water area
205	& - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
206	END DO
207	END DO
208	CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T', 1. )
209	!
210	END DO
211	!
212	END SUBROUTINE ice_dyn_adv_umx
213
214
215	SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho )
216	!!----------------------------------------------------------------------
217	!! * ROUTINE adv_umx *
218	!!
219	!! ** Purpose : Compute the now trend due to total advection of
220	!! tracers and add it to the general trend of tracer equations
221	!!
222	!! ** Method : - calculate upstream fluxes and upstream solution for tracer H
223	!! - calculate tracer H at u and v points (Ultimate)
224	!! - calculate the high order fluxes using alterning directions (Macho?)
225	!! - apply a limiter on the fluxes (nonosc_ice)
226	!! - convert this tracer flux to a tracer content flux (uH -> uV)
227	!! - calculate the high order solution for tracer content V
228	!!
229	!! ** Action : solve 2 equations => a) da/dt = -div(ua)
230	!! b) dV/dt = -div(uV) using dH/dt = -u.grad(H)
231	!! in eq. b), - fluxes uH are evaluated (with UMx) and limited (with nonosc_ice). This step is necessary to get a good H.
232	!! - then we convert this flux to a "volume" flux this way => uH*ua/u
233	!! where ua is the flux from eq. a)
234	!! - at last we estimate dV/dt = -div(uH*ua/u)
235	!!
236	!! ** Note : - this method can lead to small negative V (since we only limit H) => corrected in icedyn_adv.F90 conserving mass etc.
237	!! - negative tracers at u-v points can also occur from the Ultimate scheme (usually at the ice edge) and the solution for now
238	!! is to apply an upstream scheme when it occurs. A better solution would be to degrade the order of
239	!! the scheme automatically by applying a mask of the ice cover inside Ultimate (not done).
240	!! - Eventhough 1D tests give very good results (typically the one from Schar & Smolarkiewiecz), the 2D is less good.
241	!! Large values of H can appear for very small ice concentration, and when it does it messes the things up since we
242	!! work on H (and not V). It probably comes from the prelimiter of zalesak which is coded for 1D and not 2D.
243	!! Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice
244	!! concentration is small).
245	!! To-do: expand the prelimiter from zalesak to make it work in 2D
246	!!----------------------------------------------------------------------
247	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
248	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
249	INTEGER , INTENT(in ) :: jt ! number of sub-iteration
250	INTEGER , INTENT(in ) :: kt ! number of iteration
251	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
252	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu , pv ! 2 ice velocity components => u*e2
253	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: puc , pvc ! 2 ice velocity components => ue2 or ua*e2u
254	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pubox, pvbox ! upstream velocity
255	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pt ! tracer field
256	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: ptc ! tracer content field
257	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out), OPTIONAL :: pua_ho, pva_ho ! high order u*a fluxes
258	!
259	INTEGER :: ji, jj, jl ! dummy loop indices
260	REAL(wp) :: ztra ! local scalar
261	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zfu_ho , zfv_ho , zpt
262	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zfu_ups, zfv_ups, zt_ups
263	!!----------------------------------------------------------------------
264	!
265	! Upstream (_ups) fluxes
266	! -----------------------
267	CALL upstream( pamsk, jt, kt, pdt, pt, pu, pv, zt_ups, zfu_ups, zfv_ups )
268
269	! High order (_ho) fluxes
270	! -----------------------
271	SELECT CASE( kn_umx )
272	!
273	CASE ( 20 ) !== centered second order ==!
274	!
275	CALL cen2( pamsk, jt, kt, pdt, pt, pu, pv, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho )
276	!
277	CASE ( 1:5 ) !== 1st to 5th order ULTIMATE-MACHO scheme ==!
278	!
279	CALL macho( pamsk, kn_umx, jt, kt, pdt, pt, pu, pv, pubox, pvbox, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho )
280	!
281	END SELECT
282	!
283	! --ho --ho
284	! new fluxes = uH u*a / u
285	! ----------------------------
286	IF( pamsk == 0._wp ) THEN
287	DO jl = 1, jpl
288	DO jj = 1, jpjm1
289	DO ji = 1, fs_jpim1
290	IF( ABS( puc(ji,jj,jl) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN
291	zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj)
292	zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj)
293	ELSE
294	zfu_ho (ji,jj,jl) = 0._wp
295	zfu_ups(ji,jj,jl) = 0._wp
296	ENDIF
297	!
298	IF( ABS( pvc(ji,jj,jl) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN
299	zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj)
300	zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj)
301	ELSE
302	zfv_ho (ji,jj,jl) = 0._wp
303	zfv_ups(ji,jj,jl) = 0._wp
304	ENDIF
305	END DO
306	END DO
307	END DO
308	ENDIF
309	! --ho
310	! in case of advection of A: output u*a
311	! -------------------------------------
312	IF( PRESENT( pua_ho ) ) THEN
313	DO jl = 1, jpl
314	DO jj = 1, jpjm1
315	DO ji = 1, fs_jpim1
316	pua_ho(ji,jj,jl) = zfu_ho(ji,jj,jl)
317	pva_ho(ji,jj,jl) = zfv_ho(ji,jj,jl)
318	END DO
319	END DO
320	END DO
321	ENDIF
322	!
323	! final trend with corrected fluxes
324	! ---------------------------------
325	DO jl = 1, jpl
326	DO jj = 2, jpjm1
327	DO ji = fs_2, fs_jpim1
328	ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) )
329	!
330	ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)
331	END DO
332	END DO
333	END DO
334	CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T', 1. )
335	!
336	END SUBROUTINE adv_umx
337
338
339	SUBROUTINE upstream( pamsk, jt, kt, pdt, pt, pu, pv, pt_ups, pfu_ups, pfv_ups )
340	!!---------------------------------------------------------------------
341	!! * ROUTINE upstream *
342	!!
343	!! ** Purpose : compute the upstream fluxes and upstream guess of tracer
344	!!----------------------------------------------------------------------
345	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
346	INTEGER , INTENT(in ) :: jt ! number of sub-iteration
347	INTEGER , INTENT(in ) :: kt ! number of iteration
348	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
349	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields
350	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu, pv ! 2 ice velocity components
351	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_ups ! upstream guess of tracer
352	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ups, pfv_ups ! upstream fluxes
353	!
354	INTEGER :: ji, jj, jl ! dummy loop indices
355	REAL(wp) :: ztra ! local scalar
356	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zpt
357	!!----------------------------------------------------------------------
358
359	IF( .NOT. ll_upsxy ) THEN ! no alternate directions !
360	!
361	DO jl = 1, jpl
362	DO jj = 1, jpjm1
363	DO ji = 1, fs_jpim1
364	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)
365	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)
366	END DO
367	END DO
368	END DO
369	!
370	ELSE ! alternate directions !
371	!
372	IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==!
373	!
374	DO jl = 1, jpl !-- flux in x-direction
375	DO jj = 1, jpjm1
376	DO ji = 1, fs_jpim1
377	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)
378	END DO
379	END DO
380	END DO
381	!
382	DO jl = 1, jpl !-- first guess of tracer from u-flux
383	DO jj = 2, jpjm1
384	DO ji = fs_2, fs_jpim1
385	ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) ) &
386	& + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk)
387	!
388	zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
389	END DO
390	END DO
391	END DO
392	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
393	!
394	DO jl = 1, jpl !-- flux in y-direction
395	DO jj = 1, jpjm1
396	DO ji = 1, fs_jpim1
397	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)
398	END DO
399	END DO
400	END DO
401	!
402	ELSE !== even ice time step: adv_y then adv_x ==!
403	!
404	DO jl = 1, jpl !-- flux in y-direction
405	DO jj = 1, jpjm1
406	DO ji = 1, fs_jpim1
407	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)
408	END DO
409	END DO
410	END DO
411	!
412	DO jl = 1, jpl !-- first guess of tracer from v-flux
413	DO jj = 2, jpjm1
414	DO ji = fs_2, fs_jpim1
415	ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) ) &
416	& + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk)
417	!
418	zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
419	END DO
420	END DO
421	END DO
422	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
423	!
424	DO jl = 1, jpl !-- flux in x-direction
425	DO jj = 1, jpjm1
426	DO ji = 1, fs_jpim1
427	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)
428	END DO
429	END DO
430	END DO
431	!
432	ENDIF
433
434	ENDIF
435	!
436	DO jl = 1, jpl !-- after tracer with upstream scheme
437	DO jj = 2, jpjm1
438	DO ji = fs_2, fs_jpim1
439	ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj ,jl) &
440	& + pfv_ups(ji,jj,jl) - pfv_ups(ji ,jj-1,jl) ) &
441	& + ( pu (ji,jj ) - pu (ji-1,jj ) &
442	& + pv (ji,jj ) - pv (ji ,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk)
443	!
444	pt_ups(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
445	END DO
446	END DO
447	END DO
448	CALL lbc_lnk( 'icedyn_adv_umx', pt_ups, 'T', 1. )
449
450	END SUBROUTINE upstream
451
452
453	SUBROUTINE cen2( pamsk, jt, kt, pdt, pt, pu, pv, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
454	!!---------------------------------------------------------------------
455	!! * ROUTINE cen2 *
456	!!
457	!! ** Purpose : compute the high order fluxes using a centered
458	!! second order scheme
459	!!----------------------------------------------------------------------
460	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
461	INTEGER , INTENT(in ) :: jt ! number of sub-iteration
462	INTEGER , INTENT(in ) :: kt ! number of iteration
463	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
464	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields
465	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu, pv ! 2 ice velocity components
466	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt_ups ! upstream guess of tracer
467	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pfu_ups, pfv_ups ! upstream fluxes
468	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes
469	!
470	INTEGER :: ji, jj, jl ! dummy loop indices
471	REAL(wp) :: ztra ! local scalar
472	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zpt
473	!!----------------------------------------------------------------------
474	!
475	IF( .NOT.ll_hoxy ) THEN ! no alternate directions !
476	!
477	DO jl = 1, jpl
478	DO jj = 1, jpjm1
479	DO ji = 1, fs_jpim1
480	pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj ,jl) )
481	pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji ,jj+1,jl) )
482	END DO
483	END DO
484	END DO
485	!
486	IF ( kn_limiter == 1 ) THEN
487	CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
488	ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN
489	CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
490	CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
491	ENDIF
492	!
493	ELSE ! alternate directions !
494	!
495	IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==!
496	!
497	DO jl = 1, jpl !-- flux in x-direction
498	DO jj = 1, jpjm1
499	DO ji = 1, fs_jpim1
500	pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
501	END DO
502	END DO
503	END DO
504	IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
505
506	DO jl = 1, jpl !-- first guess of tracer from u-flux
507	DO jj = 2, jpjm1
508	DO ji = fs_2, fs_jpim1
509	ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) &
510	& + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk)
511	!
512	zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
513	END DO
514	END DO
515	END DO
516	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
517
518	DO jl = 1, jpl !-- flux in y-direction
519	DO jj = 1, jpjm1
520	DO ji = 1, fs_jpim1
521	pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) )
522	END DO
523	END DO
524	END DO
525	IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
526
527	ELSE !== even ice time step: adv_y then adv_x ==!
528	!
529	DO jl = 1, jpl !-- flux in y-direction
530	DO jj = 1, jpjm1
531	DO ji = 1, fs_jpim1
532	pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
533	END DO
534	END DO
535	END DO
536	IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
537	!
538	DO jl = 1, jpl !-- first guess of tracer from v-flux
539	DO jj = 2, jpjm1
540	DO ji = fs_2, fs_jpim1
541	ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) &
542	& + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk)
543	!
544	zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
545	END DO
546	END DO
547	END DO
548	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
549	!
550	DO jl = 1, jpl !-- flux in x-direction
551	DO jj = 1, jpjm1
552	DO ji = 1, fs_jpim1
553	pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) )
554	END DO
555	END DO
556	END DO
557	IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
558
559	ENDIF
560	IF( kn_limiter == 1 ) CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
561
562	ENDIF
563
564	END SUBROUTINE cen2
565
566
567	SUBROUTINE macho( pamsk, kn_umx, jt, kt, pdt, pt, pu, pv, pubox, pvbox, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
568	!!---------------------------------------------------------------------
569	!! * ROUTINE macho *
570	!!
571	!! ** Purpose : compute the high order fluxes using Ultimate-Macho scheme
572	!!
573	!! ** Method : ...
574	!!
575	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
576	!!----------------------------------------------------------------------
577	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
578	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
579	INTEGER , INTENT(in ) :: jt ! number of sub-iteration
580	INTEGER , INTENT(in ) :: kt ! number of iteration
581	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
582	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields
583	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu, pv ! 2 ice velocity components
584	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pubox, pvbox ! upstream velocity
585	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt_ups ! upstream guess of tracer
586	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pfu_ups, pfv_ups ! upstream fluxes
587	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes
588	!
589	INTEGER :: ji, jj, jl ! dummy loop indices
590	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zt_u, zt_v, zpt
591	!!----------------------------------------------------------------------
592	!
593	IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==!
594	!
595	! !-- ultimate interpolation of pt at u-point --!
596	CALL ultimate_x( kn_umx, pdt, pt, pu, zt_u, pfu_ho )
597	! !-- limiter in x --!
598	IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
599	! !-- advective form update in zpt --!
600	DO jl = 1, jpl
601	DO jj = 2, jpjm1
602	DO ji = fs_2, fs_jpim1
603	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) &
604	& + pt (ji,jj,jl) * ( pu (ji,jj ) - pu (ji-1,jj ) ) * r1_e1e2t(ji,jj) &
605	& * pamsk &
606	& ) * pdt ) * tmask(ji,jj,1)
607	END DO
608	END DO
609	END DO
610	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
611	!
612	! !-- ultimate interpolation of pt at v-point --!
613	IF( ll_hoxy ) THEN
614	CALL ultimate_y( kn_umx, pdt, zpt, pv, zt_v, pfv_ho )
615	ELSE
616	CALL ultimate_y( kn_umx, pdt, pt , pv, zt_v, pfv_ho )
617	ENDIF
618	! !-- limiter in y --!
619	IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
620	!
621	!
622	ELSE !== even ice time step: adv_y then adv_x ==!
623	!
624	! !-- ultimate interpolation of pt at v-point --!
625	CALL ultimate_y( kn_umx, pdt, pt, pv, zt_v, pfv_ho )
626	! !-- limiter in y --!
627	IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
628	! !-- advective form update in zpt --!
629	DO jl = 1, jpl
630	DO jj = 2, jpjm1
631	DO ji = fs_2, fs_jpim1
632	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) &
633	& + pt (ji,jj,jl) * ( pv (ji,jj ) - pv (ji,jj-1 ) ) * r1_e1e2t(ji,jj) &
634	& * pamsk &
635	& ) * pdt ) * tmask(ji,jj,1)
636	END DO
637	END DO
638	END DO
639	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
640	!
641	! !-- ultimate interpolation of pt at u-point --!
642	IF( ll_hoxy ) THEN
643	CALL ultimate_x( kn_umx, pdt, zpt, pu, zt_u, pfu_ho )
644	ELSE
645	CALL ultimate_x( kn_umx, pdt, pt , pu, zt_u, pfu_ho )
646	ENDIF
647	! !-- limiter in x --!
648	IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
649	!
650	ENDIF
651
652	IF( kn_limiter == 1 ) CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
653	!
654	END SUBROUTINE macho
655
656
657	SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho )
658	!!---------------------------------------------------------------------
659	!! * ROUTINE ultimate_x *
660	!!
661	!! ** Purpose : compute tracer at u-points
662	!!
663	!! ** Method : ...
664	!!
665	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
666	!!----------------------------------------------------------------------
667	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
668	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
669	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu ! ice i-velocity component
670	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields
671	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_u ! tracer at u-point
672	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho ! high order flux
673	!
674	INTEGER :: ji, jj, jl ! dummy loop indices
675	REAL(wp) :: zcu, zdx2, zdx4 ! - -
676	REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztu1, ztu2, ztu3, ztu4
677	!!----------------------------------------------------------------------
678	!
679	! !-- Laplacian in i-direction --!
680	DO jl = 1, jpl
681	DO jj = 2, jpjm1 ! First derivative (gradient)
682	DO ji = 1, fs_jpim1
683	ztu1(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
684	END DO
685	! ! Second derivative (Laplacian)
686	DO ji = fs_2, fs_jpim1
687	ztu2(ji,jj,jl) = ( ztu1(ji,jj,jl) - ztu1(ji-1,jj,jl) ) * r1_e1t(ji,jj)
688	END DO
689	END DO
690	END DO
691	CALL lbc_lnk( 'icedyn_adv_umx', ztu2, 'T', 1. )
692	!
693	! !-- BiLaplacian in i-direction --!
694	DO jl = 1, jpl
695	DO jj = 2, jpjm1 ! Third derivative
696	DO ji = 1, fs_jpim1
697	ztu3(ji,jj,jl) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
698	END DO
699	! ! Fourth derivative
700	DO ji = fs_2, fs_jpim1
701	ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj)
702	END DO
703	END DO
704	END DO
705	CALL lbc_lnk( 'icedyn_adv_umx', ztu4, 'T', 1. )
706	!
707	!
708	SELECT CASE (kn_umx )
709	!
710	CASE( 1 ) !== 1st order central TIM ==! (Eq. 21)
711	!
712	DO jl = 1, jpl
713	DO jj = 1, jpjm1
714	DO ji = 1, fs_jpim1 ! vector opt.
715	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) &
716	& - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
717	END DO
718	END DO
719	END DO
720	!
721	CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23)
722	!
723	DO jl = 1, jpl
724	DO jj = 1, jpjm1
725	DO ji = 1, fs_jpim1 ! vector opt.
726	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
727	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) &
728	& - zcu * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
729	END DO
730	END DO
731	END DO
732	!
733	CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24)
734	!
735	DO jl = 1, jpl
736	DO jj = 1, jpjm1
737	DO ji = 1, fs_jpim1 ! vector opt.
738	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
739	zdx2 = e1u(ji,jj) * e1u(ji,jj)
740	!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
741	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
742	& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
743	& + z1_6 * zdx2 * ( zcuzcu - 1._wp ) ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
744	& - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
745	END DO
746	END DO
747	END DO
748	!
749	CASE( 4 ) !== 4th order central TIM ==! (Eq. 27)
750	!
751	DO jl = 1, jpl
752	DO jj = 1, jpjm1
753	DO ji = 1, fs_jpim1 ! vector opt.
754	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
755	zdx2 = e1u(ji,jj) * e1u(ji,jj)
756	!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
757	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
758	& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
759	& + z1_6 * zdx2 * ( zcuzcu - 1._wp ) ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
760	& - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
761	END DO
762	END DO
763	END DO
764	!
765	CASE( 5 ) !== 5th order central TIM ==! (Eq. 29)
766	!
767	DO jl = 1, jpl
768	DO jj = 1, jpjm1
769	DO ji = 1, fs_jpim1 ! vector opt.
770	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
771	zdx2 = e1u(ji,jj) * e1u(ji,jj)
772	!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
773	zdx4 = zdx2 * zdx2
774	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
775	& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
776	& + z1_6 * zdx2 * ( zcuzcu - 1._wp ) ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
777	& - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) &
778	& + z1_120 * zdx4 * ( zcuzcu - 1._wp ) ( zcuzcu - 4._wp ) ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl) &
779	& - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) )
780	END DO
781	END DO
782	END DO
783	!
784	END SELECT
785	!
786	! if pt at u-point is negative then use the upstream value
787	! this should not be necessary if a proper sea-ice mask is set in Ultimate
788	! to degrade the order of the scheme when necessary (for ex. at the ice edge)
789	IF( ll_neg ) THEN
790	DO jl = 1, jpl
791	DO jj = 1, jpjm1
792	DO ji = 1, fs_jpim1
793	IF( pt_u(ji,jj,jl) < 0._wp ) THEN
794	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) &
795	& - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
796	ENDIF
797	END DO
798	END DO
799	END DO
800	ENDIF
801	! !-- High order flux in i-direction --!
802	DO jl = 1, jpl
803	DO jj = 1, jpjm1
804	DO ji = 1, fs_jpim1 ! vector opt.
805	pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl)
806	END DO
807	END DO
808	END DO
809	!
810	END SUBROUTINE ultimate_x
811
812
813	SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho )
814	!!---------------------------------------------------------------------
815	!! * ROUTINE ultimate_y *
816	!!
817	!! ** Purpose : compute tracer at v-points
818	!!
819	!! ** Method : ...
820	!!
821	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
822	!!----------------------------------------------------------------------
823	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
824	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
825	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pv ! ice j-velocity component
826	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields
827	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_v ! tracer at v-point
828	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfv_ho ! high order flux
829	!
830	INTEGER :: ji, jj, jl ! dummy loop indices
831	REAL(wp) :: zcv, zdy2, zdy4 ! - -
832	REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztv1, ztv2, ztv3, ztv4
833	!!----------------------------------------------------------------------
834	!
835	! !-- Laplacian in j-direction --!
836	DO jl = 1, jpl
837	DO jj = 1, jpjm1 ! First derivative (gradient)
838	DO ji = fs_2, fs_jpim1
839	ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
840	END DO
841	END DO
842	DO jj = 2, jpjm1 ! Second derivative (Laplacian)
843	DO ji = fs_2, fs_jpim1
844	ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj)
845	END DO
846	END DO
847	END DO
848	CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1. )
849	!
850	! !-- BiLaplacian in j-direction --!
851	DO jl = 1, jpl
852	DO jj = 1, jpjm1 ! First derivative
853	DO ji = fs_2, fs_jpim1
854	ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
855	END DO
856	END DO
857	DO jj = 2, jpjm1 ! Second derivative
858	DO ji = fs_2, fs_jpim1
859	ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj)
860	END DO
861	END DO
862	END DO
863	CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1. )
864	!
865	!
866	SELECT CASE (kn_umx )
867	!
868	CASE( 1 ) !== 1st order central TIM ==! (Eq. 21)
869	DO jl = 1, jpl
870	DO jj = 1, jpjm1
871	DO ji = 1, fs_jpim1
872	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) &
873	& - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
874	END DO
875	END DO
876	END DO
877	!
878	CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23)
879	DO jl = 1, jpl
880	DO jj = 1, jpjm1
881	DO ji = 1, fs_jpim1
882	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
883	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) &
884	& - zcv * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
885	END DO
886	END DO
887	END DO
888	!
889	CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24)
890	DO jl = 1, jpl
891	DO jj = 1, jpjm1
892	DO ji = 1, fs_jpim1
893	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
894	zdy2 = e2v(ji,jj) * e2v(ji,jj)
895	!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
896	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
897	& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
898	& + z1_6 * zdy2 * ( zcvzcv - 1._wp ) ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
899	& - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
900	END DO
901	END DO
902	END DO
903	!
904	CASE( 4 ) !== 4th order central TIM ==! (Eq. 27)
905	DO jl = 1, jpl
906	DO jj = 1, jpjm1
907	DO ji = 1, fs_jpim1
908	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
909	zdy2 = e2v(ji,jj) * e2v(ji,jj)
910	!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
911	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
912	& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
913	& + z1_6 * zdy2 * ( zcvzcv - 1._wp ) ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
914	& - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
915	END DO
916	END DO
917	END DO
918	!
919	CASE( 5 ) !== 5th order central TIM ==! (Eq. 29)
920	DO jl = 1, jpl
921	DO jj = 1, jpjm1
922	DO ji = 1, fs_jpim1
923	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
924	zdy2 = e2v(ji,jj) * e2v(ji,jj)
925	!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
926	zdy4 = zdy2 * zdy2
927	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
928	& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
929	& + z1_6 * zdy2 * ( zcvzcv - 1._wp ) ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
930	& - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) &
931	& + z1_120 * zdy4 * ( zcvzcv - 1._wp ) ( zcvzcv - 4._wp ) ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl) &
932	& - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) )
933	END DO
934	END DO
935	END DO
936	!
937	END SELECT
938	!
939	! if pt at v-point is negative then use the upstream value
940	! this should not be necessary if a proper sea-ice mask is set in Ultimate
941	! to degrade the order of the scheme when necessary (for ex. at the ice edge)
942	IF( ll_neg ) THEN
943	DO jl = 1, jpl
944	DO jj = 1, jpjm1
945	DO ji = 1, fs_jpim1
946	IF( pt_v(ji,jj,jl) < 0._wp ) THEN
947	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) &
948	& - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
949	ENDIF
950	END DO
951	END DO
952	END DO
953	ENDIF
954	! !-- High order flux in j-direction --!
955	DO jl = 1, jpl
956	DO jj = 1, jpjm1
957	DO ji = 1, fs_jpim1 ! vector opt.
958	pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl)
959	END DO
960	END DO
961	END DO
962	!
963	END SUBROUTINE ultimate_y
964
965
966	SUBROUTINE nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
967	!!---------------------------------------------------------------------
968	!! * ROUTINE nonosc_ice *
969	!!
970	!! ** Purpose : compute monotonic tracer fluxes from the upstream
971	!! scheme and the before field by a non-oscillatory algorithm
972	!!
973	!! ** Method : ...
974	!!----------------------------------------------------------------------
975	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
976	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
977	REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pu ! ice i-velocity => u*e2
978	REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pv ! ice j-velocity => v*e1
979	REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pt, pt_ups ! before field & upstream guess of after field
980	REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pfv_ups, pfu_ups ! upstream flux
981	REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfv_ho, pfu_ho ! monotonic flux
982	!
983	INTEGER :: ji, jj, jl ! dummy loop indices
984	REAL(wp) :: zpos, zneg, zbig, zup, zdo, z1_dt ! local scalars
985	REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zcoef, zzt ! - -
986	REAL(wp), DIMENSION(jpi,jpj ) :: zbup, zbdo
987	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zbetup, zbetdo, zti_ups, ztj_ups
988	!!----------------------------------------------------------------------
989	zbig = 1.e+40_wp
990
991	! antidiffusive flux : high order minus low order
992	! --------------------------------------------------
993	DO jl = 1, jpl
994	DO jj = 1, jpjm1
995	DO ji = 1, fs_jpim1 ! vector opt.
996	pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)
997	pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)
998	END DO
999	END DO
1000	END DO
1001
1002	! extreme case where pfu_ho has to be zero
1003	! ----------------------------------------
1004	! pfu_ho
1005	! * --->
1006	! \| \| * \| \|
1007	! \| \| \| * \|
1008	! \| \| \| \| *
1009	! t_ups : i-1 i i+1 i+2
1010	IF( ll_prelimiter_zalesak ) THEN
1011
1012	DO jl = 1, jpl
1013	DO jj = 2, jpjm1
1014	DO ji = fs_2, fs_jpim1
1015	zti_ups(ji,jj,jl)= pt_ups(ji+1,jj ,jl)
1016	ztj_ups(ji,jj,jl)= pt_ups(ji ,jj+1,jl)
1017	END DO
1018	END DO
1019	END DO
1020	CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_ups, 'T', 1., ztj_ups, 'T', 1. )
1021
1022	DO jl = 1, jpl
1023	DO jj = 2, jpjm1
1024	DO ji = fs_2, fs_jpim1
1025	IF ( pfu_ho(ji,jj,jl) * ( pt_ups(ji+1,jj ,jl) - pt_ups(ji,jj,jl) ) <= 0._wp .AND. &
1026	& pfv_ho(ji,jj,jl) * ( pt_ups(ji ,jj+1,jl) - pt_ups(ji,jj,jl) ) <= 0._wp ) THEN
1027	!
1028	IF( pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj ,jl) - zti_ups(ji,jj,jl) ) <= 0._wp .AND. &
1029	& pfv_ho(ji,jj,jl) * ( ztj_ups(ji ,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0._wp ) THEN
1030	pfu_ho(ji,jj,jl)=0._wp
1031	pfv_ho(ji,jj,jl)=0._wp
1032	ENDIF
1033	!
1034	IF( pfu_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji-1,jj ,jl) ) <= 0._wp .AND. &
1035	& pfv_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji ,jj-1,jl) ) <= 0._wp ) THEN
1036	pfu_ho(ji,jj,jl)=0._wp
1037	pfv_ho(ji,jj,jl)=0._wp
1038	ENDIF
1039	!
1040	ENDIF
1041	END DO
1042	END DO
1043	END DO
1044	CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond.
1045
1046	ENDIF
1047
1048	! Search local extrema
1049	! --------------------
1050	! max/min of pt & pt_ups with large negative/positive value (-/+zbig) outside ice cover
1051	z1_dt = 1._wp / pdt
1052	DO jl = 1, jpl
1053
1054	DO jj = 1, jpj
1055	DO ji = 1, jpi
1056	IF ( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN
1057	zbup(ji,jj) = -zbig
1058	zbdo(ji,jj) = zbig
1059	ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) > 0._wp ) THEN
1060	zbup(ji,jj) = pt_ups(ji,jj,jl)
1061	zbdo(ji,jj) = pt_ups(ji,jj,jl)
1062	ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN
1063	zbup(ji,jj) = pt(ji,jj,jl)
1064	zbdo(ji,jj) = pt(ji,jj,jl)
1065	ELSE
1066	zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_ups(ji,jj,jl) )
1067	zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_ups(ji,jj,jl) )
1068	ENDIF
1069	END DO
1070	END DO
1071
1072	DO jj = 2, jpjm1
1073	DO ji = fs_2, fs_jpim1 ! vector opt.
1074	!
1075	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
1076	zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj), zbdo(ji+1,jj), zbdo(ji,jj-1), zbdo(ji,jj+1) )
1077	!
1078	zpos = MAX( 0._wp, pfu_ho(ji-1,jj ,jl) ) - MIN( 0._wp, pfu_ho(ji ,jj ,jl) ) & ! positive/negative part of the flux
1079	& + MAX( 0._wp, pfv_ho(ji ,jj-1,jl) ) - MIN( 0._wp, pfv_ho(ji ,jj ,jl) )
1080	zneg = MAX( 0._wp, pfu_ho(ji ,jj ,jl) ) - MIN( 0._wp, pfu_ho(ji-1,jj ,jl) ) &
1081	& + MAX( 0._wp, pfv_ho(ji ,jj ,jl) ) - MIN( 0._wp, pfv_ho(ji ,jj-1,jl) )
1082	!
1083	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) ) &
1084	& ) * ( 1. - pamsk )
1085	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) ) &
1086	& ) * ( 1. - pamsk )
1087	!
1088	! ! up & down beta terms
1089	IF( zpos > 0._wp ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt
1090	ELSE ; zbetup(ji,jj,jl) = 0._wp ! zbig
1091	ENDIF
1092	!
1093	IF( zneg > 0._wp ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
1094	ELSE ; zbetdo(ji,jj,jl) = 0._wp ! zbig
1095	ENDIF
1096	!
1097	! if all the points are outside ice cover
1098	IF( zup == -zbig ) zbetup(ji,jj,jl) = 0._wp ! zbig
1099	IF( zdo == zbig ) zbetdo(ji,jj,jl) = 0._wp ! zbig
1100	!
1101	END DO
1102	END DO
1103	END DO
1104	CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)
1105
1106
1107	! monotonic flux in the y direction
1108	! ---------------------------------
1109	DO jl = 1, jpl
1110	DO jj = 1, jpjm1
1111	DO ji = 1, fs_jpim1 ! vector opt.
1112	zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) )
1113	zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) )
1114	zcu = 0.5_wp + SIGN( 0.5_wp , pfu_ho(ji,jj,jl) )
1115	!
1116	zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu )
1117	!
1118	pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_ups(ji,jj,jl)
1119	!
1120	END DO
1121	END DO
1122
1123	DO jj = 1, jpjm1
1124	DO ji = 1, fs_jpim1 ! vector opt.
1125	zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) )
1126	zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) )
1127	zcv = 0.5_wp + SIGN( 0.5_wp , pfv_ho(ji,jj,jl) )
1128	!
1129	zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv )
1130	!
1131	pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_ups(ji,jj,jl)
1132	!
1133	END DO
1134	END DO
1135
1136	! clem test
1137	!! DO jj = 2, jpjm1
1138	!! DO ji = 2, fs_jpim1 ! vector opt.
1139	!! zzt = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &
1140	!! & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &
1141	!! & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1142	!! & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1143	!! & ) * tmask(ji,jj,1)
1144	!! IF( zzt < -epsi20 ) THEN
1145	!! WRITE(numout,*) 'T<0 nonosc_ice',zzt
1146	!! ENDIF
1147	!! END DO
1148	!! END DO
1149
1150	END DO
1151	!
1152	END SUBROUTINE nonosc_ice
1153
1154
1155	SUBROUTINE limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
1156	!!---------------------------------------------------------------------
1157	!! * ROUTINE limiter_x *
1158	!!
1159	!! ** Purpose : compute flux limiter
1160	!!----------------------------------------------------------------------
1161	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
1162	REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pu ! ice i-velocity => u*e2
1163	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt ! ice tracer
1164	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pfu_ups ! upstream flux
1165	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pfu_ho ! high order flux
1166	!
1167	REAL(wp) :: Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr
1168	INTEGER :: ji, jj, jl ! dummy loop indices
1169	REAL(wp), DIMENSION (jpi,jpj,jpl) :: zslpx ! tracer slopes
1170	!!----------------------------------------------------------------------
1171	!
1172	DO jl = 1, jpl
1173	DO jj = 2, jpjm1
1174	DO ji = fs_2, fs_jpim1 ! vector opt.
1175	zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1)
1176	END DO
1177	END DO
1178	END DO
1179	CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.) ! lateral boundary cond.
1180
1181	DO jl = 1, jpl
1182	DO jj = 2, jpjm1
1183	DO ji = fs_2, fs_jpim1 ! vector opt.
1184	uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj)
1185
1186	Rjm = zslpx(ji-1,jj,jl)
1187	Rj = zslpx(ji ,jj,jl)
1188	Rjp = zslpx(ji+1,jj,jl)
1189
1190	IF( kn_limiter == 3 ) THEN
1191
1192	IF( pu(ji,jj) > 0. ) THEN ; Rr = Rjm
1193	ELSE ; Rr = Rjp
1194	ENDIF
1195
1196	zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)
1197	IF( Rj > 0. ) THEN
1198	zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pu(ji,jj)), &
1199	& MIN( 2. * Rr * 0.5 * ABS(pu(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
1200	ELSE
1201	zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)), &
1202	& MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
1203	ENDIF
1204	pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter
1205
1206	ELSEIF( kn_limiter == 2 ) THEN
1207	IF( Rj /= 0. ) THEN
1208	IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj
1209	ELSE ; Cr = Rjp / Rj
1210	ENDIF
1211	ELSE
1212	Cr = 0.
1213	ENDIF
1214
1215	! -- superbee --
1216	zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
1217	! -- van albada 2 --
1218	!!zpsi = 2.Cr / (CrCr+1.)
1219	! -- sweby (with beta=1) --
1220	!!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
1221	! -- van Leer --
1222	!!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
1223	! -- ospre --
1224	!!zpsi = 1.5 * ( CrCr + Cr ) / ( CrCr + Cr + 1. )
1225	! -- koren --
1226	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( (1.+2Cr)/3., 2. ) ) )
1227	! -- charm --
1228	!IF( Cr > 0. ) THEN ; zpsi = Cr * (3.Cr + 1.) / ( (Cr + 1.) (Cr + 1.) )
1229	!ELSE ; zpsi = 0.
1230	!ENDIF
1231	! -- van albada 1 --
1232	!!zpsi = (CrCr + Cr) / (CrCr +1)
1233	! -- smart --
1234	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, 4. ) ) )
1235	! -- umist --
1236	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
1237
1238	! high order flux corrected by the limiter
1239	pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFLzpsi ) Rj * 0.5
1240
1241	ENDIF
1242	END DO
1243	END DO
1244	END DO
1245	CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.) ! lateral boundary cond.
1246	!
1247	END SUBROUTINE limiter_x
1248
1249
1250	SUBROUTINE limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
1251	!!---------------------------------------------------------------------
1252	!! * ROUTINE limiter_y *
1253	!!
1254	!! ** Purpose : compute flux limiter
1255	!!----------------------------------------------------------------------
1256	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
1257	REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pv ! ice i-velocity => u*e2
1258	REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pt ! ice tracer
1259	REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pfv_ups ! upstream flux
1260	REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfv_ho ! high order flux
1261	!
1262	REAL(wp) :: Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr
1263	INTEGER :: ji, jj, jl ! dummy loop indices
1264	REAL(wp), DIMENSION (jpi,jpj,jpl) :: zslpy ! tracer slopes
1265	!!----------------------------------------------------------------------
1266	!
1267	DO jl = 1, jpl
1268	DO jj = 2, jpjm1
1269	DO ji = fs_2, fs_jpim1 ! vector opt.
1270	zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1)
1271	END DO
1272	END DO
1273	END DO
1274	CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.) ! lateral boundary cond.
1275
1276	DO jl = 1, jpl
1277	DO jj = 2, jpjm1
1278	DO ji = fs_2, fs_jpim1 ! vector opt.
1279	vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj)
1280
1281	Rjm = zslpy(ji,jj-1,jl)
1282	Rj = zslpy(ji,jj ,jl)
1283	Rjp = zslpy(ji,jj+1,jl)
1284
1285	IF( kn_limiter == 3 ) THEN
1286
1287	IF( pv(ji,jj) > 0. ) THEN ; Rr = Rjm
1288	ELSE ; Rr = Rjp
1289	ENDIF
1290
1291	zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)
1292	IF( Rj > 0. ) THEN
1293	zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pv(ji,jj)), &
1294	& MIN( 2. * Rr * 0.5 * ABS(pv(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
1295	ELSE
1296	zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)), &
1297	& MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
1298	ENDIF
1299	pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter
1300
1301	ELSEIF( kn_limiter == 2 ) THEN
1302
1303	IF( Rj /= 0. ) THEN
1304	IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj
1305	ELSE ; Cr = Rjp / Rj
1306	ENDIF
1307	ELSE
1308	Cr = 0.
1309	ENDIF
1310
1311	! -- superbee --
1312	zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
1313	! -- van albada 2 --
1314	!!zpsi = 2.Cr / (CrCr+1.)
1315	! -- sweby (with beta=1) --
1316	!!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
1317	! -- van Leer --
1318	!!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
1319	! -- ospre --
1320	!!zpsi = 1.5 * ( CrCr + Cr ) / ( CrCr + Cr + 1. )
1321	! -- koren --
1322	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( (1.+2Cr)/3., 2. ) ) )
1323	! -- charm --
1324	!IF( Cr > 0. ) THEN ; zpsi = Cr * (3.Cr + 1.) / ( (Cr + 1.) (Cr + 1.) )
1325	!ELSE ; zpsi = 0.
1326	!ENDIF
1327	! -- van albada 1 --
1328	!!zpsi = (CrCr + Cr) / (CrCr +1)
1329	! -- smart --
1330	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, 4. ) ) )
1331	! -- umist --
1332	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
1333
1334	! high order flux corrected by the limiter
1335	pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFLzpsi ) Rj * 0.5
1336
1337	ENDIF
1338	END DO
1339	END DO
1340	END DO
1341	CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.) ! lateral boundary cond.
1342	!
1343	END SUBROUTINE limiter_y
1344
1345	#else
1346	!!----------------------------------------------------------------------
1347	!! Default option Dummy module NO SI3 sea-ice model
1348	!!----------------------------------------------------------------------
1349	#endif
1350
1351	!!======================================================================
1352	END MODULE icedyn_adv_umx

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

Download in other formats: