New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

icedyn_adv_umx.F90 in NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_adv_umx.F90 @ 10315

Last change on this file since 10315 was 10315, checked in by clem, 5 years ago
still playing around with nonosc limiter to get a good ice thickness at the ice edge. There is still a problem when the ice concentration is very small
File size: 63.4 KB

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

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

Download in other formats: