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

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

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

Download in other formats: