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

Context Navigation

source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_adv_umx.F90 @ 10419

Last change on this file since 10419 was 10419, checked in by smasson, 5 years ago
dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10418, see #2133
Property svn:keywords set to `Id`
File size: 94.2 KB

Rev	Line
[8586]	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
[9604]	7	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
[8586]	8	!!----------------------------------------------------------------------
[9570]	9	#if defined key_si3
[8586]	10	!!----------------------------------------------------------------------
[9570]	11	!! 'key_si3' SI3 sea-ice model
[8586]	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 : ???
[10419]	16	!! nonosc : compute monotonic tracer fluxes by a non-oscillatory algorithm
[8586]	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
[10419]	22	USE icevar ! sea-ice: operations
[8586]	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
[10419]	37	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: amaxu, amaxv
	38
	39	! advect H all the way (and get V=H*A at the end)
	40	LOGICAL :: ll_thickness = .FALSE.
	41
	42	! look for 9 points around in nonosc limiter
	43	LOGICAL :: ll_9points = .FALSE. ! false=better h?
	44
	45	! use HgradU in nonosc limiter
	46	LOGICAL :: ll_HgradU = .TRUE. ! no effect?
	47
	48	! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme
	49	LOGICAL :: ll_neg = .TRUE. ! keep TRUE
	50
	51	! limit the fluxes
	52	LOGICAL :: ll_zeroup1 = .FALSE. ! false ok if Hbig otherwise needed for 2D sinon on a des valeurs de H trop fortes !!
	53	LOGICAL :: ll_zeroup2 = .FALSE. ! false ok for 1D, 2D, 3D
	54	LOGICAL :: ll_zeroup4 = .FALSE. ! false ok for 1D, 2D, 3D
	55	LOGICAL :: ll_zeroup5 = .FALSE. ! false ok for 1D, 2D
	56
	57	! fluxes that are limited are uH, then (uH)*(ua/u) is used for V (only for nonosc)
	58	LOGICAL :: ll_clem = .TRUE. ! simpler than rachid and works
	59
	60	! First advect H as H=Hdiv(u), then use H for H(n+1)=H(n)-div(uH*)
	61	LOGICAL :: ll_gurvan = .FALSE. ! must be false for 1D case !!
	62
	63	! First guess as div(uH) (-true-) or Hdiv(u)+ugradH (-false-)
	64	LOGICAL :: ll_1stguess_clem = .FALSE. ! better negative values but less good h
	65
	66	! advect (or not) open water. If not, retrieve it from advection of A
	67	LOGICAL :: ll_ADVopw = .FALSE. !
	68
	69	! alternate directions for upstream
	70	LOGICAL :: ll_upsxy = .TRUE.
	71
	72	! alternate directions for high order
	73	LOGICAL :: ll_hoxy = .TRUE.
	74
	75	! prelimiter: use it to avoid overshoot in H
	76	LOGICAL :: ll_prelimiter_zalesak = .TRUE. ! from: Zalesak(1979) eq. 14 => true is better for 1D but false is better in 3D (for h and negative values) => pb in x-y?
	77	LOGICAL :: ll_prelimiter_devore = .FALSE. ! from: Devore eq. 11 (worth than zalesak)
	78
	79	! iterate on the limiter (only nonosc)
	80	LOGICAL :: ll_limiter_it2 = .FALSE.
	81
	82
[8586]	83	!! * Substitutions
	84	# include "vectopt_loop_substitute.h90"
	85	!!----------------------------------------------------------------------
[9598]	86	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[10069]	87	!! $Id$
[10419]	88	!! Software governed by the CeCILL licence (./LICENSE)
[8586]	89	!!----------------------------------------------------------------------
	90	CONTAINS
	91
[10419]	92	SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, &
	93	& pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
[8586]	94	!!----------------------------------------------------------------------
	95	!! * ROUTINE ice_dyn_adv_umx *
	96	!!
	97	!! ** Purpose : Compute the now trend due to total advection of
	98	!! tracers and add it to the general trend of tracer equations
	99	!! using an "Ultimate-Macho" scheme
	100	!!
	101	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
	102	!!----------------------------------------------------------------------
[10419]	103	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
[8586]	104	INTEGER , INTENT(in ) :: kt ! time step
	105	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity
	106	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pv_ice ! ice j-velocity
	107	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area
	108	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume
	109	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_s ! snw volume
	110	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psv_i ! salt content
	111	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content
	112	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration
	113	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
	114	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume
	115	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content
	116	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content
	117	!
	118	INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices
[10419]	119	INTEGER :: icycle ! number of sub-timestep for the advection
	120	REAL(wp) :: zamsk ! 1 if advection of concentration, 0 if advection of other tracers
	121	REAL(wp) :: zdt
	122	REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! send zcflnow and receive zcflprv
	123	REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx
	124	REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2
[10180]	125
[10419]	126
	127
	128	REAL(wp), DIMENSION(jpi,jpj) :: zcu_box, zcv_box
	129	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zua_ho, zva_ho
	130	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_ai, z1_aip
	131	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhvar
	132
[8586]	133	!!----------------------------------------------------------------------
	134	!
	135	IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
	136	!
	137	!
[10292]	138	! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !
	139	! When needed, the advection split is applied at the next time-step in order to avoid blocking global comm.
	140	! ...this should not affect too much the stability... Was ok on the tests we did...
[10397]	141	zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
	142	zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
[10292]	143
[10397]	144	! non-blocking global communication send zcflnow and receive zcflprv
	145	CALL mpp_delay_max( 'icedyn_adv_umx', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 )
[8586]	146
[10419]	147	IF( zcflprv(1) > .5 ) THEN ; icycle = 2
	148	ELSE ; icycle = 1
[8586]	149	ENDIF
[10419]	150
	151	zdt = rdt_ice / REAL(icycle)
[8586]	152
	153	! --- transport --- !
	154	zudy(:,:) = pu_ice(:,:) * e2u(:,:)
	155	zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
	156
	157	! --- define velocity for advection: u*grad(H) --- !
	158	DO jj = 2, jpjm1
	159	DO ji = fs_2, fs_jpim1
	160	IF ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN ; zcu_box(ji,jj) = 0._wp
	161	ELSEIF( pu_ice(ji,jj) > 0._wp ) THEN ; zcu_box(ji,jj) = pu_ice(ji-1,jj)
	162	ELSE ; zcu_box(ji,jj) = pu_ice(ji ,jj)
	163	ENDIF
	164
	165	IF ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN ; zcv_box(ji,jj) = 0._wp
	166	ELSEIF( pv_ice(ji,jj) > 0._wp ) THEN ; zcv_box(ji,jj) = pv_ice(ji,jj-1)
	167	ELSE ; zcv_box(ji,jj) = pv_ice(ji,jj )
	168	ENDIF
	169	END DO
	170	END DO
	171
[10419]	172	IF( ll_zeroup2 ) THEN
	173	IF(.NOT. ALLOCATED(amaxu)) ALLOCATE(amaxu (jpi,jpj,jpl))
	174	IF(.NOT. ALLOCATED(amaxv)) ALLOCATE(amaxv (jpi,jpj,jpl))
	175	ENDIF
[8586]	176	!---------------!
	177	!== advection ==!
	178	!---------------!
[10419]	179	DO jt = 1, icycle
	180
	181	!!$ IF( ll_ADVopw ) THEN
	182	!!$ zamsk = 1._wp
	183	!!$ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) ) ! Open water area
	184	!!$ zamsk = 0._wp
	185	!!$ ELSE
	186	zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
	187	!!$ ENDIF
[10386]	188
[10419]	189	WHERE( pa_i(:,:,:) >= epsi20 ) ; z1_ai(:,:,:) = 1._wp / pa_i(:,:,:)
	190	ELSEWHERE ; z1_ai(:,:,:) = 0.
	191	END WHERE
	192	!
	193	WHERE( pa_ip(:,:,:) >= epsi20 ) ; z1_aip(:,:,:) = 1._wp / pa_ip(:,:,:)
	194	ELSEWHERE ; z1_aip(:,:,:) = 0.
	195	END WHERE
	196	!
	197	IF( ll_zeroup2 ) THEN
	198	DO jl = 1, jpl
	199	DO jj = 2, jpjm1
	200	DO ji = fs_2, fs_jpim1
	201	amaxu(ji,jj,jl)=MAX( pa_i(ji,jj,jl), pa_i(ji,jj-1,jl), pa_i(ji,jj+1,jl), &
	202	& pa_i(ji+1,jj,jl), pa_i(ji+1,jj-1,jl), pa_i(ji+1,jj+1,jl) )
	203	amaxv(ji,jj,jl)=MAX( pa_i(ji,jj,jl), pa_i(ji-1,jj,jl), pa_i(ji+1,jj,jl), &
	204	& pa_i(ji,jj+1,jl), pa_i(ji-1,jj+1,jl), pa_i(ji+1,jj+1,jl) )
	205	END DO
	206	END DO
	207	END DO
	208	CALL lbc_lnk_multi('icedyn_adv_umx', amaxu, 'T', 1., amaxv, 'T', 1.)
	209	ENDIF
	210	!
	211	DO jl = 1, jpl
	212	zua_ho(:,:,jl) = zudy(:,:)
	213	zva_ho(:,:,jl) = zvdx(:,:)
	214	END DO
	215
	216	zamsk = 1._wp
	217	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
	218	zamsk = 0._wp
	219	!
	220	IF( ll_thickness ) THEN
	221	DO jl = 1, jpl
	222	zua_ho(:,:,jl) = zudy(:,:)
	223	zva_ho(:,:,jl) = zvdx(:,:)
	224	END DO
	225	ENDIF
	226	!
	227	zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
	228	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_i ) ! Ice volume
	229	IF( ll_thickness ) pv_i(:,:,:) = zhvar(:,:,:) * pa_i(:,:,:)
	230	!
	231	zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
	232	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_s ) ! Snw volume
	233	IF( ll_thickness ) pv_s(:,:,:) = zhvar(:,:,:) * pa_i(:,:,:)
	234	!
	235	zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:)
	236	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, psv_i ) ! Salt content
	237	!
	238	zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:)
	239	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, poa_i ) ! Age content
	240	!
[10180]	241	DO jk = 1, nlay_i
[10419]	242	zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:)
	243	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
[8586]	244	END DO
[10419]	245	!
[10180]	246	DO jk = 1, nlay_s
[10419]	247	zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:)
	248	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
[10180]	249	END DO
[10419]	250	!
[10180]	251	IF ( ln_pnd_H12 ) THEN
[10419]	252	!
	253	DO jl = 1, jpl
	254	zua_ho(:,:,jl) = zudy(:,:)
	255	zva_ho(:,:,jl) = zvdx(:,:)
	256	END DO
	257
	258	zamsk = 1._wp
	259	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
	260	zamsk = 0._wp
	261	!
	262	zhvar(:,:,:) = pv_ip(:,:,:) * z1_ai(:,:,:)
	263	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_ip ) ! mp volume
[10386]	264	ENDIF
[10419]	265	!
	266	!
	267	!!$ IF( .NOT. ll_ADVopw ) THEN
	268	zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
	269	DO jj = 2, jpjm1
	270	DO ji = fs_2, fs_jpim1
	271	pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & ! Open water area
	272	& - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )r1_e1e2t(ji,jj)zdt
	273	END DO
	274	END DO
	275	CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T', 1. )
	276	!!$ ENDIF
	277	!
[8586]	278	END DO
	279	!
	280	END SUBROUTINE ice_dyn_adv_umx
[9929]	281
[8586]	282
[10419]	283	SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho )
[8586]	284	!!----------------------------------------------------------------------
	285	!! * ROUTINE adv_umx *
	286	!!
	287	!! ** Purpose : Compute the now trend due to total advection of
	288	!! tracers and add it to the general trend of tracer equations
	289	!!
	290	!! ** Method : TVD scheme, i.e. 2nd order centered scheme with
	291	!! corrected flux (monotonic correction)
	292	!! note: - this advection scheme needs a leap-frog time scheme
	293	!!
	294	!! ** Action : - pt the after advective tracer
	295	!!----------------------------------------------------------------------
[10419]	296	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
	297	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
	298	INTEGER , INTENT(in ) :: jt ! number of sub-iteration
	299	INTEGER , INTENT(in ) :: kt ! number of iteration
	300	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
	301	REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pu , pv ! 2 ice velocity components => u*e2
	302	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: puc , pvc ! 2 ice velocity components => ue2 or ua*e2u
	303	REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pubox, pvbox ! upstream velocity
	304	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pt ! tracer field
	305	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptc ! tracer content field
	306	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out), OPTIONAL :: pua_ho, pva_ho ! high order u*a fluxes
[8586]	307	!
[10180]	308	INTEGER :: ji, jj, jl ! dummy loop indices
[8586]	309	REAL(wp) :: ztra ! local scalar
[10419]	310	INTEGER :: kn_limiter = 1 ! 1=nonosc ; 2=superbee ; 3=h3(rachid)
	311	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zfu_ho , zfv_ho , zt_u, zt_v, zpt
	312	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zfu_ups, zfv_ups, zt_ups ! only for nonosc
[8586]	313	!!----------------------------------------------------------------------
	314	!
[10419]	315	! upstream (_ups) advection with initial mass fluxes
	316	! ---------------------------------------------------
	317
	318	IF( ll_gurvan .AND. pamsk==0. ) THEN
	319	DO jl = 1, jpl
	320	DO jj = 2, jpjm1
	321	DO ji = fs_2, fs_jpim1
	322	pt(ji,jj,jl) = ( pt (ji,jj,jl) + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) &
	323	& + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) ) &
	324	& * tmask(ji,jj,1)
	325	END DO
[10180]	326	END DO
[8586]	327	END DO
[10419]	328	CALL lbc_lnk( 'icedyn_adv_umx', pt, 'T', 1. )
	329	ENDIF
	330
[8586]	331
[10419]	332	IF( .NOT. ll_upsxy ) THEN
	333
	334	! fluxes in both x-y directions
	335	DO jl = 1, jpl
[10180]	336	DO jj = 1, jpjm1
[10419]	337	DO ji = 1, fs_jpim1
	338	IF( ll_clem ) THEN
	339	zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl)
	340	zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl)
	341	ELSE
	342	zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * pt(ji+1,jj,jl)
	343	zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj+1,jl)
	344	ENDIF
[10180]	345	END DO
[8586]	346	END DO
	347	END DO
[10419]	348
	349	ELSE
[8586]	350	!
[10419]	351	IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==!
	352	! flux in x-direction
	353	DO jl = 1, jpl
	354	DO jj = 1, jpjm1
	355	DO ji = 1, fs_jpim1
	356	IF( ll_clem ) THEN
	357	zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl)
	358	ELSE
	359	zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * pt(ji+1,jj,jl)
	360	ENDIF
	361	END DO
	362	END DO
	363	END DO
	364
	365	! first guess of tracer content from u-flux
	366	DO jl = 1, jpl
	367	DO jj = 2, jpjm1
	368	DO ji = fs_2, fs_jpim1 ! vector opt.
	369	IF( ll_clem ) THEN
	370	IF( ll_gurvan ) THEN
	371	zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
	372	ELSE
	373	zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &
	374	& + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
	375	& ) * tmask(ji,jj,1)
	376	ENDIF
	377	ELSE
	378	zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) &
	379	& * tmask(ji,jj,1)
	380	ENDIF
	381	!! IF( ji==26 .AND. jj==86) THEN
	382	!! WRITE(numout,) '***********************'
	383	!! WRITE(numout,*) 'zpt upstream',zpt(ji,jj)
	384	!! ENDIF
	385	END DO
	386	END DO
	387	END DO
	388	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
	389	!
	390	! flux in y-direction
	391	DO jl = 1, jpl
	392	DO jj = 1, jpjm1
	393	DO ji = 1, fs_jpim1
	394	IF( ll_clem ) THEN
	395	zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl)
	396	ELSE
	397	zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * zpt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * zpt(ji,jj+1,jl)
	398	ENDIF
	399	END DO
	400	END DO
	401	END DO
[8586]	402	!
[10419]	403	ELSE !== even ice time step: adv_y then adv_x ==!
	404	! flux in y-direction
	405	DO jl = 1, jpl
	406	DO jj = 1, jpjm1
	407	DO ji = 1, fs_jpim1
	408	IF( ll_clem ) THEN
	409	zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl)
	410	ELSE
	411	zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj+1,jl)
	412	ENDIF
	413	END DO
[10345]	414	END DO
	415	END DO
[10419]	416
	417	! first guess of tracer content from v-flux
	418	DO jl = 1, jpl
	419	DO jj = 2, jpjm1
	420	DO ji = fs_2, fs_jpim1 ! vector opt.
	421	IF( ll_clem ) THEN
	422	IF( ll_gurvan ) THEN
	423	zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
	424	ELSE
	425	zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &
	426	& + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) &
	427	& * tmask(ji,jj,1)
	428	ENDIF
	429	ELSE
	430	zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) &
	431	& * tmask(ji,jj,1)
	432	ENDIF
	433	!! IF( ji==26 .AND. jj==86) THEN
	434	!! WRITE(numout,) '***********************'
	435	!! WRITE(numout,*) 'zpt upstream',zpt(ji,jj)
	436	!! ENDIF
	437	END DO
	438	END DO
	439	END DO
	440	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
	441	!
	442	! flux in x-direction
	443	DO jl = 1, jpl
	444	DO jj = 1, jpjm1
	445	DO ji = 1, fs_jpim1
	446	IF( ll_clem ) THEN
	447	zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl)
	448	ELSE
	449	zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * zpt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * zpt(ji+1,jj,jl)
	450	ENDIF
	451	END DO
	452	END DO
	453	END DO
	454	!
	455	ENDIF
	456
	457	ENDIF
	458
	459	IF( ll_clem .AND. kn_limiter /= 1 ) &
	460	& CALL ctl_stop( 'STOP', 'icedyn_adv_umx: ll_clem incompatible with limiters other than nonosc' )
	461
	462	IF( ll_zeroup2 ) THEN
	463	DO jl = 1, jpl
[10345]	464	DO jj = 1, jpjm1
[10419]	465	DO ji = 1, fs_jpim1 ! vector opt.
	466	IF( amaxu(ji,jj,jl) == 0._wp ) zfu_ups(ji,jj,jl) = 0._wp
	467	IF( amaxv(ji,jj,jl) == 0._wp ) zfv_ups(ji,jj,jl) = 0._wp
[10180]	468	END DO
[8586]	469	END DO
	470	END DO
[10419]	471	ENDIF
	472
	473	! guess after content field with upstream scheme
	474	DO jl = 1, jpl
[10345]	475	DO jj = 2, jpjm1
[10419]	476	DO ji = fs_2, fs_jpim1
	477	ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj ,jl) &
	478	& + zfv_ups(ji,jj,jl) - zfv_ups(ji ,jj-1,jl) ) * r1_e1e2t(ji,jj)
	479	IF( ll_clem ) THEN
	480	IF( ll_gurvan ) THEN
	481	zt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1)
	482	ELSE
	483	zt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + pdt * ztra + ( pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) &
	484	& + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) ) &
	485	& * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
	486	ENDIF
	487	ELSE
	488	zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1)
	489	ENDIF
	490	!! IF( ji==26 .AND. jj==86) THEN
	491	!! WRITE(numout,) '*************************'
	492	!! WRITE(numout,*) 'zt upstream',zt_ups(ji,jj)
	493	!! ENDIF
[10345]	494	END DO
	495	END DO
	496	END DO
[10419]	497	CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. )
	498
	499	! High order (_ho) fluxes
	500	! -----------------------
	501	SELECT CASE( kn_umx )
	502	!
	503	CASE ( 20 ) !== centered second order ==!
	504	!
	505	CALL cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho, &
	506	& zt_ups, zfu_ups, zfv_ups )
	507	!
	508	CASE ( 1:5 ) !== 1st to 5th order ULTIMATE-MACHO scheme ==!
	509	!
	510	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, &
	511	& zt_ups, zfu_ups, zfv_ups )
	512	!
	513	END SELECT
	514
	515	IF( ll_thickness ) THEN
	516	! final trend with corrected fluxes
	517	! ------------------------------------
	518	DO jl = 1, jpl
	519	DO jj = 2, jpjm1
	520	DO ji = fs_2, fs_jpim1
	521	IF( ll_gurvan ) THEN
	522	ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj)
	523	ELSE
	524	ztra = ( - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) &
	525	& + ( pt(ji,jj,jl) * ( pu(ji,jj) - pu(ji-1,jj) ) * (1.-pamsk) ) &
	526	& + ( pt(ji,jj,jl) * ( pv(ji,jj) - pv(ji,jj-1) ) * (1.-pamsk) ) ) * r1_e1e2t(ji,jj)
	527	ENDIF
	528	pt(ji,jj,jl) = ( pt(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1)
	529
	530	IF( pt(ji,jj,jl) < -epsi20 ) THEN
	531	WRITE(numout,*) 'T<0 ',pt(ji,jj,jl)
	532	ENDIF
	533
	534	IF( pt(ji,jj,jl) < 0._wp .AND. pt(ji,jj,jl) >= -epsi20 ) pt(ji,jj,jl) = 0._wp
	535
	536	!! IF( ji==26 .AND. jj==86) THEN
	537	!! WRITE(numout,*) 'zt high order',pt(ji,jj)
	538	!! ENDIF
	539	END DO
[10180]	540	END DO
[8586]	541	END DO
[10419]	542	CALL lbc_lnk( 'icedyn_adv_umx', pt, 'T', 1. )
	543	ENDIF
	544
	545	! Rachid trick
	546	! ------------
	547	IF( ll_clem ) THEN
	548	IF( pamsk == 0. ) THEN
	549	DO jl = 1, jpl
	550	DO jj = 1, jpjm1
	551	DO ji = 1, fs_jpim1
	552	IF( ABS( puc(ji,jj,jl) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN
	553	zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj)
	554	zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj)
	555	ELSE
	556	zfu_ho (ji,jj,jl) = 0._wp
	557	zfu_ups(ji,jj,jl) = 0._wp
	558	ENDIF
	559	!
	560	IF( ABS( pvc(ji,jj,jl) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN
	561	zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj)
	562	zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj)
	563	ELSE
	564	zfv_ho (ji,jj,jl) = 0._wp
	565	zfv_ups(ji,jj,jl) = 0._wp
	566	ENDIF
	567	END DO
	568	END DO
	569	END DO
	570	ENDIF
	571	ENDIF
[10180]	572
[10419]	573	IF( ll_zeroup5 ) THEN
	574	DO jl = 1, jpl
	575	DO jj = 2, jpjm1
	576	DO ji = 2, fs_jpim1 ! vector opt.
	577	zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &
	578	& - ( zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
	579	IF( zpt(ji,jj,jl) < 0. ) THEN
	580	zfu_ho(ji ,jj,jl) = zfu_ups(ji ,jj,jl)
	581	zfu_ho(ji-1,jj,jl) = zfu_ups(ji-1,jj,jl)
	582	zfv_ho(ji ,jj,jl) = zfv_ups(ji ,jj,jl)
	583	zfv_ho(ji,jj-1,jl) = zfv_ups(ji,jj-1,jl)
	584	ENDIF
	585	END DO
[10180]	586	END DO
[8586]	587	END DO
[10419]	588	CALL lbc_lnk_multi( 'icedyn_adv_umx', zfu_ho, 'U', -1., zfv_ho, 'V', -1. )
	589	ENDIF
	590
	591	! output high order fluxes u*a
	592	! ----------------------------
	593	IF( PRESENT( pua_ho ) ) THEN
	594	DO jl = 1, jpl
	595	DO jj = 1, jpjm1
	596	DO ji = 1, fs_jpim1
	597	pua_ho(ji,jj,jl) = zfu_ho(ji,jj,jl)
	598	pva_ho(ji,jj,jl) = zfv_ho(ji,jj,jl)
	599	END DO
	600	END DO
	601	END DO
	602	ENDIF
	603
	604
	605	IF( .NOT.ll_thickness ) THEN
	606	! final trend with corrected fluxes
	607	! ------------------------------------
	608	DO jl = 1, jpl
	609	DO jj = 2, jpjm1
	610	DO ji = fs_2, fs_jpim1
	611	ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj) * pdt
	612
	613	ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra ) * tmask(ji,jj,1)
	614
	615	!! IF( ji==26 .AND. jj==86) THEN
	616	!! WRITE(numout,*) 'ztc high order',ptc(ji,jj)
	617	!! ENDIF
	618
	619	END DO
	620	END DO
	621	END DO
	622	CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T', 1. )
	623	ENDIF
	624
[8586]	625	!
	626	END SUBROUTINE adv_umx
	627
[10419]	628	SUBROUTINE cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, &
	629	& pt_ups, pfu_ups, pfv_ups )
[8586]	630	!!---------------------------------------------------------------------
[10419]	631	!! * ROUTINE macho *
[8586]	632	!!
	633	!! ** Purpose : compute
	634	!!
	635	!! ** Method : ... ???
	636	!! TIM = transient interpolation Modeling
	637	!!
	638	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
	639	!!----------------------------------------------------------------------
[10419]	640	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
	641	INTEGER , INTENT(in ) :: kn_limiter ! limiter
	642	INTEGER , INTENT(in ) :: jt ! number of sub-iteration
	643	INTEGER , INTENT(in ) :: kt ! number of iteration
	644	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
	645	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt ! tracer fields
	646	REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pu, pv ! 2 ice velocity components
	647	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: puc, pvc ! 2 ice velocity * A components
	648	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptc ! tracer content at before time step
	649	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes
	650	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt_ups ! upstream guess of tracer content
	651	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pfu_ups, pfv_ups ! upstream fluxes
[8586]	652	!
[10180]	653	INTEGER :: ji, jj, jl ! dummy loop indices
[10419]	654	LOGICAL :: ll_xy = .TRUE.
	655	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zzt
[8586]	656	!!----------------------------------------------------------------------
	657	!
[10419]	658	IF( .NOT.ll_xy ) THEN !-- no alternate directions --!
[8586]	659	!
[10419]	660	DO jl = 1, jpl
	661	DO jj = 1, jpjm1
	662	DO ji = 1, fs_jpim1
	663	IF( ll_clem ) THEN
	664	pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
	665	pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
	666	ELSE
	667	pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
	668	pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
	669	ENDIF
[10180]	670	END DO
[8586]	671	END DO
	672	END DO
[10419]	673	IF ( kn_limiter == 1 ) THEN
	674	IF( ll_clem ) THEN
	675	CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
	676	ELSE
	677	CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
	678	ENDIF
	679	ELSEIF( kn_limiter == 2 ) THEN
	680	CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
	681	CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
	682	ELSEIF( kn_limiter == 3 ) THEN
	683	CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
	684	CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
	685	ENDIF
[8586]	686	!
[10419]	687	ELSE !-- alternate directions --!
[8586]	688	!
[10419]	689	IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==!
	690	!
	691	! flux in x-direction
	692	DO jl = 1, jpl
	693	DO jj = 1, jpjm1
	694	DO ji = 1, fs_jpim1
	695	IF( ll_clem ) THEN
	696	pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
	697	ELSE
	698	pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
	699	ENDIF
	700	END DO
	701	END DO
	702	END DO
	703	IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
	704	IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
	705
	706	! first guess of tracer content from u-flux
	707	DO jl = 1, jpl
	708	DO jj = 2, jpjm1
	709	DO ji = fs_2, fs_jpim1 ! vector opt.
	710	IF( ll_clem ) THEN
	711	zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &
	712	& + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) &
	713	& * tmask(ji,jj,1)
	714	ELSE
	715	zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
	716	ENDIF
	717	END DO
	718	END DO
	719	END DO
	720	CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
	721
	722	! flux in y-direction
	723	DO jl = 1, jpl
	724	DO jj = 1, jpjm1
	725	DO ji = 1, fs_jpim1
	726	IF( ll_clem ) THEN
	727	pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji,jj+1,jl) )
	728	ELSE
	729	pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( zzt(ji,jj,jl) + zzt(ji,jj+1,jl) )
	730	ENDIF
	731	END DO
	732	END DO
	733	END DO
	734	IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
	735	IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
	736
	737	ELSE !== even ice time step: adv_y then adv_x ==!
	738	!
	739	! flux in y-direction
	740	DO jl = 1, jpl
	741	DO jj = 1, jpjm1
	742	DO ji = 1, fs_jpim1
	743	IF( ll_clem ) THEN
	744	pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
	745	ELSE
	746	pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
	747	ENDIF
	748	END DO
	749	END DO
	750	END DO
	751	IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
	752	IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
	753	!
	754	! first guess of tracer content from v-flux
	755	DO jl = 1, jpl
	756	DO jj = 2, jpjm1
	757	DO ji = fs_2, fs_jpim1 ! vector opt.
	758	IF( ll_clem ) THEN
	759	zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &
	760	& + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) &
	761	& * tmask(ji,jj,1)
	762	ELSE
	763	zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
	764	ENDIF
	765	END DO
	766	END DO
	767	END DO
	768	CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
	769	!
	770	! flux in x-direction
	771	DO jl = 1, jpl
	772	DO jj = 1, jpjm1
	773	DO ji = 1, fs_jpim1
	774	IF( ll_clem ) THEN
	775	pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji+1,jj,jl) )
	776	ELSE
	777	pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( zzt(ji,jj,jl) + zzt(ji+1,jj,jl) )
	778	ENDIF
	779	END DO
	780	END DO
	781	END DO
	782	IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
	783	IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
	784
	785	ENDIF
	786	IF( ll_clem ) THEN
	787	IF( kn_limiter == 1 ) CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
	788	ELSE
	789	IF( kn_limiter == 1 ) CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
	790	ENDIF
	791
	792	ENDIF
	793
	794	END SUBROUTINE cen2
	795
	796
	797	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, &
	798	& pt_ups, pfu_ups, pfv_ups )
	799	!!---------------------------------------------------------------------
	800	!! * ROUTINE macho *
	801	!!
	802	!! ** Purpose : compute
	803	!!
	804	!! ** Method : ... ???
	805	!! TIM = transient interpolation Modeling
	806	!!
	807	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
	808	!!----------------------------------------------------------------------
	809	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
	810	INTEGER , INTENT(in ) :: kn_limiter ! limiter
	811	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
	812	INTEGER , INTENT(in ) :: jt ! number of sub-iteration
	813	INTEGER , INTENT(in ) :: kt ! number of iteration
	814	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
	815	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt ! tracer fields
	816	REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pu, pv ! 2 ice velocity components
	817	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: puc, pvc ! 2 ice velocity * A components
	818	REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pubox, pvbox ! upstream velocity
	819	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptc ! tracer content at before time step
	820	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_u, pt_v ! tracer at u- and v-points
	821	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes
	822	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt_ups ! upstream guess of tracer content
	823	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pfu_ups, pfv_ups ! upstream fluxes
	824	!
	825	INTEGER :: ji, jj, jl ! dummy loop indices
	826	REAL(wp) :: ztra
	827	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zzt, zzfu_ho, zzfv_ho
	828	!!----------------------------------------------------------------------
	829	!
	830	IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==!
[8586]	831	!
[10419]	832	! !-- ultimate interpolation of pt at u-point --!
	833	CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
	834	! !-- limiter in x --!
	835	IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
	836	IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
	837	! !-- advective form update in zzt --!
	838
	839	IF( ll_1stguess_clem ) THEN
	840
	841	! first guess of tracer content from u-flux
	842	DO jl = 1, jpl
	843	DO jj = 2, jpjm1
	844	DO ji = fs_2, fs_jpim1 ! vector opt.
	845	IF( ll_clem ) THEN
	846	IF( ll_gurvan ) THEN
	847	zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
	848	ELSE
	849	zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &
	850	& + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
	851	& ) * tmask(ji,jj,1)
	852	ENDIF
	853	ELSE
	854	zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
	855	ENDIF
	856	END DO
	857	END DO
	858	END DO
	859	CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
	860
	861	ELSE
	862
	863	DO jl = 1, jpl
	864	DO jj = 2, jpjm1
	865	DO ji = fs_2, fs_jpim1 ! vector opt.
	866	IF( ll_gurvan ) THEN
	867	zzt(ji,jj,jl) = pt(ji,jj,jl) - pubox(ji,jj ) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj) &
	868	& - pt (ji,jj,jl) * pdt * ( pu (ji,jj) - pu (ji-1,jj) ) * r1_e1e2t(ji,jj)
	869	ELSE
	870	zzt(ji,jj,jl) = pt(ji,jj,jl) - pubox(ji,jj ) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj) &
	871	& - pt (ji,jj,jl) * pdt * ( pu (ji,jj) - pu (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk
	872	ENDIF
	873	zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1)
	874	END DO
	875	END DO
	876	END DO
	877	CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
	878	ENDIF
[8586]	879	!
[10419]	880	! !-- ultimate interpolation of pt at v-point --!
	881	IF( ll_hoxy ) THEN
	882	CALL ultimate_y( kn_umx, pdt, zzt, pv, pvc, pt_v, pfv_ho )
	883	ELSE
	884	CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
	885	ENDIF
	886	! !-- limiter in y --!
	887	IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
	888	IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
	889	!
	890	!
	891	ELSE !== even ice time step: adv_y then adv_x ==!
	892	!
	893	! !-- ultimate interpolation of pt at v-point --!
	894	CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
	895	! !-- limiter in y --!
	896	IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
	897	IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
	898	! !-- advective form update in zzt --!
	899	IF( ll_1stguess_clem ) THEN
	900
	901	! first guess of tracer content from v-flux
	902	DO jl = 1, jpl
	903	DO jj = 2, jpjm1
	904	DO ji = fs_2, fs_jpim1 ! vector opt.
	905	IF( ll_clem ) THEN
	906	IF( ll_gurvan ) THEN
	907	zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
	908	ELSE
	909	zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &
	910	& + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
	911	& ) * tmask(ji,jj,1)
	912	ENDIF
	913	ELSE
	914	zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) &
	915	& * tmask(ji,jj,1)
	916	ENDIF
	917	END DO
[10180]	918	END DO
[8586]	919	END DO
[10419]	920	CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
	921
	922	ELSE
	923
	924	DO jl = 1, jpl
	925	DO jj = 2, jpjm1
	926	DO ji = fs_2, fs_jpim1
	927	IF( ll_gurvan ) THEN
	928	zzt(ji,jj,jl) = pt(ji,jj,jl) - pvbox(ji,jj ) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj) &
	929	& - pt (ji,jj,jl) * pdt * ( pv (ji,jj) - pv (ji,jj-1) ) * r1_e1e2t(ji,jj)
	930	ELSE
	931	zzt(ji,jj,jl) = pt(ji,jj,jl) - pvbox(ji,jj ) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj) &
	932	& - pt (ji,jj,jl) * pdt * ( pv (ji,jj) - pv (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk
	933	ENDIF
	934	zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1)
	935	END DO
	936	END DO
	937	END DO
	938	CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
	939	ENDIF
[8586]	940	!
[10419]	941	! !-- ultimate interpolation of pt at u-point --!
	942	IF( ll_hoxy ) THEN
	943	CALL ultimate_x( kn_umx, pdt, zzt, pu, puc, pt_u, pfu_ho )
	944	ELSE
	945	CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
	946	ENDIF
	947	! !-- limiter in x --!
	948	IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
	949	IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
	950	!
	951	!
	952	ENDIF
	953
	954
	955	IF( kn_limiter == 1 ) THEN
	956	IF( .NOT. ll_limiter_it2 ) THEN
	957	IF( ll_clem ) THEN
	958	CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
	959	ELSE
	960	CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
	961	ENDIF
	962	ELSE
	963	zzfu_ho(:,:,:) = pfu_ho(:,:,:)
	964	zzfv_ho(:,:,:) = pfv_ho(:,:,:)
	965	! 1st iteration of nonosc (limit the flux with the upstream solution)
	966	IF( ll_clem ) THEN
	967	CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho )
	968	ELSE
	969	CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho )
	970	ENDIF
	971	! guess after content field with high order
	972	DO jl = 1, jpl
	973	DO jj = 2, jpjm1
	974	DO ji = fs_2, fs_jpim1
	975	ztra = - ( zzfu_ho(ji,jj,jl) - zzfu_ho(ji-1,jj,jl) + zzfv_ho(ji,jj,jl) - zzfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj)
	976	zzt(ji,jj,jl) = ( ptc(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1)
	977	END DO
	978	END DO
	979	END DO
	980	CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
	981	! 2nd iteration of nonosc (limit the flux with the limited high order solution)
	982	IF( ll_clem ) THEN
	983	CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho )
	984	ELSE
	985	CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho )
	986	ENDIF
	987	ENDIF
	988	ENDIF
[8586]	989	!
	990	END SUBROUTINE macho
	991
	992
[10419]	993	SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
[8586]	994	!!---------------------------------------------------------------------
	995	!! * ROUTINE ultimate_x *
	996	!!
	997	!! ** Purpose : compute
	998	!!
	999	!! ** Method : ... ???
	1000	!! TIM = transient interpolation Modeling
	1001	!!
	1002	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
	1003	!!----------------------------------------------------------------------
[10419]	1004	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
	1005	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
	1006	REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pu ! ice i-velocity component
	1007	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: puc ! ice i-velocity * A component
	1008	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt ! tracer fields
	1009	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_u ! tracer at u-point
	1010	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho ! high order flux
[8586]	1011	!
[10419]	1012	INTEGER :: ji, jj, jl ! dummy loop indices
[8586]	1013	REAL(wp) :: zcu, zdx2, zdx4 ! - -
[10419]	1014	REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztu1, ztu2, ztu3, ztu4
[8586]	1015	!!----------------------------------------------------------------------
	1016	!
	1017	! !-- Laplacian in i-direction --!
[10419]	1018	DO jl = 1, jpl
[10180]	1019	DO jj = 2, jpjm1 ! First derivative (gradient)
	1020	DO ji = 1, fs_jpim1
[10419]	1021	ztu1(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
[10180]	1022	END DO
	1023	! ! Second derivative (Laplacian)
	1024	DO ji = fs_2, fs_jpim1
[10419]	1025	ztu2(ji,jj,jl) = ( ztu1(ji,jj,jl) - ztu1(ji-1,jj,jl) ) * r1_e1t(ji,jj)
[10180]	1026	END DO
[8586]	1027	END DO
	1028	END DO
[10170]	1029	CALL lbc_lnk( 'icedyn_adv_umx', ztu2, 'T', 1. )
[8586]	1030	!
	1031	! !-- BiLaplacian in i-direction --!
[10419]	1032	DO jl = 1, jpl
[10180]	1033	DO jj = 2, jpjm1 ! Third derivative
	1034	DO ji = 1, fs_jpim1
[10419]	1035	ztu3(ji,jj,jl) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
[10180]	1036	END DO
	1037	! ! Fourth derivative
	1038	DO ji = fs_2, fs_jpim1
[10419]	1039	ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj)
[10180]	1040	END DO
[8586]	1041	END DO
	1042	END DO
[10170]	1043	CALL lbc_lnk( 'icedyn_adv_umx', ztu4, 'T', 1. )
[8586]	1044	!
	1045	!
[10419]	1046	SELECT CASE (kn_umx )
[8586]	1047	!
	1048	CASE( 1 ) !== 1st order central TIM ==! (Eq. 21)
	1049	!
[10419]	1050	DO jl = 1, jpl
	1051	DO jj = 1, jpjm1
[10180]	1052	DO ji = 1, fs_jpim1 ! vector opt.
[10419]	1053	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) &
	1054	& - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
[10180]	1055	END DO
[8586]	1056	END DO
	1057	END DO
	1058	!
	1059	CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23)
	1060	!
[10419]	1061	DO jl = 1, jpl
	1062	DO jj = 1, jpjm1
[10180]	1063	DO ji = 1, fs_jpim1 ! vector opt.
[10419]	1064	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
	1065	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) &
	1066	& - zcu * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
[10180]	1067	END DO
[8586]	1068	END DO
	1069	END DO
	1070	!
	1071	CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24)
	1072	!
[10419]	1073	DO jl = 1, jpl
	1074	DO jj = 1, jpjm1
[10180]	1075	DO ji = 1, fs_jpim1 ! vector opt.
[10419]	1076	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
[10180]	1077	zdx2 = e1u(ji,jj) * e1u(ji,jj)
[8586]	1078	!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
[10419]	1079	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
	1080	& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
	1081	& + z1_6 * zdx2 * ( zcuzcu - 1._wp ) ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
	1082	& - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
[10180]	1083	END DO
[8586]	1084	END DO
	1085	END DO
	1086	!
	1087	CASE( 4 ) !== 4th order central TIM ==! (Eq. 27)
	1088	!
[10419]	1089	DO jl = 1, jpl
	1090	DO jj = 1, jpjm1
[10180]	1091	DO ji = 1, fs_jpim1 ! vector opt.
[10419]	1092	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
[10180]	1093	zdx2 = e1u(ji,jj) * e1u(ji,jj)
[8586]	1094	!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
[10419]	1095	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
	1096	& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
	1097	& + z1_6 * zdx2 * ( zcuzcu - 1._wp ) ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
	1098	& - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
[10180]	1099	END DO
[8586]	1100	END DO
	1101	END DO
	1102	!
	1103	CASE( 5 ) !== 5th order central TIM ==! (Eq. 29)
	1104	!
[10419]	1105	DO jl = 1, jpl
	1106	DO jj = 1, jpjm1
[10180]	1107	DO ji = 1, fs_jpim1 ! vector opt.
[10419]	1108	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
[10180]	1109	zdx2 = e1u(ji,jj) * e1u(ji,jj)
[10419]	1110	!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
[10180]	1111	zdx4 = zdx2 * zdx2
[10419]	1112	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
	1113	& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
	1114	& + z1_6 * zdx2 * ( zcuzcu - 1._wp ) ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
	1115	& - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) &
	1116	& + z1_120 * zdx4 * ( zcuzcu - 1._wp ) ( zcuzcu - 4._wp ) ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl) &
	1117	& - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) )
[10180]	1118	END DO
[8586]	1119	END DO
	1120	END DO
	1121	!
	1122	END SELECT
[10419]	1123	! !-- High order flux in i-direction --!
	1124	IF( ll_neg ) THEN
	1125	DO jl = 1, jpl
	1126	DO jj = 1, jpjm1
	1127	DO ji = 1, fs_jpim1
	1128	IF( pt_u(ji,jj,jl) < 0._wp ) THEN
	1129	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) &
	1130	& - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
	1131	ENDIF
	1132	END DO
	1133	END DO
	1134	END DO
	1135	ENDIF
	1136
	1137	DO jl = 1, jpl
	1138	DO jj = 1, jpjm1
	1139	DO ji = 1, fs_jpim1 ! vector opt.
	1140	IF( ll_clem ) THEN
	1141	pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl)
	1142	ELSE
	1143	pfu_ho(ji,jj,jl) = puc(ji,jj,jl) * pt_u(ji,jj,jl)
	1144	ENDIF
	1145	END DO
	1146	END DO
	1147	END DO
[8586]	1148	!
	1149	END SUBROUTINE ultimate_x
	1150
	1151
[10419]	1152	SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
[8586]	1153	!!---------------------------------------------------------------------
	1154	!! * ROUTINE ultimate_y *
	1155	!!
	1156	!! ** Purpose : compute
	1157	!!
	1158	!! ** Method : ... ???
	1159	!! TIM = transient interpolation Modeling
	1160	!!
	1161	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
	1162	!!----------------------------------------------------------------------
[10419]	1163	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
	1164	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
	1165	REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pv ! ice j-velocity component
	1166	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pvc ! ice j-velocity*A component
	1167	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt ! tracer fields
	1168	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_v ! tracer at v-point
	1169	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfv_ho ! high order flux
[8586]	1170	!
[10180]	1171	INTEGER :: ji, jj, jl ! dummy loop indices
[8586]	1172	REAL(wp) :: zcv, zdy2, zdy4 ! - -
[10419]	1173	REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztv1, ztv2, ztv3, ztv4
[8586]	1174	!!----------------------------------------------------------------------
	1175	!
	1176	! !-- Laplacian in j-direction --!
[10419]	1177	DO jl = 1, jpl
[10180]	1178	DO jj = 1, jpjm1 ! First derivative (gradient)
	1179	DO ji = fs_2, fs_jpim1
[10419]	1180	ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
[10180]	1181	END DO
[8586]	1182	END DO
[10180]	1183	DO jj = 2, jpjm1 ! Second derivative (Laplacian)
	1184	DO ji = fs_2, fs_jpim1
[10419]	1185	ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj)
[10180]	1186	END DO
[8586]	1187	END DO
	1188	END DO
[10170]	1189	CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1. )
[8586]	1190	!
	1191	! !-- BiLaplacian in j-direction --!
[10419]	1192	DO jl = 1, jpl
[10180]	1193	DO jj = 1, jpjm1 ! First derivative
	1194	DO ji = fs_2, fs_jpim1
[10419]	1195	ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
[10180]	1196	END DO
[8586]	1197	END DO
[10180]	1198	DO jj = 2, jpjm1 ! Second derivative
	1199	DO ji = fs_2, fs_jpim1
[10419]	1200	ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj)
[10180]	1201	END DO
[8586]	1202	END DO
	1203	END DO
[10170]	1204	CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1. )
[8586]	1205	!
	1206	!
[10419]	1207	SELECT CASE (kn_umx )
	1208	!
[8586]	1209	CASE( 1 ) !== 1st order central TIM ==! (Eq. 21)
[10419]	1210	DO jl = 1, jpl
[10180]	1211	DO jj = 1, jpjm1
[10419]	1212	DO ji = 1, fs_jpim1
	1213	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) &
	1214	& - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
[10180]	1215	END DO
[8586]	1216	END DO
	1217	END DO
	1218	!
	1219	CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23)
[10419]	1220	DO jl = 1, jpl
[10180]	1221	DO jj = 1, jpjm1
[10419]	1222	DO ji = 1, fs_jpim1
	1223	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
[10180]	1224	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) &
	1225	& - zcv * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
	1226	END DO
[8586]	1227	END DO
	1228	END DO
[10170]	1229	CALL lbc_lnk( 'icedyn_adv_umx', pt_v, 'V', 1. )
[8586]	1230	!
	1231	CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24)
[10419]	1232	DO jl = 1, jpl
[10180]	1233	DO jj = 1, jpjm1
[10419]	1234	DO ji = 1, fs_jpim1
	1235	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
[10180]	1236	zdy2 = e2v(ji,jj) * e2v(ji,jj)
[8586]	1237	!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
[10180]	1238	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
	1239	& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
	1240	& + z1_6 * zdy2 * ( zcvzcv - 1._wp ) ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
	1241	& - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
	1242	END DO
[8586]	1243	END DO
	1244	END DO
	1245	!
	1246	CASE( 4 ) !== 4th order central TIM ==! (Eq. 27)
[10419]	1247	DO jl = 1, jpl
[10180]	1248	DO jj = 1, jpjm1
[10419]	1249	DO ji = 1, fs_jpim1
	1250	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
[10180]	1251	zdy2 = e2v(ji,jj) * e2v(ji,jj)
[8586]	1252	!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
[10180]	1253	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
	1254	& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
	1255	& + z1_6 * zdy2 * ( zcvzcv - 1._wp ) ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
	1256	& - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
	1257	END DO
[8586]	1258	END DO
	1259	END DO
	1260	!
	1261	CASE( 5 ) !== 5th order central TIM ==! (Eq. 29)
[10419]	1262	DO jl = 1, jpl
[10180]	1263	DO jj = 1, jpjm1
[10419]	1264	DO ji = 1, fs_jpim1
	1265	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
[10180]	1266	zdy2 = e2v(ji,jj) * e2v(ji,jj)
[8586]	1267	!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
[10180]	1268	zdy4 = zdy2 * zdy2
	1269	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
	1270	& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
	1271	& + z1_6 * zdy2 * ( zcvzcv - 1._wp ) ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
	1272	& - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) &
	1273	& + z1_120 * zdy4 * ( zcvzcv - 1._wp ) ( zcvzcv - 4._wp ) ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl) &
	1274	& - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) )
	1275	END DO
[8586]	1276	END DO
	1277	END DO
	1278	!
	1279	END SELECT
[10419]	1280	! !-- High order flux in j-direction --!
	1281	IF( ll_neg ) THEN
	1282	DO jl = 1, jpl
	1283	DO jj = 1, jpjm1
	1284	DO ji = 1, fs_jpim1
	1285	IF( pt_v(ji,jj,jl) < 0._wp ) THEN
	1286	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) &
	1287	& - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
	1288	ENDIF
	1289	END DO
	1290	END DO
	1291	END DO
	1292	ENDIF
	1293
	1294	DO jl = 1, jpl
	1295	DO jj = 1, jpjm1
	1296	DO ji = 1, fs_jpim1 ! vector opt.
	1297	IF( ll_clem ) THEN
	1298	pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl)
	1299	ELSE
	1300	pfv_ho(ji,jj,jl) = pvc(ji,jj,jl) * pt_v(ji,jj,jl)
	1301	ENDIF
	1302	END DO
	1303	END DO
	1304	END DO
[8586]	1305	!
	1306	END SUBROUTINE ultimate_y
[10419]	1307
	1308
	1309	SUBROUTINE nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho )
[8586]	1310	!!---------------------------------------------------------------------
	1311	!! * ROUTINE nonosc *
	1312	!!
[10419]	1313	!! ** Purpose : compute monotonic tracer fluxes from the upstream
[8586]	1314	!! scheme and the before field by a nonoscillatory algorithm
	1315	!!
	1316	!! ** Method : ... ???
[10419]	1317	!! warning : pt and pt_low must be masked, but the boundaries
[8586]	1318	!! conditions on the fluxes are not necessary zalezak (1979)
	1319	!! drange (1995) multi-dimensional forward-in-time and upstream-
	1320	!! in-space based differencing for fluid
	1321	!!----------------------------------------------------------------------
[10419]	1322	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
	1323	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
	1324	REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pu ! ice i-velocity => u*e2
	1325	REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: puc ! ice i-velocity A => ue2*a
	1326	REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pv ! ice j-velocity => v*e1
	1327	REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pvc ! ice j-velocity A => ve1*a
	1328	REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: ptc, pt, pt_low ! before field & upstream guess of after field
	1329	REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfv_low, pfu_low ! upstream flux
	1330	REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfv_ho, pfu_ho ! monotonic flux
[8586]	1331	!
[10180]	1332	INTEGER :: ji, jj, jl ! dummy loop indices
[10419]	1333	REAL(wp) :: zpos, zneg, zbig, zsml, z1_dt, zpos2, zneg2 ! local scalars
	1334	REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign, zcoef ! - -
	1335	REAL(wp), DIMENSION(jpi,jpj ) :: zbup, zbdo
	1336	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zbetup, zbetdo, zti_low, ztj_low, zzt
[8586]	1337	!!----------------------------------------------------------------------
	1338	zbig = 1.e+40_wp
[10419]	1339	zsml = epsi20
	1340
	1341	IF( ll_zeroup2 ) THEN
	1342	DO jl = 1, jpl
	1343	DO jj = 1, jpjm1
	1344	DO ji = 1, fs_jpim1 ! vector opt.
	1345	IF( amaxu(ji,jj,jl) == 0._wp ) pfu_ho(ji,jj,jl) = 0._wp
	1346	IF( amaxv(ji,jj,jl) == 0._wp ) pfv_ho(ji,jj,jl) = 0._wp
	1347	END DO
	1348	END DO
	1349	END DO
	1350	ENDIF
[8586]	1351
[10419]	1352	IF( ll_zeroup4 ) THEN
	1353	DO jl = 1, jpl
	1354	DO jj = 1, jpjm1
	1355	DO ji = 1, fs_jpim1 ! vector opt.
	1356	IF( pfu_low(ji,jj,jl) == 0._wp ) pfu_ho(ji,jj,jl) = 0._wp
	1357	IF( pfv_low(ji,jj,jl) == 0._wp ) pfv_ho(ji,jj,jl) = 0._wp
	1358	END DO
[10180]	1359	END DO
[8586]	1360	END DO
[10419]	1361	ENDIF
	1362
	1363
	1364	IF( ll_zeroup1 ) THEN
	1365	DO jl = 1, jpl
	1366	DO jj = 2, jpjm1
	1367	DO ji = fs_2, fs_jpim1
	1368	IF( ll_gurvan ) THEN
	1369	zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &
	1370	& - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
	1371	ELSE
	1372	zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &
	1373	& - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &
	1374	& + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
	1375	& + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
	1376	& ) * tmask(ji,jj,1)
	1377	ENDIF
	1378	IF( zzt(ji,jj,jl) < 0._wp ) THEN
	1379	pfu_ho(ji,jj,jl) = pfu_low(ji,jj,jl)
	1380	pfv_ho(ji,jj,jl) = pfv_low(ji,jj,jl)
	1381	WRITE(numout,) ' 1 negative high order zzt *',ji,jj,zzt(ji,jj,jl)
	1382	ENDIF
	1383	!! IF( ji==26 .AND. jj==86) THEN
	1384	!! WRITE(numout,*) 'zzt high order',zzt(ji,jj)
	1385	!! WRITE(numout,) 'pfu_ho',(pfu_ho(ji,jj,jl)) r1_e1e2t(ji,jj) * pdt
	1386	!! WRITE(numout,) 'pfv_ho',(pfv_ho(ji,jj,jl)) r1_e1e2t(ji,jj) * pdt
	1387	!! WRITE(numout,) 'pfu_hom1',(pfu_ho(ji-1,jj,jl)) r1_e1e2t(ji,jj) * pdt
	1388	!! WRITE(numout,) 'pfv_hom1',(pfv_ho(ji,jj-1,jl)) r1_e1e2t(ji,jj) * pdt
	1389	!! ENDIF
	1390	IF( ll_gurvan ) THEN
	1391	zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &
	1392	& - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
	1393	ELSE
	1394	zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &
	1395	& - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &
	1396	& + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
	1397	& + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
	1398	& ) * tmask(ji,jj,1)
	1399	ENDIF
	1400	IF( zzt(ji,jj,jl) < 0._wp ) THEN
	1401	pfu_ho(ji-1,jj,jl) = pfu_low(ji-1,jj,jl)
	1402	pfv_ho(ji,jj-1,jl) = pfv_low(ji,jj-1,jl)
	1403	WRITE(numout,) ' 2 negative high order zzt *',ji,jj,zzt(ji,jj,jl)
	1404	ENDIF
	1405	IF( ll_gurvan ) THEN
	1406	zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &
	1407	& - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
	1408	ELSE
	1409	zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &
	1410	& - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &
	1411	& + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
	1412	& + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
	1413	& ) * tmask(ji,jj,1)
	1414	ENDIF
	1415	IF( zzt(ji,jj,jl) < 0._wp ) THEN
	1416	WRITE(numout,) ' 3 negative high order zzt *',ji,jj,zzt(ji,jj,jl)
	1417	ENDIF
	1418	END DO
	1419	END DO
	1420	END DO
	1421	CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. )
	1422	ENDIF
	1423
	1424
	1425	! antidiffusive flux : high order minus low order
	1426	! --------------------------------------------------
	1427	DO jl = 1, jpl
	1428	DO jj = 1, jpjm1
	1429	DO ji = 1, fs_jpim1 ! vector opt.
	1430	pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_low(ji,jj,jl)
	1431	pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_low(ji,jj,jl)
	1432	END DO
	1433	END DO
[8586]	1434	END DO
	1435
[10419]	1436	! extreme case where pfu_ho has to be zero
	1437	! ----------------------------------------
	1438	! pfu_ho
	1439	! * --->
	1440	! \| \| * \| \|
	1441	! \| \| \| * \|
	1442	! \| \| \| \| *
	1443	! t_low : i-1 i i+1 i+2
	1444	IF( ll_prelimiter_zalesak ) THEN
	1445
	1446	DO jl = 1, jpl
	1447	DO jj = 2, jpjm1
	1448	DO ji = fs_2, fs_jpim1
	1449	zti_low(ji,jj,jl)= pt_low(ji+1,jj ,jl)
	1450	ztj_low(ji,jj,jl)= pt_low(ji ,jj+1,jl)
	1451	END DO
	1452	END DO
	1453	END DO
	1454	CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_low, 'T', 1., ztj_low, 'T', 1. )
[8586]	1455
[10419]	1456	!! this does not work ??
	1457	!! DO jj = 2, jpjm1
	1458	!! DO ji = fs_2, fs_jpim1
	1459	!! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji+1,jj ) - pt_low (ji ,jj) ) .AND. &
	1460	!! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji ,jj+1) - pt_low (ji ,jj) ) &
	1461	!! & ) THEN
	1462	!! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_low(ji+1,jj ) - zti_low(ji ,jj) ) .AND. &
	1463	!! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_low(ji,jj+1 ) - ztj_low(ji ,jj) ) &
	1464	!! & ) THEN
	1465	!! pfu_ho(ji,jj) = 0. ; pfv_ho(ji,jj) = 0.
	1466	!! ENDIF
	1467	!! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji ,jj) - pt_low (ji-1,jj ) ) .AND. &
	1468	!! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji ,jj) - pt_low (ji ,jj-1) ) &
	1469	!! & ) THEN
	1470	!! pfu_ho(ji,jj) = 0. ; pfv_ho(ji,jj) = 0.
	1471	!! ENDIF
	1472	!! ENDIF
	1473	!! END DO
	1474	!! END DO
[8586]	1475
[10419]	1476	DO jl = 1, jpl
	1477	DO jj = 2, jpjm1
	1478	DO ji = fs_2, fs_jpim1
	1479	IF ( pfu_ho(ji,jj,jl) * ( pt_low(ji+1,jj,jl) - pt_low(ji,jj,jl) ) <= 0. .AND. &
	1480	& pfv_ho(ji,jj,jl) * ( pt_low(ji,jj+1,jl) - pt_low(ji,jj,jl) ) <= 0. ) THEN
	1481	!
	1482	IF( pfu_ho(ji,jj,jl) * ( zti_low(ji+1,jj,jl) - zti_low(ji,jj,jl) ) <= 0 .AND. &
	1483	& pfv_ho(ji,jj,jl) * ( ztj_low(ji,jj+1,jl) - ztj_low(ji,jj,jl) ) <= 0) pfu_ho(ji,jj,jl)=0. ; pfv_ho(ji,jj,jl)=0.
	1484	!
	1485	IF( pfu_ho(ji,jj,jl) * ( pt_low(ji ,jj,jl) - pt_low(ji-1,jj,jl) ) <= 0 .AND. &
	1486	& pfv_ho(ji,jj,jl) * ( pt_low(ji ,jj,jl) - pt_low(ji,jj-1,jl) ) <= 0) pfu_ho(ji,jj,jl)=0. ; pfv_ho(ji,jj,jl)=0.
	1487	!
	1488	ENDIF
	1489	END DO
	1490	END DO
	1491	END DO
	1492	CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond.
	1493
	1494	ELSEIF( ll_prelimiter_devore ) THEN
	1495	DO jl = 1, jpl
	1496	DO jj = 2, jpjm1
	1497	DO ji = fs_2, fs_jpim1
	1498	zti_low(ji,jj,jl)= pt_low(ji+1,jj ,jl)
	1499	ztj_low(ji,jj,jl)= pt_low(ji ,jj+1,jl)
	1500	END DO
	1501	END DO
	1502	END DO
	1503	CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_low, 'T', 1., ztj_low, 'T', 1. )
	1504
	1505	z1_dt = 1._wp / pdt
	1506	DO jl = 1, jpl
	1507	DO jj = 2, jpjm1
	1508	DO ji = fs_2, fs_jpim1
	1509	zsign = SIGN( 1., pt_low(ji+1,jj,jl) - pt_low(ji,jj,jl) )
	1510	pfu_ho(ji,jj,jl) = zsign * MAX( 0. , MIN( ABS(pfu_ho(ji,jj,jl)) , &
	1511	& zsign * ( pt_low (ji ,jj,jl) - pt_low (ji-1,jj,jl) ) * e1e2t(ji ,jj) * z1_dt , &
	1512	& zsign * ( zti_low(ji+1,jj,jl) - zti_low(ji ,jj,jl) ) * e1e2t(ji+1,jj) * z1_dt ) )
	1513
	1514	zsign = SIGN( 1., pt_low(ji,jj+1,jl) - pt_low(ji,jj,jl) )
	1515	pfv_ho(ji,jj,jl) = zsign * MAX( 0. , MIN( ABS(pfv_ho(ji,jj,jl)) , &
	1516	& zsign * ( pt_low (ji,jj ,jl) - pt_low (ji,jj-1,jl) ) * e1e2t(ji,jj ) * z1_dt , &
	1517	& zsign * ( ztj_low(ji,jj+1,jl) - ztj_low(ji,jj ,jl) ) * e1e2t(ji,jj+1) * z1_dt ) )
	1518	END DO
	1519	END DO
	1520	END DO
	1521	CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond.
	1522
	1523	ENDIF
	1524
	1525
	1526	! Search local extrema
	1527	! --------------------
	1528	! max/min of pt & pt_low with large negative/positive value (-/+zbig) outside ice cover
[8586]	1529	z1_dt = 1._wp / pdt
[10419]	1530	DO jl = 1, jpl
	1531
	1532	DO jj = 1, jpj
	1533	DO ji = 1, jpi
	1534	IF ( pt(ji,jj,jl) <= 0._wp .AND. pt_low(ji,jj,jl) <= 0._wp ) THEN
	1535	zbup(ji,jj) = -zbig
	1536	zbdo(ji,jj) = zbig
	1537	ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_low(ji,jj,jl) > 0._wp ) THEN
	1538	zbup(ji,jj) = pt_low(ji,jj,jl)
	1539	zbdo(ji,jj) = pt_low(ji,jj,jl)
	1540	ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_low(ji,jj,jl) <= 0._wp ) THEN
	1541	zbup(ji,jj) = pt(ji,jj,jl)
	1542	zbdo(ji,jj) = pt(ji,jj,jl)
	1543	ELSE
	1544	zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_low(ji,jj,jl) )
	1545	zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_low(ji,jj,jl) )
	1546	ENDIF
	1547	END DO
	1548	END DO
	1549
	1550	DO jj = 2, jpjm1
	1551	DO ji = fs_2, fs_jpim1 ! vector opt.
[8586]	1552	!
[10419]	1553	IF( .NOT. ll_9points ) THEN
	1554	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
	1555	zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj ), zbdo(ji+1,jj ), zbdo(ji ,jj-1), zbdo(ji ,jj+1) )
	1556	!
	1557	ELSE
	1558	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
	1559	& zbup(ji-1,jj-1), zbup(ji+1,jj+1), zbup(ji+1,jj-1), zbup(ji-1,jj+1) )
	1560	zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj ), zbdo(ji+1,jj ), zbdo(ji ,jj-1), zbdo(ji ,jj+1), &
	1561	& zbdo(ji-1,jj-1), zbdo(ji+1,jj+1), zbdo(ji+1,jj-1), zbdo(ji-1,jj+1) )
	1562	ENDIF
[8586]	1563	!
[10419]	1564	zpos = MAX( 0., pfu_ho(ji-1,jj,jl) ) - MIN( 0., pfu_ho(ji ,jj,jl) ) & ! positive/negative part of the flux
	1565	& + MAX( 0., pfv_ho(ji,jj-1,jl) ) - MIN( 0., pfv_ho(ji,jj ,jl) )
	1566	zneg = MAX( 0., pfu_ho(ji ,jj,jl) ) - MIN( 0., pfu_ho(ji-1,jj,jl) ) &
	1567	& + MAX( 0., pfv_ho(ji,jj ,jl) ) - MIN( 0., pfv_ho(ji,jj-1,jl) )
	1568	!
	1569	IF( ll_HgradU .AND. .NOT.ll_gurvan ) THEN
	1570	zneg2 = ( 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) ) &
	1571	& ) * ( 1. - pamsk )
	1572	zpos2 = ( - 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) ) &
	1573	& ) * ( 1. - pamsk )
	1574	ELSE
	1575	zneg2 = 0. ; zpos2 = 0.
	1576	ENDIF
	1577	!
	1578	! ! up & down beta terms
	1579	IF( (zpos+zpos2) > 0. ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_low(ji,jj,jl) ) / (zpos+zpos2) * e1e2t(ji,jj) * z1_dt
	1580	ELSE ; zbetup(ji,jj,jl) = 0. ! zbig
	1581	ENDIF
	1582	!
	1583	IF( (zneg+zneg2) > 0. ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_low(ji,jj,jl) - zdo ) / (zneg+zneg2) * e1e2t(ji,jj) * z1_dt
	1584	ELSE ; zbetdo(ji,jj,jl) = 0. ! zbig
	1585	ENDIF
	1586	!
	1587	! if all the points are outside ice cover
	1588	IF( zup == -zbig ) zbetup(ji,jj,jl) = 0. ! zbig
	1589	IF( zdo == zbig ) zbetdo(ji,jj,jl) = 0. ! zbig
	1590	!
	1591
	1592	!! IF( ji==26 .AND. jj==86) THEN
	1593	! WRITE(numout,*) '-----------------'
	1594	! WRITE(numout,*) 'zpos',zpos,zpos2
	1595	! WRITE(numout,*) 'zneg',zneg,zneg2
	1596	! WRITE(numout,*) 'puc/pu',ABS(puc(ji,jj))/MAX(epsi20, ABS(pu(ji,jj)))
	1597	! WRITE(numout,*) 'pvc/pv',ABS(pvc(ji,jj))/MAX(epsi20, ABS(pv(ji,jj)))
	1598	! WRITE(numout,*) 'pucm1/pu',ABS(puc(ji-1,jj))/MAX(epsi20, ABS(pu(ji-1,jj)))
	1599	! WRITE(numout,*) 'pvcm1/pv',ABS(pvc(ji,jj-1))/MAX(epsi20, ABS(pv(ji,jj-1)))
	1600	! WRITE(numout,) 'pfu_ho',(pfu_ho(ji,jj)+pfu_low(ji,jj)) r1_e1e2t(ji,jj) * pdt
	1601	! WRITE(numout,) 'pfv_ho',(pfv_ho(ji,jj)+pfv_low(ji,jj)) r1_e1e2t(ji,jj) * pdt
	1602	! WRITE(numout,) 'pfu_hom1',(pfu_ho(ji-1,jj)+pfu_low(ji-1,jj)) r1_e1e2t(ji,jj) * pdt
	1603	! WRITE(numout,) 'pfv_hom1',(pfv_ho(ji,jj-1)+pfv_low(ji,jj-1)) r1_e1e2t(ji,jj) * pdt
	1604	! WRITE(numout,) 'pfu_low',pfu_low(ji,jj) r1_e1e2t(ji,jj) * pdt
	1605	! WRITE(numout,) 'pfv_low',pfv_low(ji,jj) r1_e1e2t(ji,jj) * pdt
	1606	! WRITE(numout,) 'pfu_lowm1',pfu_low(ji-1,jj) r1_e1e2t(ji,jj) * pdt
	1607	! WRITE(numout,) 'pfv_lowm1',pfv_low(ji,jj-1) r1_e1e2t(ji,jj) * pdt
	1608	!
	1609	! WRITE(numout,*) 'pt',pt(ji,jj)
	1610	! WRITE(numout,*) 'ptim1',pt(ji-1,jj)
	1611	! WRITE(numout,*) 'ptjm1',pt(ji,jj-1)
	1612	! WRITE(numout,*) 'pt_low',pt_low(ji,jj)
	1613	! WRITE(numout,*) 'zbetup',zbetup(ji,jj)
	1614	! WRITE(numout,*) 'zbetdo',zbetdo(ji,jj)
	1615	! WRITE(numout,*) 'zup',zup
	1616	! WRITE(numout,*) 'zdo',zdo
	1617	! ENDIF
	1618	!
[10180]	1619	END DO
[8586]	1620	END DO
	1621	END DO
[10170]	1622	CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)
[8586]	1623
[10419]	1624
	1625	! monotonic flux in the y direction
	1626	! ---------------------------------
	1627	DO jl = 1, jpl
	1628	DO jj = 1, jpjm1
[10180]	1629	DO ji = 1, fs_jpim1 ! vector opt.
	1630	zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) )
	1631	zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) )
[10419]	1632	zcu = 0.5 + SIGN( 0.5 , pfu_ho(ji,jj,jl) )
[10180]	1633	!
[10419]	1634	zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu )
	1635
	1636	pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_low(ji,jj,jl)
	1637
	1638	!! IF( ji==26 .AND. jj==86) THEN
	1639	!! WRITE(numout,*) 'coefU',zcoef
	1640	!! WRITE(numout,) 'pfu_ho',(pfu_ho(ji,jj,jl)) r1_e1e2t(ji,jj) * pdt
	1641	!! WRITE(numout,) 'pfu_hom1',(pfu_ho(ji-1,jj,jl)) r1_e1e2t(ji,jj) * pdt
	1642	!! ENDIF
	1643
[10180]	1644	END DO
[8637]	1645	END DO
[10419]	1646
[10180]	1647	DO jj = 1, jpjm1
[10419]	1648	DO ji = 1, fs_jpim1 ! vector opt.
[10180]	1649	zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) )
	1650	zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) )
[10419]	1651	zcv = 0.5 + SIGN( 0.5 , pfv_ho(ji,jj,jl) )
[10180]	1652	!
[10419]	1653	zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv )
	1654
	1655	pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_low(ji,jj,jl)
	1656
	1657	!! IF( ji==26 .AND. jj==86) THEN
	1658	!! WRITE(numout,*) 'coefV',zcoef
	1659	!! WRITE(numout,) 'pfv_ho',(pfv_ho(ji,jj,jl)) r1_e1e2t(ji,jj) * pdt
	1660	!! WRITE(numout,) 'pfv_hom1',(pfv_ho(ji,jj-1,jl)) r1_e1e2t(ji,jj) * pdt
	1661	!! ENDIF
[10180]	1662	END DO
[8586]	1663	END DO
[10419]	1664
	1665	! clem test
	1666	DO jj = 2, jpjm1
	1667	DO ji = 2, fs_jpim1 ! vector opt.
	1668	IF( ll_gurvan ) THEN
	1669	zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &
	1670	& - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
	1671	ELSE
	1672	zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &
	1673	& - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &
	1674	& + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
	1675	& + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
	1676	& ) * tmask(ji,jj,1)
	1677	ENDIF
	1678	IF( zzt(ji,jj,jl) < -epsi20 ) THEN
	1679	WRITE(numout,*) 'T<0 nonosc',zzt(ji,jj,jl)
	1680	ENDIF
	1681	END DO
	1682	END DO
	1683
[8586]	1684	END DO
[10419]	1685
[8586]	1686	!
[10419]	1687	!
[8586]	1688	END SUBROUTINE nonosc_2d
	1689
[10419]	1690	SUBROUTINE limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
	1691	!!---------------------------------------------------------------------
	1692	!! * ROUTINE limiter_x *
	1693	!!
	1694	!! ** Purpose : compute flux limiter
	1695	!!----------------------------------------------------------------------
	1696	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
	1697	REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pu ! ice i-velocity => u*e2
	1698	REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: puc ! ice i-velocity A => ue2*a
	1699	REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pt ! ice tracer
	1700	REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfu_ho ! high order flux
	1701	REAL(wp), DIMENSION (:,:,:), INTENT(in ), OPTIONAL :: pfu_ups ! upstream flux
	1702	!
	1703	REAL(wp) :: Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr
	1704	INTEGER :: ji, jj, jl ! dummy loop indices
	1705	REAL(wp), DIMENSION (jpi,jpj,jpl) :: zslpx ! tracer slopes
	1706	!!----------------------------------------------------------------------
	1707	!
	1708	DO jl = 1, jpl
	1709	DO jj = 2, jpjm1
	1710	DO ji = fs_2, fs_jpim1 ! vector opt.
	1711	zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1)
	1712	END DO
	1713	END DO
	1714	END DO
	1715	CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.) ! lateral boundary cond.
	1716
	1717	DO jl = 1, jpl
	1718	DO jj = 2, jpjm1
	1719	DO ji = fs_2, fs_jpim1 ! vector opt.
	1720	uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj)
	1721
	1722	Rjm = zslpx(ji-1,jj,jl)
	1723	Rj = zslpx(ji ,jj,jl)
	1724	Rjp = zslpx(ji+1,jj,jl)
	1725
	1726	IF( PRESENT(pfu_ups) ) THEN
	1727
	1728	IF( pu(ji,jj) > 0. ) THEN ; Rr = Rjm
	1729	ELSE ; Rr = Rjp
	1730	ENDIF
	1731
	1732	zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)
	1733	IF( Rj > 0. ) THEN
	1734	zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pu(ji,jj)), &
	1735	& MIN( 2. * Rr * 0.5 * ABS(pu(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
	1736	ELSE
	1737	zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)), &
	1738	& MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
	1739	ENDIF
	1740	pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter
	1741
	1742	ELSE
	1743	IF( Rj /= 0. ) THEN
	1744	IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj
	1745	ELSE ; Cr = Rjp / Rj
	1746	ENDIF
	1747	ELSE
	1748	Cr = 0.
	1749	!IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm * 1.e20
	1750	!ELSE ; Cr = Rjp * 1.e20
	1751	!ENDIF
	1752	ENDIF
	1753
	1754	! -- superbee --
	1755	zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
	1756	! -- van albada 2 --
	1757	!!zpsi = 2.Cr / (CrCr+1.)
	1758
	1759	! -- sweby (with beta=1) --
	1760	!!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
	1761	! -- van Leer --
	1762	!!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
	1763	! -- ospre --
	1764	!!zpsi = 1.5 * ( CrCr + Cr ) / ( CrCr + Cr + 1. )
	1765	! -- koren --
	1766	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( (1.+2Cr)/3., 2. ) ) )
	1767	! -- charm --
	1768	!IF( Cr > 0. ) THEN ; zpsi = Cr * (3.Cr + 1.) / ( (Cr + 1.) (Cr + 1.) )
	1769	!ELSE ; zpsi = 0.
	1770	!ENDIF
	1771	! -- van albada 1 --
	1772	!!zpsi = (CrCr + Cr) / (CrCr +1)
	1773	! -- smart --
	1774	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, 4. ) ) )
	1775	! -- umist --
	1776	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
	1777
	1778	! high order flux corrected by the limiter
	1779	pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFLzpsi ) Rj * 0.5
	1780
	1781	ENDIF
	1782	END DO
	1783	END DO
	1784	END DO
	1785	CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.) ! lateral boundary cond.
	1786	!
	1787	END SUBROUTINE limiter_x
	1788
	1789	SUBROUTINE limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
	1790	!!---------------------------------------------------------------------
	1791	!! * ROUTINE limiter_y *
	1792	!!
	1793	!! ** Purpose : compute flux limiter
	1794	!!----------------------------------------------------------------------
	1795	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
	1796	REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pv ! ice i-velocity => u*e2
	1797	REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pvc ! ice i-velocity A => ue2*a
	1798	REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pt ! ice tracer
	1799	REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfv_ho ! high order flux
	1800	REAL(wp), DIMENSION (:,:,:), INTENT(in ), OPTIONAL :: pfv_ups ! upstream flux
	1801	!
	1802	REAL(wp) :: Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr
	1803	INTEGER :: ji, jj, jl ! dummy loop indices
	1804	REAL(wp), DIMENSION (jpi,jpj,jpl) :: zslpy ! tracer slopes
	1805	!!----------------------------------------------------------------------
	1806	!
	1807	DO jl = 1, jpl
	1808	DO jj = 2, jpjm1
	1809	DO ji = fs_2, fs_jpim1 ! vector opt.
	1810	zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1)
	1811	END DO
	1812	END DO
	1813	END DO
	1814	CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.) ! lateral boundary cond.
	1815
	1816	DO jl = 1, jpl
	1817	DO jj = 2, jpjm1
	1818	DO ji = fs_2, fs_jpim1 ! vector opt.
	1819	vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj)
	1820
	1821	Rjm = zslpy(ji,jj-1,jl)
	1822	Rj = zslpy(ji,jj ,jl)
	1823	Rjp = zslpy(ji,jj+1,jl)
	1824
	1825	IF( PRESENT(pfv_ups) ) THEN
	1826
	1827	IF( pv(ji,jj) > 0. ) THEN ; Rr = Rjm
	1828	ELSE ; Rr = Rjp
	1829	ENDIF
	1830
	1831	zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)
	1832	IF( Rj > 0. ) THEN
	1833	zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pv(ji,jj)), &
	1834	& MIN( 2. * Rr * 0.5 * ABS(pv(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
	1835	ELSE
	1836	zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)), &
	1837	& MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
	1838	ENDIF
	1839	pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter
	1840
	1841	ELSE
	1842
	1843	IF( Rj /= 0. ) THEN
	1844	IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj
	1845	ELSE ; Cr = Rjp / Rj
	1846	ENDIF
	1847	ELSE
	1848	Cr = 0.
	1849	!IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm * 1.e20
	1850	!ELSE ; Cr = Rjp * 1.e20
	1851	!ENDIF
	1852	ENDIF
	1853
	1854	! -- superbee --
	1855	zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
	1856	! -- van albada 2 --
	1857	!!zpsi = 2.Cr / (CrCr+1.)
	1858
	1859	! -- sweby (with beta=1) --
	1860	!!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
	1861	! -- van Leer --
	1862	!!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
	1863	! -- ospre --
	1864	!!zpsi = 1.5 * ( CrCr + Cr ) / ( CrCr + Cr + 1. )
	1865	! -- koren --
	1866	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( (1.+2Cr)/3., 2. ) ) )
	1867	! -- charm --
	1868	!IF( Cr > 0. ) THEN ; zpsi = Cr * (3.Cr + 1.) / ( (Cr + 1.) (Cr + 1.) )
	1869	!ELSE ; zpsi = 0.
	1870	!ENDIF
	1871	! -- van albada 1 --
	1872	!!zpsi = (CrCr + Cr) / (CrCr +1)
	1873	! -- smart --
	1874	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, 4. ) ) )
	1875	! -- umist --
	1876	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
	1877
	1878	! high order flux corrected by the limiter
	1879	pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFLzpsi ) Rj * 0.5
	1880
	1881	ENDIF
	1882	END DO
	1883	END DO
	1884	END DO
	1885	CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.) ! lateral boundary cond.
	1886	!
	1887	END SUBROUTINE limiter_y
	1888
[8586]	1889	#else
	1890	!!----------------------------------------------------------------------
[9570]	1891	!! Default option Dummy module NO SI3 sea-ice model
[8586]	1892	!!----------------------------------------------------------------------
	1893	#endif
	1894
	1895	!!======================================================================
	1896	END MODULE icedyn_adv_umx

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

Download in other formats: