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/UKMO/dev_r10448_bdyvol/src/ICE – NEMO

Context Navigation

source: NEMO/branches/UKMO/dev_r10448_bdyvol/src/ICE/icedyn_adv_umx.F90 @ 10455

Last change on this file since 10455 was 10446, checked in by clem, 5 years ago
clean the routine icedyn_adv_umx
Property svn:keywords set to `Id`
File size: 68.0 KB

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

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

Download in other formats: