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

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

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

Download in other formats: