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

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

icedyn_adv_umx.F90 in NEMO/trunk/src/ICE – NEMO

Context Navigation

source: NEMO/trunk/src/ICE/icedyn_adv_umx.F90 @ 10579

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

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

Download in other formats: