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/releases/release-4.0/src/ICE – NEMO

Context Navigation

source: NEMO/releases/release-4.0/src/ICE/icedyn_adv_umx.F90 @ 10785

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

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

Download in other formats: