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

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

icedyn_adv_umx.F90 in NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icedyn_adv_umx.F90 @ 12379

Last change on this file since 12379 was 12379, checked in by dancopsey, 4 years ago
Add meltpond lid thickness as a new prognostic.
File size: 84.6 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	!!----------------------------------------------------------------------
[11081]	13	!! ice_dyn_adv_umx : update the tracer fields
[8586]	14	!! ultimate_x(_y) : compute a tracer value at velocity points using ULTIMATE scheme at various orders
[11081]	15	!! macho : compute the fluxes
	16	!! nonosc_ice : limit the fluxes using 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
[11081]	34	!
	35	INTEGER, PARAMETER :: np_advS = 1 ! advection for S and T: dVS/dt = -div( uVS ) => np_advS = 1
	36	! or dVS/dt = -div( uA * uHS / u ) => np_advS = 2
	37	! or dVS/dt = -div( uV * uS / u ) => np_advS = 3
	38	INTEGER, PARAMETER :: np_limiter = 1 ! limiter: 1 = nonosc
	39	! 2 = superbee
	40	! 3 = h3
	41	LOGICAL :: ll_upsxy = .TRUE. ! alternate directions for upstream
	42	LOGICAL :: ll_hoxy = .TRUE. ! alternate directions for high order
	43	LOGICAL :: ll_neg = .TRUE. ! if T interpolated at u/v points is negative or v_i < 1.e-6
	44	! then interpolate T at u/v points using the upstream scheme
	45	LOGICAL :: ll_prelim = .FALSE. ! prelimiter from: Zalesak(1979) eq. 14 => not well defined in 2D
	46	!
	47	REAL(wp) :: z1_6 = 1._wp / 6._wp ! =1/6
	48	REAL(wp) :: z1_120 = 1._wp / 120._wp ! =1/120
	49	!
	50	INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: imsk_small, jmsk_small
	51	!
[8586]	52	!! * Substitutions
	53	# include "vectopt_loop_substitute.h90"
	54	!!----------------------------------------------------------------------
[9598]	55	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[10069]	56	!! $Id$
[10413]	57	!! Software governed by the CeCILL licence (./LICENSE)
[8586]	58	!!----------------------------------------------------------------------
	59	CONTAINS
	60
[11081]	61	SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip, &
[12379]	62	& pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i )
[8586]	63	!!----------------------------------------------------------------------
	64	!! * ROUTINE ice_dyn_adv_umx *
	65	!!
	66	!! ** Purpose : Compute the now trend due to total advection of
	67	!! tracers and add it to the general trend of tracer equations
	68	!! using an "Ultimate-Macho" scheme
	69	!!
	70	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
	71	!!----------------------------------------------------------------------
[10413]	72	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
[8586]	73	INTEGER , INTENT(in ) :: kt ! time step
	74	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity
	75	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pv_ice ! ice j-velocity
[11081]	76	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_i ! ice thickness
	77	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_s ! snw thickness
	78	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_ip ! ice pond thickness
[8586]	79	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area
	80	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume
	81	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_s ! snw volume
	82	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psv_i ! salt content
	83	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content
	84	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration
	85	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
	86	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume
[12379]	87	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: plh_ip ! melt pond lid thickness
[8586]	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
[11081]	94	REAL(wp) :: zdt, zvi_cen
	95	REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication
	96	REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box
[10439]	97	REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2
[11081]	98	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zu_cat, zv_cat
	99	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zua_ho, zva_ho, zua_ups, zva_ups
	100	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_ai , z1_aip, zhvar
	101	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max
	102	!
	103	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs
[8586]	104	!!----------------------------------------------------------------------
	105	!
	106	IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
	107	!
[11081]	108	! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- !
	109	DO jl = 1, jpl
	110	DO jj = 2, jpjm1
	111	DO ji = fs_2, fs_jpim1
	112	zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), &
	113	& ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), &
	114	& ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), &
	115	& ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) )
	116	zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), &
	117	& ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), &
	118	& ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), &
	119	& ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) )
	120	zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), &
	121	& ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), &
	122	& ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), &
	123	& ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) )
	124	END DO
	125	END DO
	126	END DO
	127	CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. )
	128	!
	129	!
	130	! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- !
	131	! Note: the advection split is applied at the next time-step in order to avoid blocking global comm.
	132	! this should not affect too much the stability
[10425]	133	zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
	134	zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
	135
	136	! non-blocking global communication send zcflnow and receive zcflprv
	137	CALL mpp_delay_max( 'icedyn_adv_umx', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 )
[8586]	138
[10425]	139	IF( zcflprv(1) > .5 ) THEN ; icycle = 2
	140	ELSE ; icycle = 1
[8586]	141	ENDIF
[10413]	142	zdt = rdt_ice / REAL(icycle)
[8586]	143
	144	! --- transport --- !
	145	zudy(:,:) = pu_ice(:,:) * e2u(:,:)
	146	zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
[11081]	147	!
	148	! setup transport for each ice cat
	149	DO jl = 1, jpl
	150	zu_cat(:,:,jl) = zudy(:,:)
	151	zv_cat(:,:,jl) = zvdx(:,:)
	152	END DO
	153	!
[8586]	154	! --- define velocity for advection: u*grad(H) --- !
	155	DO jj = 2, jpjm1
	156	DO ji = fs_2, fs_jpim1
	157	IF ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN ; zcu_box(ji,jj) = 0._wp
	158	ELSEIF( pu_ice(ji,jj) > 0._wp ) THEN ; zcu_box(ji,jj) = pu_ice(ji-1,jj)
	159	ELSE ; zcu_box(ji,jj) = pu_ice(ji ,jj)
	160	ENDIF
	161
	162	IF ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN ; zcv_box(ji,jj) = 0._wp
	163	ELSEIF( pv_ice(ji,jj) > 0._wp ) THEN ; zcv_box(ji,jj) = pv_ice(ji,jj-1)
	164	ELSE ; zcv_box(ji,jj) = pv_ice(ji,jj )
	165	ENDIF
	166	END DO
	167	END DO
	168
	169	!---------------!
	170	!== advection ==!
	171	!---------------!
[10413]	172	DO jt = 1, icycle
	173
[10439]	174	! record at_i before advection (for open water)
	175	zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
[10413]	176
[10439]	177	! inverse of A and Ap
[10425]	178	WHERE( pa_i(:,:,:) >= epsi20 ) ; z1_ai(:,:,:) = 1._wp / pa_i(:,:,:)
	179	ELSEWHERE ; z1_ai(:,:,:) = 0.
	180	END WHERE
	181	WHERE( pa_ip(:,:,:) >= epsi20 ) ; z1_aip(:,:,:) = 1._wp / pa_ip(:,:,:)
	182	ELSEWHERE ; z1_aip(:,:,:) = 0.
	183	END WHERE
	184	!
[11081]	185	! setup a mask where advection will be upstream
	186	IF( ll_neg ) THEN
	187	IF( .NOT. ALLOCATED(imsk_small) ) ALLOCATE( imsk_small(jpi,jpj,jpl) )
	188	IF( .NOT. ALLOCATED(jmsk_small) ) ALLOCATE( jmsk_small(jpi,jpj,jpl) )
	189	DO jl = 1, jpl
	190	DO jj = 1, jpjm1
	191	DO ji = 1, jpim1
	192	zvi_cen = 0.5_wp * ( pv_i(ji+1,jj,jl) + pv_i(ji,jj,jl) )
	193	IF( zvi_cen < epsi06) THEN ; imsk_small(ji,jj,jl) = 0
	194	ELSE ; imsk_small(ji,jj,jl) = 1 ; ENDIF
	195	zvi_cen = 0.5_wp * ( pv_i(ji,jj+1,jl) + pv_i(ji,jj,jl) )
	196	IF( zvi_cen < epsi06) THEN ; jmsk_small(ji,jj,jl) = 0
	197	ELSE ; jmsk_small(ji,jj,jl) = 1 ; ENDIF
	198	END DO
	199	END DO
	200	END DO
	201	ENDIF
[10425]	202	!
[11081]	203	! ----------------------- !
	204	! ==> start advection <== !
	205	! ----------------------- !
[10425]	206	!
[11081]	207	!== Ice area ==!
	208	zamsk = 1._wp
	209	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zu_cat , zv_cat , zcu_box, zcv_box, &
	210	& pa_i, pa_i, zua_ups, zva_ups, zua_ho , zva_ho )
[10425]	211	!
[11081]	212	! ! --------------------------------- !
	213	IF( np_advS == 1 ) THEN ! -- advection form: -div( uVS ) -- !
	214	! ! --------------------------------- !
	215	zamsk = 0._wp
	216	!== Ice volume ==!
	217	zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
	218	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
	219	& zhvar, pv_i, zua_ups, zva_ups )
	220	!== Snw volume ==!
	221	zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
	222	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
	223	& zhvar, pv_s, zua_ups, zva_ups )
	224	!
	225	zamsk = 1._wp
	226	!== Salt content ==!
	227	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
	228	& psv_i, psv_i )
	229	!== Ice heat content ==!
	230	DO jk = 1, nlay_i
	231	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
	232	& pe_i(:,:,jk,:), pe_i(:,:,jk,:) )
[10785]	233	END DO
[11081]	234	!== Snw heat content ==!
	235	DO jk = 1, nlay_s
	236	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
	237	& pe_s(:,:,jk,:), pe_s(:,:,jk,:) )
	238	END DO
	239	!
	240	! ! ------------------------------------------ !
	241	ELSEIF( np_advS == 2 ) THEN ! -- advection form: -div( uA * uHS / u ) -- !
	242	! ! ------------------------------------------ !
[10785]	243	zamsk = 0._wp
[11081]	244	!== Ice volume ==!
	245	zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
	246	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
	247	& zhvar, pv_i, zua_ups, zva_ups )
	248	!== Snw volume ==!
	249	zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
	250	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
	251	& zhvar, pv_s, zua_ups, zva_ups )
	252	!== Salt content ==!
	253	zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:)
	254	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
	255	& zhvar, psv_i, zua_ups, zva_ups )
	256	!== Ice heat content ==!
	257	DO jk = 1, nlay_i
	258	zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:)
	259	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, &
	260	& zhvar, pe_i(:,:,jk,:), zua_ups, zva_ups )
	261	END DO
	262	!== Snw heat content ==!
	263	DO jk = 1, nlay_s
	264	zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:)
	265	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, &
	266	& zhvar, pe_s(:,:,jk,:), zua_ups, zva_ups )
	267	END DO
	268	!
	269	! ! ----------------------------------------- !
	270	ELSEIF( np_advS == 3 ) THEN ! -- advection form: -div( uV * uS / u ) -- !
	271	! ! ----------------------------------------- !
	272	zamsk = 0._wp
	273	!
	274	ALLOCATE( zuv_ho (jpi,jpj,jpl), zvv_ho (jpi,jpj,jpl), &
	275	& zuv_ups(jpi,jpj,jpl), zvv_ups(jpi,jpj,jpl), z1_vi(jpi,jpj,jpl), z1_vs(jpi,jpj,jpl) )
	276	!
	277	! inverse of Vi
	278	WHERE( pv_i(:,:,:) >= epsi20 ) ; z1_vi(:,:,:) = 1._wp / pv_i(:,:,:)
	279	ELSEWHERE ; z1_vi(:,:,:) = 0.
	280	END WHERE
	281	! inverse of Vs
	282	WHERE( pv_s(:,:,:) >= epsi20 ) ; z1_vs(:,:,:) = 1._wp / pv_s(:,:,:)
	283	ELSEWHERE ; z1_vs(:,:,:) = 0.
	284	END WHERE
	285	!
	286	! It is important to first calculate the ice fields and then the snow fields (because we use the same arrays)
	287	!
	288	!== Ice volume ==!
	289	zuv_ups = zua_ups
	290	zvv_ups = zva_ups
	291	zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
	292	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
	293	& zhvar, pv_i, zuv_ups, zvv_ups, zuv_ho , zvv_ho )
	294	!== Salt content ==!
	295	zhvar(:,:,:) = psv_i(:,:,:) * z1_vi(:,:,:)
	296	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zuv_ho , zvv_ho , zcu_box, zcv_box, &
	297	& zhvar, psv_i, zuv_ups, zvv_ups )
	298	!== Ice heat content ==!
	299	DO jk = 1, nlay_i
	300	zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_vi(:,:,:)
	301	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, &
	302	& zhvar, pe_i(:,:,jk,:), zuv_ups, zvv_ups )
	303	END DO
	304	!== Snow volume ==!
	305	zuv_ups = zua_ups
	306	zvv_ups = zva_ups
	307	zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
	308	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, &
	309	& zhvar, pv_s, zuv_ups, zvv_ups, zuv_ho , zvv_ho )
	310	!== Snw heat content ==!
	311	DO jk = 1, nlay_s
	312	zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_vs(:,:,:)
	313	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, &
	314	& zhvar, pe_s(:,:,jk,:), zuv_ups, zvv_ups )
	315	END DO
	316	!
	317	DEALLOCATE( zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs )
	318	!
[10785]	319	ENDIF
	320	!
[11081]	321	!== Ice age ==!
	322	IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN
[10425]	323	zamsk = 1._wp
[11081]	324	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
	325	& poa_i, poa_i )
	326	ENDIF
	327	!
	328	!== melt ponds ==!
	329	IF ( ln_pnd_H12 ) THEN
	330	! fraction
	331	zamsk = 1._wp
	332	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat , zv_cat , zcu_box, zcv_box, &
	333	& pa_ip, pa_ip, zua_ups, zva_ups, zua_ho , zva_ho )
	334	! volume
[10425]	335	zamsk = 0._wp
[10475]	336	zhvar(:,:,:) = pv_ip(:,:,:) * z1_aip(:,:,:)
[11081]	337	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
	338	& zhvar, pv_ip, zua_ups, zva_ups )
[12379]	339	! lid thickness
	340	zamsk = 0._wp
	341	zhvar(:,:,:) = plh_ip(:,:,:) * z1_aip(:,:,:)
	342	CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
	343	& zhvar, plh_ip, zua_ups, zva_ups )
	344
[10425]	345	ENDIF
[10418]	346	!
[11081]	347	!== Open water area ==!
[10439]	348	zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
	349	DO jj = 2, jpjm1
	350	DO ji = fs_2, fs_jpim1
[11081]	351	pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &
[10439]	352	& - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
[10413]	353	END DO
[10439]	354	END DO
[11081]	355	CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1. )
[10418]	356	!
[11081]	357	!
	358	! --- Ensure non-negative fields and in-bound thicknesses --- !
	359	! Remove negative values (conservation is ensured)
	360	! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
[12379]	361	CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i )
[11081]	362	!
	363	! Make sure ice thickness is not too big
	364	! (because ice thickness can be too large where ice concentration is very small)
	365	CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
	366
[8586]	367	END DO
	368	!
	369	END SUBROUTINE ice_dyn_adv_umx
[9929]	370
[8586]	371
[11081]	372	SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, &
	373	& pt, ptc, pua_ups, pva_ups, pua_ho, pva_ho )
[8586]	374	!!----------------------------------------------------------------------
	375	!! * ROUTINE adv_umx *
	376	!!
	377	!! ** Purpose : Compute the now trend due to total advection of
[10446]	378	!! tracers and add it to the general trend of tracer equations
[8586]	379	!!
[11081]	380	!! ** Method : - calculate upstream fluxes and upstream solution for tracers V/A(=H) etc
[10446]	381	!! - calculate tracer H at u and v points (Ultimate)
[11081]	382	!! - calculate the high order fluxes using alterning directions (Macho)
[10519]	383	!! - apply a limiter on the fluxes (nonosc_ice)
[11081]	384	!! - convert this tracer flux to a "volume" flux (uH -> uV)
	385	!! - apply a limiter a second time on the volumes fluxes (nonosc_ice)
	386	!! - calculate the high order solution for V
[8586]	387	!!
[11081]	388	!! ** Action : solve 3 equations => a) dA/dt = -div(uA)
	389	!! b) dV/dt = -div(uV) using dH/dt = -u.grad(H)
	390	!! c) dVS/dt = -div(uVS) using either dHS/dt = -u.grad(HS) or dS/dt = -u.grad(S)
[10446]	391	!!
[11081]	392	!! in eq. b), - fluxes uH are evaluated (with UMx) and limited with nonosc_ice. This step is necessary to get a good H.
	393	!! - then we convert this flux to a "volume" flux this way => uH * uA / u
	394	!! where uA is the flux from eq. a)
	395	!! this "volume" flux is also limited with nonosc_ice (otherwise overshoots can occur)
	396	!! - at last we estimate dV/dt = -div(uH * uA / u)
	397	!!
	398	!! in eq. c), one can solve the equation for S (ln_advS=T), then dVS/dt = -div(uV * uS / u)
	399	!! or for HS (ln_advS=F), then dVS/dt = -div(uA * uHS / u)
	400	!!
	401	!! ** Note : - this method can lead to tiny negative V (-1.e-20) => set it to 0 while conserving mass etc.
	402	!! - At the ice edge, Ultimate scheme can lead to:
	403	!! 1) negative interpolated tracers at u-v points
	404	!! 2) non-zero interpolated tracers at u-v points eventhough there is no ice and velocity is outward
	405	!! Solution for 1): apply an upstream scheme when it occurs. A better solution would be to degrade the order of
	406	!! the scheme automatically by applying a mask of the ice cover inside Ultimate (not done).
	407	!! Solution for 2): we set it to 0 in this case
[10446]	408	!! - Eventhough 1D tests give very good results (typically the one from Schar & Smolarkiewiecz), the 2D is less good.
	409	!! Large values of H can appear for very small ice concentration, and when it does it messes the things up since we
[11081]	410	!! work on H (and not V). It is partly related to the multi-category approach
[10446]	411	!! Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice
[11081]	412	!! concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter
	413	!! since sv_i and e_i are still good.
[8586]	414	!!----------------------------------------------------------------------
[11081]	415	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
	416	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
	417	INTEGER , INTENT(in ) :: jt ! number of sub-iteration
	418	INTEGER , INTENT(in ) :: kt ! number of iteration
	419	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
	420	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu , pv ! 2 ice velocity components => u*e2
	421	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: puc , pvc ! 2 ice velocity components => ue2 or ua*e2u
	422	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pubox, pvbox ! upstream velocity
	423	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pt ! tracer field
	424	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: ptc ! tracer content field
	425	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout), OPTIONAL :: pua_ups, pva_ups ! upstream u*a fluxes
	426	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out), OPTIONAL :: pua_ho, pva_ho ! high order u*a fluxes
[8586]	427	!
[10425]	428	INTEGER :: ji, jj, jl ! dummy loop indices
[8586]	429	REAL(wp) :: ztra ! local scalar
[10446]	430	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zfu_ho , zfv_ho , zpt
[10439]	431	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zfu_ups, zfv_ups, zt_ups
[8586]	432	!!----------------------------------------------------------------------
	433	!
[10446]	434	! Upstream (_ups) fluxes
	435	! -----------------------
	436	CALL upstream( pamsk, jt, kt, pdt, pt, pu, pv, zt_ups, zfu_ups, zfv_ups )
	437
	438	! High order (_ho) fluxes
	439	! -----------------------
	440	SELECT CASE( kn_umx )
	441	!
	442	CASE ( 20 ) !== centered second order ==!
	443	!
[10475]	444	CALL cen2( pamsk, jt, kt, pdt, pt, pu, pv, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho )
[10446]	445	!
	446	CASE ( 1:5 ) !== 1st to 5th order ULTIMATE-MACHO scheme ==!
	447	!
[10475]	448	CALL macho( pamsk, kn_umx, jt, kt, pdt, pt, pu, pv, pubox, pvbox, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho )
[10446]	449	!
	450	END SELECT
[10439]	451	!
[10446]	452	! --ho --ho
	453	! new fluxes = uH u*a / u
	454	! ----------------------------
[10475]	455	IF( pamsk == 0._wp ) THEN
[10446]	456	DO jl = 1, jpl
	457	DO jj = 1, jpjm1
	458	DO ji = 1, fs_jpim1
[11081]	459	IF( ABS( pu(ji,jj) ) > epsi10 ) THEN
	460	zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc (ji,jj,jl) / pu(ji,jj)
	461	zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj)
[10446]	462	ELSE
	463	zfu_ho (ji,jj,jl) = 0._wp
	464	zfu_ups(ji,jj,jl) = 0._wp
	465	ENDIF
	466	!
[11081]	467	IF( ABS( pv(ji,jj) ) > epsi10 ) THEN
	468	zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc (ji,jj,jl) / pv(ji,jj)
	469	zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj)
[10446]	470	ELSE
	471	zfv_ho (ji,jj,jl) = 0._wp
	472	zfv_ups(ji,jj,jl) = 0._wp
	473	ENDIF
	474	END DO
	475	END DO
	476	END DO
[11081]	477
	478	! the new "volume" fluxes must also be "flux corrected"
	479	! thus we calculate the upstream solution and apply a limiter again
	480	DO jl = 1, jpl
	481	DO jj = 2, jpjm1
	482	DO ji = fs_2, fs_jpim1
	483	ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) )
	484	!
	485	zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)
	486	END DO
	487	END DO
	488	END DO
	489	CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. )
	490	!
	491	IF ( np_limiter == 1 ) THEN
	492	CALL nonosc_ice( 1._wp, pdt, pu, pv, ptc, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho )
	493	ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN
	494	CALL limiter_x( pdt, pu, ptc, zfu_ups, zfu_ho )
	495	CALL limiter_y( pdt, pv, ptc, zfv_ups, zfv_ho )
	496	ENDIF
	497	!
[10446]	498	ENDIF
[11081]	499	! --ho --ups
	500	! in case of advection of A: output ua and ua
	501	! -----------------------------------------------
[10446]	502	IF( PRESENT( pua_ho ) ) THEN
	503	DO jl = 1, jpl
	504	DO jj = 1, jpjm1
	505	DO ji = 1, fs_jpim1
[11081]	506	pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl)
	507	pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl)
	508	END DO
[10446]	509	END DO
	510	END DO
	511	ENDIF
	512	!
	513	! final trend with corrected fluxes
	514	! ---------------------------------
	515	DO jl = 1, jpl
	516	DO jj = 2, jpjm1
	517	DO ji = fs_2, fs_jpim1
	518	ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) )
	519	!
	520	ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)
	521	END DO
	522	END DO
	523	END DO
	524	CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T', 1. )
	525	!
	526	END SUBROUTINE adv_umx
	527
	528
	529	SUBROUTINE upstream( pamsk, jt, kt, pdt, pt, pu, pv, pt_ups, pfu_ups, pfv_ups )
	530	!!---------------------------------------------------------------------
	531	!! * ROUTINE upstream *
	532	!!
	533	!! ** Purpose : compute the upstream fluxes and upstream guess of tracer
	534	!!----------------------------------------------------------------------
	535	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
	536	INTEGER , INTENT(in ) :: jt ! number of sub-iteration
	537	INTEGER , INTENT(in ) :: kt ! number of iteration
	538	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
	539	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields
	540	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu, pv ! 2 ice velocity components
	541	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_ups ! upstream guess of tracer
	542	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ups, pfv_ups ! upstream fluxes
	543	!
	544	INTEGER :: ji, jj, jl ! dummy loop indices
	545	REAL(wp) :: ztra ! local scalar
	546	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zpt
	547	!!----------------------------------------------------------------------
	548
[10439]	549	IF( .NOT. ll_upsxy ) THEN ! no alternate directions !
[10446]	550	!
[10425]	551	DO jl = 1, jpl
	552	DO jj = 1, jpjm1
	553	DO ji = 1, fs_jpim1
[10446]	554	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)
	555	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]	556	END DO
[10413]	557	END DO
	558	END DO
[10446]	559	!
[10439]	560	ELSE ! alternate directions !
[10413]	561	!
	562	IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==!
[10439]	563	!
	564	DO jl = 1, jpl !-- flux in x-direction
[10425]	565	DO jj = 1, jpjm1
	566	DO ji = 1, fs_jpim1
[10446]	567	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]	568	END DO
[10413]	569	END DO
	570	END DO
[10439]	571	!
	572	DO jl = 1, jpl !-- first guess of tracer from u-flux
[10425]	573	DO jj = 2, jpjm1
[10439]	574	DO ji = fs_2, fs_jpim1
[10446]	575	ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) ) &
	576	& + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk)
	577	!
	578	zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
[10425]	579	END DO
[10413]	580	END DO
	581	END DO
[10425]	582	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
[8586]	583	!
[10439]	584	DO jl = 1, jpl !-- flux in y-direction
[10425]	585	DO jj = 1, jpjm1
	586	DO ji = 1, fs_jpim1
[10446]	587	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]	588	END DO
[10413]	589	END DO
	590	END DO
[10439]	591	!
[10413]	592	ELSE !== even ice time step: adv_y then adv_x ==!
[10439]	593	!
	594	DO jl = 1, jpl !-- flux in y-direction
[10425]	595	DO jj = 1, jpjm1
	596	DO ji = 1, fs_jpim1
[10446]	597	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]	598	END DO
[10413]	599	END DO
	600	END DO
[10439]	601	!
	602	DO jl = 1, jpl !-- first guess of tracer from v-flux
[10425]	603	DO jj = 2, jpjm1
[10439]	604	DO ji = fs_2, fs_jpim1
[10446]	605	ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) ) &
	606	& + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk)
	607	!
	608	zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
[10425]	609	END DO
	610	END DO
[10413]	611	END DO
[10425]	612	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
[10413]	613	!
[10439]	614	DO jl = 1, jpl !-- flux in x-direction
[10425]	615	DO jj = 1, jpjm1
	616	DO ji = 1, fs_jpim1
[10446]	617	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]	618	END DO
[10413]	619	END DO
	620	END DO
	621	!
	622	ENDIF
	623
	624	ENDIF
[10439]	625	!
	626	DO jl = 1, jpl !-- after tracer with upstream scheme
[10425]	627	DO jj = 2, jpjm1
	628	DO ji = fs_2, fs_jpim1
[10446]	629	ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj ,jl) &
	630	& + pfv_ups(ji,jj,jl) - pfv_ups(ji ,jj-1,jl) ) &
	631	& + ( pu (ji,jj ) - pu (ji-1,jj ) &
	632	& + pv (ji,jj ) - pv (ji ,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk)
	633	!
[10475]	634	pt_ups(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
[10425]	635	END DO
[10413]	636	END DO
[8586]	637	END DO
[10446]	638	CALL lbc_lnk( 'icedyn_adv_umx', pt_ups, 'T', 1. )
[10413]	639
[10446]	640	END SUBROUTINE upstream
[8586]	641
[10446]	642
[10475]	643	SUBROUTINE cen2( pamsk, jt, kt, pdt, pt, pu, pv, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
[8586]	644	!!---------------------------------------------------------------------
[10446]	645	!! * ROUTINE cen2 *
[8586]	646	!!
[10446]	647	!! ** Purpose : compute the high order fluxes using a centered
	648	!! second order scheme
[8586]	649	!!----------------------------------------------------------------------
[10439]	650	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
	651	INTEGER , INTENT(in ) :: jt ! number of sub-iteration
	652	INTEGER , INTENT(in ) :: kt ! number of iteration
	653	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
	654	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields
	655	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu, pv ! 2 ice velocity components
[10446]	656	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt_ups ! upstream guess of tracer
	657	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pfu_ups, pfv_ups ! upstream fluxes
[10425]	658	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes
[8586]	659	!
[10425]	660	INTEGER :: ji, jj, jl ! dummy loop indices
[10446]	661	REAL(wp) :: ztra ! local scalar
	662	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zpt
[8586]	663	!!----------------------------------------------------------------------
	664	!
[10439]	665	IF( .NOT.ll_hoxy ) THEN ! no alternate directions !
[8586]	666	!
[10425]	667	DO jl = 1, jpl
	668	DO jj = 1, jpjm1
	669	DO ji = 1, fs_jpim1
[10475]	670	pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj ,jl) )
	671	pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji ,jj+1,jl) )
[10425]	672	END DO
[8586]	673	END DO
	674	END DO
[10475]	675	!
[11081]	676	IF ( np_limiter == 1 ) THEN
[10519]	677	CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
[11081]	678	ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN
[10446]	679	CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
	680	CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
[10413]	681	ENDIF
[8586]	682	!
[10439]	683	ELSE ! alternate directions !
[8586]	684	!
[10413]	685	IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==!
	686	!
[10439]	687	DO jl = 1, jpl !-- flux in x-direction
[10425]	688	DO jj = 1, jpjm1
	689	DO ji = 1, fs_jpim1
[10475]	690	pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
[10425]	691	END DO
[10413]	692	END DO
	693	END DO
[11081]	694	IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
[10413]	695
[10439]	696	DO jl = 1, jpl !-- first guess of tracer from u-flux
[10425]	697	DO jj = 2, jpjm1
[10439]	698	DO ji = fs_2, fs_jpim1
[10446]	699	ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) &
	700	& + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk)
	701	!
	702	zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
[10425]	703	END DO
[10413]	704	END DO
	705	END DO
[10446]	706	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
[10413]	707
[10439]	708	DO jl = 1, jpl !-- flux in y-direction
[10425]	709	DO jj = 1, jpjm1
	710	DO ji = 1, fs_jpim1
[10475]	711	pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) )
[10425]	712	END DO
[10413]	713	END DO
	714	END DO
[11081]	715	IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
[10413]	716
	717	ELSE !== even ice time step: adv_y then adv_x ==!
	718	!
[10439]	719	DO jl = 1, jpl !-- flux in y-direction
[10425]	720	DO jj = 1, jpjm1
	721	DO ji = 1, fs_jpim1
[10475]	722	pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
[10425]	723	END DO
[10413]	724	END DO
	725	END DO
[11081]	726	IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
[10413]	727	!
[10439]	728	DO jl = 1, jpl !-- first guess of tracer from v-flux
[10425]	729	DO jj = 2, jpjm1
[10439]	730	DO ji = fs_2, fs_jpim1
[10446]	731	ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) &
	732	& + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk)
	733	!
	734	zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
[10425]	735	END DO
[10413]	736	END DO
	737	END DO
[10446]	738	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
[10413]	739	!
[10439]	740	DO jl = 1, jpl !-- flux in x-direction
[10425]	741	DO jj = 1, jpjm1
	742	DO ji = 1, fs_jpim1
[10475]	743	pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) )
[10425]	744	END DO
[10413]	745	END DO
	746	END DO
[11081]	747	IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
[10413]	748
	749	ENDIF
[11081]	750	IF( np_limiter == 1 ) CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
[10413]	751
	752	ENDIF
	753
	754	END SUBROUTINE cen2
	755
	756
[10475]	757	SUBROUTINE macho( pamsk, kn_umx, jt, kt, pdt, pt, pu, pv, pubox, pvbox, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
[10413]	758	!!---------------------------------------------------------------------
	759	!! * ROUTINE macho *
	760	!!
[10446]	761	!! ** Purpose : compute the high order fluxes using Ultimate-Macho scheme
[10413]	762	!!
[10446]	763	!! ** Method : ...
[10413]	764	!!
	765	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
	766	!!----------------------------------------------------------------------
[10439]	767	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
	768	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
	769	INTEGER , INTENT(in ) :: jt ! number of sub-iteration
	770	INTEGER , INTENT(in ) :: kt ! number of iteration
	771	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
	772	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields
	773	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu, pv ! 2 ice velocity components
	774	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pubox, pvbox ! upstream velocity
[10446]	775	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt_ups ! upstream guess of tracer
	776	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pfu_ups, pfv_ups ! upstream fluxes
[10425]	777	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes
[10413]	778	!
[10425]	779	INTEGER :: ji, jj, jl ! dummy loop indices
[10446]	780	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zt_u, zt_v, zpt
[10413]	781	!!----------------------------------------------------------------------
	782	!
	783	IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==!
[8586]	784	!
[10413]	785	! !-- ultimate interpolation of pt at u-point --!
[11081]	786	CALL ultimate_x( pamsk, kn_umx, pdt, pt, pu, zt_u, pfu_ho )
[10413]	787	! !-- limiter in x --!
[11081]	788	IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
[10446]	789	! !-- advective form update in zpt --!
[10439]	790	DO jl = 1, jpl
	791	DO jj = 2, jpjm1
	792	DO ji = fs_2, fs_jpim1
[10446]	793	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) &
	794	& + pt (ji,jj,jl) * ( pu (ji,jj ) - pu (ji-1,jj ) ) * r1_e1e2t(ji,jj) &
	795	& * pamsk &
[11081]	796	& ) * pdt ) * tmask(ji,jj,1)
[10413]	797	END DO
[8586]	798	END DO
[10439]	799	END DO
[10446]	800	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
[8586]	801	!
[10413]	802	! !-- ultimate interpolation of pt at v-point --!
	803	IF( ll_hoxy ) THEN
[11081]	804	CALL ultimate_y( pamsk, kn_umx, pdt, zpt, pv, zt_v, pfv_ho )
[10413]	805	ELSE
[11081]	806	CALL ultimate_y( pamsk, kn_umx, pdt, pt , pv, zt_v, pfv_ho )
[10413]	807	ENDIF
	808	! !-- limiter in y --!
[11081]	809	IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
[10413]	810	!
	811	!
	812	ELSE !== even ice time step: adv_y then adv_x ==!
	813	!
	814	! !-- ultimate interpolation of pt at v-point --!
[11081]	815	CALL ultimate_y( pamsk, kn_umx, pdt, pt, pv, zt_v, pfv_ho )
[10413]	816	! !-- limiter in y --!
[11081]	817	IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
[10446]	818	! !-- advective form update in zpt --!
[10439]	819	DO jl = 1, jpl
	820	DO jj = 2, jpjm1
	821	DO ji = fs_2, fs_jpim1
[10446]	822	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) &
	823	& + pt (ji,jj,jl) * ( pv (ji,jj ) - pv (ji,jj-1 ) ) * r1_e1e2t(ji,jj) &
	824	& * pamsk &
	825	& ) * pdt ) * tmask(ji,jj,1)
[10425]	826	END DO
[10413]	827	END DO
[10439]	828	END DO
[10446]	829	CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
[10413]	830	!
	831	! !-- ultimate interpolation of pt at u-point --!
	832	IF( ll_hoxy ) THEN
[11081]	833	CALL ultimate_x( pamsk, kn_umx, pdt, zpt, pu, zt_u, pfu_ho )
[10413]	834	ELSE
[11081]	835	CALL ultimate_x( pamsk, kn_umx, pdt, pt , pu, zt_u, pfu_ho )
[10413]	836	ENDIF
	837	! !-- limiter in x --!
[11081]	838	IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
[10413]	839	!
	840	ENDIF
	841
[11081]	842	IF( np_limiter == 1 ) CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
[8586]	843	!
	844	END SUBROUTINE macho
	845
	846
[11081]	847	SUBROUTINE ultimate_x( pamsk, kn_umx, pdt, pt, pu, pt_u, pfu_ho )
[8586]	848	!!---------------------------------------------------------------------
	849	!! * ROUTINE ultimate_x *
	850	!!
[10446]	851	!! ** Purpose : compute tracer at u-points
[8586]	852	!!
[10446]	853	!! ** Method : ...
[8586]	854	!!
	855	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
	856	!!----------------------------------------------------------------------
[11081]	857	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
[10439]	858	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
	859	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
	860	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu ! ice i-velocity component
	861	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields
[10425]	862	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_u ! tracer at u-point
	863	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho ! high order flux
[8586]	864	!
[10425]	865	INTEGER :: ji, jj, jl ! dummy loop indices
[10439]	866	REAL(wp) :: zcu, zdx2, zdx4 ! - -
[10425]	867	REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztu1, ztu2, ztu3, ztu4
[8586]	868	!!----------------------------------------------------------------------
	869	!
	870	! !-- Laplacian in i-direction --!
[10425]	871	DO jl = 1, jpl
	872	DO jj = 2, jpjm1 ! First derivative (gradient)
	873	DO ji = 1, fs_jpim1
	874	ztu1(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
	875	END DO
	876	! ! Second derivative (Laplacian)
	877	DO ji = fs_2, fs_jpim1
	878	ztu2(ji,jj,jl) = ( ztu1(ji,jj,jl) - ztu1(ji-1,jj,jl) ) * r1_e1t(ji,jj)
	879	END DO
[8586]	880	END DO
	881	END DO
[10425]	882	CALL lbc_lnk( 'icedyn_adv_umx', ztu2, 'T', 1. )
[8586]	883	!
	884	! !-- BiLaplacian in i-direction --!
[10425]	885	DO jl = 1, jpl
	886	DO jj = 2, jpjm1 ! Third derivative
	887	DO ji = 1, fs_jpim1
	888	ztu3(ji,jj,jl) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
	889	END DO
	890	! ! Fourth derivative
	891	DO ji = fs_2, fs_jpim1
	892	ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj)
	893	END DO
[8586]	894	END DO
	895	END DO
[10425]	896	CALL lbc_lnk( 'icedyn_adv_umx', ztu4, 'T', 1. )
[8586]	897	!
	898	!
[10413]	899	SELECT CASE (kn_umx )
[8586]	900	!
	901	CASE( 1 ) !== 1st order central TIM ==! (Eq. 21)
	902	!
[10425]	903	DO jl = 1, jpl
	904	DO jj = 1, jpjm1
	905	DO ji = 1, fs_jpim1 ! vector opt.
[10475]	906	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) &
	907	& - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
[10425]	908	END DO
[8586]	909	END DO
	910	END DO
	911	!
	912	CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23)
	913	!
[10425]	914	DO jl = 1, jpl
	915	DO jj = 1, jpjm1
	916	DO ji = 1, fs_jpim1 ! vector opt.
	917	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
	918	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) &
[10475]	919	& - zcu * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
[10425]	920	END DO
[8586]	921	END DO
	922	END DO
	923	!
	924	CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24)
	925	!
[10425]	926	DO jl = 1, jpl
	927	DO jj = 1, jpjm1
	928	DO ji = 1, fs_jpim1 ! vector opt.
	929	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
	930	zdx2 = e1u(ji,jj) * e1u(ji,jj)
[10439]	931	!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
[10446]	932	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
[10475]	933	& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
	934	& + z1_6 * zdx2 * ( zcuzcu - 1._wp ) ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
[10446]	935	& - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
[10425]	936	END DO
[8586]	937	END DO
	938	END DO
	939	!
	940	CASE( 4 ) !== 4th order central TIM ==! (Eq. 27)
	941	!
[10425]	942	DO jl = 1, jpl
	943	DO jj = 1, jpjm1
	944	DO ji = 1, fs_jpim1 ! vector opt.
	945	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
	946	zdx2 = e1u(ji,jj) * e1u(ji,jj)
[10439]	947	!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
[10475]	948	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
	949	& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
	950	& + z1_6 * zdx2 * ( zcuzcu - 1._wp ) ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
	951	& - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
[10425]	952	END DO
[8586]	953	END DO
	954	END DO
	955	!
	956	CASE( 5 ) !== 5th order central TIM ==! (Eq. 29)
	957	!
[10425]	958	DO jl = 1, jpl
	959	DO jj = 1, jpjm1
	960	DO ji = 1, fs_jpim1 ! vector opt.
	961	zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
	962	zdx2 = e1u(ji,jj) * e1u(ji,jj)
[10439]	963	!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
[10425]	964	zdx4 = zdx2 * zdx2
[10475]	965	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
	966	& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
	967	& + z1_6 * zdx2 * ( zcuzcu - 1._wp ) ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
	968	& - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) &
[10446]	969	& + z1_120 * zdx4 * ( zcuzcu - 1._wp ) ( zcuzcu - 4._wp ) ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl) &
[10425]	970	& - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) )
	971	END DO
[8586]	972	END DO
	973	END DO
	974	!
	975	END SELECT
[10439]	976	!
	977	! if pt at u-point is negative then use the upstream value
	978	! this should not be necessary if a proper sea-ice mask is set in Ultimate
	979	! to degrade the order of the scheme when necessary (for ex. at the ice edge)
[10413]	980	IF( ll_neg ) THEN
[10425]	981	DO jl = 1, jpl
	982	DO jj = 1, jpjm1
	983	DO ji = 1, fs_jpim1
[11081]	984	IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN
[10475]	985	pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) &
	986	& - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
[10425]	987	ENDIF
	988	END DO
[10413]	989	END DO
	990	END DO
	991	ENDIF
[10439]	992	! !-- High order flux in i-direction --!
[10425]	993	DO jl = 1, jpl
	994	DO jj = 1, jpjm1
	995	DO ji = 1, fs_jpim1 ! vector opt.
[10439]	996	pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl)
[10425]	997	END DO
[10413]	998	END DO
	999	END DO
[8586]	1000	!
	1001	END SUBROUTINE ultimate_x
	1002
	1003
[11081]	1004	SUBROUTINE ultimate_y( pamsk, kn_umx, pdt, pt, pv, pt_v, pfv_ho )
[8586]	1005	!!---------------------------------------------------------------------
	1006	!! * ROUTINE ultimate_y *
	1007	!!
[10446]	1008	!! ** Purpose : compute tracer at v-points
[8586]	1009	!!
[10446]	1010	!! ** Method : ...
[8586]	1011	!!
	1012	!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
	1013	!!----------------------------------------------------------------------
[11081]	1014	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
[10439]	1015	INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
	1016	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
	1017	REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pv ! ice j-velocity component
	1018	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields
[10425]	1019	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_v ! tracer at v-point
	1020	REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfv_ho ! high order flux
[8586]	1021	!
[10439]	1022	INTEGER :: ji, jj, jl ! dummy loop indices
[8586]	1023	REAL(wp) :: zcv, zdy2, zdy4 ! - -
[10425]	1024	REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztv1, ztv2, ztv3, ztv4
[8586]	1025	!!----------------------------------------------------------------------
	1026	!
	1027	! !-- Laplacian in j-direction --!
[10425]	1028	DO jl = 1, jpl
	1029	DO jj = 1, jpjm1 ! First derivative (gradient)
	1030	DO ji = fs_2, fs_jpim1
	1031	ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
	1032	END DO
[8586]	1033	END DO
[10425]	1034	DO jj = 2, jpjm1 ! Second derivative (Laplacian)
	1035	DO ji = fs_2, fs_jpim1
	1036	ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj)
	1037	END DO
[8586]	1038	END DO
	1039	END DO
[10425]	1040	CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1. )
[8586]	1041	!
	1042	! !-- BiLaplacian in j-direction --!
[10425]	1043	DO jl = 1, jpl
	1044	DO jj = 1, jpjm1 ! First derivative
	1045	DO ji = fs_2, fs_jpim1
	1046	ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
	1047	END DO
[8586]	1048	END DO
[10425]	1049	DO jj = 2, jpjm1 ! Second derivative
	1050	DO ji = fs_2, fs_jpim1
	1051	ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj)
	1052	END DO
[8586]	1053	END DO
	1054	END DO
[10425]	1055	CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1. )
[8586]	1056	!
	1057	!
[10413]	1058	SELECT CASE (kn_umx )
[10425]	1059	!
[8586]	1060	CASE( 1 ) !== 1st order central TIM ==! (Eq. 21)
[10425]	1061	DO jl = 1, jpl
	1062	DO jj = 1, jpjm1
	1063	DO ji = 1, fs_jpim1
[10475]	1064	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) &
	1065	& - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
[10425]	1066	END DO
[8586]	1067	END DO
	1068	END DO
	1069	!
	1070	CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23)
[10425]	1071	DO jl = 1, jpl
	1072	DO jj = 1, jpjm1
	1073	DO ji = 1, fs_jpim1
	1074	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
[10475]	1075	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) &
	1076	& - zcv * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
[10425]	1077	END DO
[8586]	1078	END DO
	1079	END DO
	1080	!
	1081	CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24)
[10425]	1082	DO jl = 1, jpl
	1083	DO jj = 1, jpjm1
	1084	DO ji = 1, fs_jpim1
	1085	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
	1086	zdy2 = e2v(ji,jj) * e2v(ji,jj)
[10439]	1087	!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
[10475]	1088	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
	1089	& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
[10446]	1090	& + z1_6 * zdy2 * ( zcvzcv - 1._wp ) ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
[10425]	1091	& - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
	1092	END DO
[8586]	1093	END DO
	1094	END DO
	1095	!
	1096	CASE( 4 ) !== 4th order central TIM ==! (Eq. 27)
[10425]	1097	DO jl = 1, jpl
	1098	DO jj = 1, jpjm1
	1099	DO ji = 1, fs_jpim1
	1100	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
	1101	zdy2 = e2v(ji,jj) * e2v(ji,jj)
[10439]	1102	!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
[10475]	1103	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
	1104	& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
	1105	& + z1_6 * zdy2 * ( zcvzcv - 1._wp ) ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
	1106	& - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
[10425]	1107	END DO
[8586]	1108	END DO
	1109	END DO
	1110	!
	1111	CASE( 5 ) !== 5th order central TIM ==! (Eq. 29)
[10425]	1112	DO jl = 1, jpl
	1113	DO jj = 1, jpjm1
	1114	DO ji = 1, fs_jpim1
	1115	zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
	1116	zdy2 = e2v(ji,jj) * e2v(ji,jj)
[10439]	1117	!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
[10425]	1118	zdy4 = zdy2 * zdy2
[10446]	1119	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
[10475]	1120	& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
	1121	& + z1_6 * zdy2 * ( zcvzcv - 1._wp ) ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
	1122	& - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) &
[10446]	1123	& + z1_120 * zdy4 * ( zcvzcv - 1._wp ) ( zcvzcv - 4._wp ) ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl) &
[10425]	1124	& - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) )
	1125	END DO
[8586]	1126	END DO
	1127	END DO
	1128	!
	1129	END SELECT
[10439]	1130	!
	1131	! if pt at v-point is negative then use the upstream value
	1132	! this should not be necessary if a proper sea-ice mask is set in Ultimate
	1133	! to degrade the order of the scheme when necessary (for ex. at the ice edge)
[10413]	1134	IF( ll_neg ) THEN
[10425]	1135	DO jl = 1, jpl
	1136	DO jj = 1, jpjm1
	1137	DO ji = 1, fs_jpim1
[11081]	1138	IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN
[10475]	1139	pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) &
	1140	& - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
[10425]	1141	ENDIF
	1142	END DO
[10413]	1143	END DO
	1144	END DO
	1145	ENDIF
[10439]	1146	! !-- High order flux in j-direction --!
[10425]	1147	DO jl = 1, jpl
	1148	DO jj = 1, jpjm1
	1149	DO ji = 1, fs_jpim1 ! vector opt.
[10439]	1150	pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl)
[10425]	1151	END DO
[10413]	1152	END DO
	1153	END DO
[8586]	1154	!
	1155	END SUBROUTINE ultimate_y
[10413]	1156
	1157
[10519]	1158	SUBROUTINE nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
[8586]	1159	!!---------------------------------------------------------------------
[10519]	1160	!! * ROUTINE nonosc_ice *
[8586]	1161	!!
[10446]	1162	!! ** Purpose : compute monotonic tracer fluxes from the upstream
[10519]	1163	!! scheme and the before field by a non-oscillatory algorithm
[8586]	1164	!!
[10446]	1165	!! ** Method : ...
[8586]	1166	!!----------------------------------------------------------------------
[10439]	1167	REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)
	1168	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
[10425]	1169	REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pu ! ice i-velocity => u*e2
	1170	REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pv ! ice j-velocity => v*e1
[10446]	1171	REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pt, pt_ups ! before field & upstream guess of after field
	1172	REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pfv_ups, pfu_ups ! upstream flux
[10425]	1173	REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfv_ho, pfu_ho ! monotonic flux
[8586]	1174	!
[10425]	1175	INTEGER :: ji, jj, jl ! dummy loop indices
[10475]	1176	REAL(wp) :: zpos, zneg, zbig, zup, zdo, z1_dt ! local scalars
	1177	REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zcoef, zzt ! - -
[10425]	1178	REAL(wp), DIMENSION(jpi,jpj ) :: zbup, zbdo
[10439]	1179	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zbetup, zbetdo, zti_ups, ztj_ups
[8586]	1180	!!----------------------------------------------------------------------
	1181	zbig = 1.e+40_wp
[10425]	1182
[10413]	1183	! antidiffusive flux : high order minus low order
	1184	! --------------------------------------------------
[10425]	1185	DO jl = 1, jpl
	1186	DO jj = 1, jpjm1
	1187	DO ji = 1, fs_jpim1 ! vector opt.
[10439]	1188	pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)
	1189	pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)
[10425]	1190	END DO
	1191	END DO
[8586]	1192	END DO
	1193
[10413]	1194	! extreme case where pfu_ho has to be zero
	1195	! ----------------------------------------
	1196	! pfu_ho
	1197	! * --->
	1198	! \| \| * \| \|
	1199	! \| \| \| * \|
	1200	! \| \| \| \| *
[10439]	1201	! t_ups : i-1 i i+1 i+2
[11081]	1202	IF( ll_prelim ) THEN
[10413]	1203
[10425]	1204	DO jl = 1, jpl
	1205	DO jj = 2, jpjm1
	1206	DO ji = fs_2, fs_jpim1
[10439]	1207	zti_ups(ji,jj,jl)= pt_ups(ji+1,jj ,jl)
	1208	ztj_ups(ji,jj,jl)= pt_ups(ji ,jj+1,jl)
[10425]	1209	END DO
[10413]	1210	END DO
	1211	END DO
[10439]	1212	CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_ups, 'T', 1., ztj_ups, 'T', 1. )
[8586]	1213
[10425]	1214	DO jl = 1, jpl
	1215	DO jj = 2, jpjm1
	1216	DO ji = fs_2, fs_jpim1
[10475]	1217	IF ( pfu_ho(ji,jj,jl) * ( pt_ups(ji+1,jj ,jl) - pt_ups(ji,jj,jl) ) <= 0._wp .AND. &
	1218	& pfv_ho(ji,jj,jl) * ( pt_ups(ji ,jj+1,jl) - pt_ups(ji,jj,jl) ) <= 0._wp ) THEN
[10425]	1219	!
[10475]	1220	IF( pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj ,jl) - zti_ups(ji,jj,jl) ) <= 0._wp .AND. &
	1221	& pfv_ho(ji,jj,jl) * ( ztj_ups(ji ,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0._wp ) THEN
	1222	pfu_ho(ji,jj,jl)=0._wp
	1223	pfv_ho(ji,jj,jl)=0._wp
[10439]	1224	ENDIF
[10425]	1225	!
[10475]	1226	IF( pfu_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji-1,jj ,jl) ) <= 0._wp .AND. &
	1227	& pfv_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji ,jj-1,jl) ) <= 0._wp ) THEN
	1228	pfu_ho(ji,jj,jl)=0._wp
	1229	pfv_ho(ji,jj,jl)=0._wp
[10439]	1230	ENDIF
[10425]	1231	!
	1232	ENDIF
	1233	END DO
[10413]	1234	END DO
	1235	END DO
[10425]	1236	CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond.
[10413]	1237
	1238	ENDIF
[10425]	1239
[8586]	1240	! Search local extrema
	1241	! --------------------
[10439]	1242	! max/min of pt & pt_ups with large negative/positive value (-/+zbig) outside ice cover
[10425]	1243	z1_dt = 1._wp / pdt
	1244	DO jl = 1, jpl
	1245
	1246	DO jj = 1, jpj
	1247	DO ji = 1, jpi
[10439]	1248	IF ( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN
[10425]	1249	zbup(ji,jj) = -zbig
	1250	zbdo(ji,jj) = zbig
[10439]	1251	ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) > 0._wp ) THEN
	1252	zbup(ji,jj) = pt_ups(ji,jj,jl)
	1253	zbdo(ji,jj) = pt_ups(ji,jj,jl)
	1254	ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN
[10425]	1255	zbup(ji,jj) = pt(ji,jj,jl)
	1256	zbdo(ji,jj) = pt(ji,jj,jl)
	1257	ELSE
[10439]	1258	zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_ups(ji,jj,jl) )
	1259	zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_ups(ji,jj,jl) )
[10425]	1260	ENDIF
	1261	END DO
[10413]	1262	END DO
[8586]	1263
[10425]	1264	DO jj = 2, jpjm1
	1265	DO ji = fs_2, fs_jpim1 ! vector opt.
	1266	!
[10475]	1267	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
	1268	zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj), zbdo(ji+1,jj), zbdo(ji,jj-1), zbdo(ji,jj+1) )
[10425]	1269	!
[10475]	1270	zpos = MAX( 0._wp, pfu_ho(ji-1,jj ,jl) ) - MIN( 0._wp, pfu_ho(ji ,jj ,jl) ) & ! positive/negative part of the flux
	1271	& + MAX( 0._wp, pfv_ho(ji ,jj-1,jl) ) - MIN( 0._wp, pfv_ho(ji ,jj ,jl) )
	1272	zneg = MAX( 0._wp, pfu_ho(ji ,jj ,jl) ) - MIN( 0._wp, pfu_ho(ji-1,jj ,jl) ) &
	1273	& + MAX( 0._wp, pfv_ho(ji ,jj ,jl) ) - MIN( 0._wp, pfv_ho(ji ,jj-1,jl) )
[10425]	1274	!
[10475]	1275	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]	1276	& ) * ( 1. - pamsk )
[10475]	1277	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]	1278	& ) * ( 1. - pamsk )
[10425]	1279	!
	1280	! ! up & down beta terms
[11081]	1281	! clem: zbetup and zbetdo must be 0 for zpos>1.e-10 & zneg>1.e-10 (do not put 0 instead of 1.e-10 !!!)
	1282	IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt
	1283	ELSE ; zbetup(ji,jj,jl) = 0._wp ! zbig
[10425]	1284	ENDIF
	1285	!
[11081]	1286	IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
	1287	ELSE ; zbetdo(ji,jj,jl) = 0._wp ! zbig
[10425]	1288	ENDIF
	1289	!
	1290	! if all the points are outside ice cover
[10475]	1291	IF( zup == -zbig ) zbetup(ji,jj,jl) = 0._wp ! zbig
	1292	IF( zdo == zbig ) zbetdo(ji,jj,jl) = 0._wp ! zbig
[10425]	1293	!
	1294	END DO
[8586]	1295	END DO
	1296	END DO
[10425]	1297	CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)
[8586]	1298
[10413]	1299
	1300	! monotonic flux in the y direction
	1301	! ---------------------------------
[10425]	1302	DO jl = 1, jpl
	1303	DO jj = 1, jpjm1
	1304	DO ji = 1, fs_jpim1 ! vector opt.
	1305	zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) )
	1306	zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) )
[10475]	1307	zcu = 0.5_wp + SIGN( 0.5_wp , pfu_ho(ji,jj,jl) )
[10425]	1308	!
	1309	zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu )
[10439]	1310	!
	1311	pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_ups(ji,jj,jl)
	1312	!
[10425]	1313	END DO
[8637]	1314	END DO
[10413]	1315
[10425]	1316	DO jj = 1, jpjm1
	1317	DO ji = 1, fs_jpim1 ! vector opt.
	1318	zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) )
	1319	zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) )
[10475]	1320	zcv = 0.5_wp + SIGN( 0.5_wp , pfv_ho(ji,jj,jl) )
[10425]	1321	!
	1322	zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv )
[10439]	1323	!
	1324	pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_ups(ji,jj,jl)
	1325	!
[10425]	1326	END DO
[8586]	1327	END DO
[10413]	1328
	1329	END DO
[8586]	1330	!
[10519]	1331	END SUBROUTINE nonosc_ice
[8586]	1332
[10446]	1333
	1334	SUBROUTINE limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )
[10413]	1335	!!---------------------------------------------------------------------
	1336	!! * ROUTINE limiter_x *
	1337	!!
	1338	!! ** Purpose : compute flux limiter
	1339	!!----------------------------------------------------------------------
[10446]	1340	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
	1341	REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pu ! ice i-velocity => u*e2
	1342	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt ! ice tracer
	1343	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pfu_ups ! upstream flux
	1344	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pfu_ho ! high order flux
[10413]	1345	!
	1346	REAL(wp) :: Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr
[10425]	1347	INTEGER :: ji, jj, jl ! dummy loop indices
	1348	REAL(wp), DIMENSION (jpi,jpj,jpl) :: zslpx ! tracer slopes
[10413]	1349	!!----------------------------------------------------------------------
	1350	!
[10425]	1351	DO jl = 1, jpl
	1352	DO jj = 2, jpjm1
	1353	DO ji = fs_2, fs_jpim1 ! vector opt.
	1354	zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1)
	1355	END DO
[10413]	1356	END DO
	1357	END DO
[10425]	1358	CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.) ! lateral boundary cond.
[10413]	1359
[10425]	1360	DO jl = 1, jpl
	1361	DO jj = 2, jpjm1
	1362	DO ji = fs_2, fs_jpim1 ! vector opt.
	1363	uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj)
	1364
	1365	Rjm = zslpx(ji-1,jj,jl)
	1366	Rj = zslpx(ji ,jj,jl)
	1367	Rjp = zslpx(ji+1,jj,jl)
[10413]	1368
[11081]	1369	IF( np_limiter == 3 ) THEN
[10413]	1370
[10425]	1371	IF( pu(ji,jj) > 0. ) THEN ; Rr = Rjm
	1372	ELSE ; Rr = Rjp
	1373	ENDIF
[10413]	1374
[10425]	1375	zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)
	1376	IF( Rj > 0. ) THEN
	1377	zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pu(ji,jj)), &
	1378	& MIN( 2. * Rr * 0.5 * ABS(pu(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
	1379	ELSE
	1380	zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)), &
	1381	& MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
	1382	ENDIF
	1383	pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter
[10413]	1384
[11081]	1385	ELSEIF( np_limiter == 2 ) THEN
[10425]	1386	IF( Rj /= 0. ) THEN
	1387	IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj
	1388	ELSE ; Cr = Rjp / Rj
	1389	ENDIF
	1390	ELSE
	1391	Cr = 0.
[10413]	1392	ENDIF
[10425]	1393
	1394	! -- superbee --
	1395	zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
	1396	! -- van albada 2 --
	1397	!!zpsi = 2.Cr / (CrCr+1.)
	1398	! -- sweby (with beta=1) --
	1399	!!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
	1400	! -- van Leer --
	1401	!!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
	1402	! -- ospre --
	1403	!!zpsi = 1.5 * ( CrCr + Cr ) / ( CrCr + Cr + 1. )
	1404	! -- koren --
	1405	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( (1.+2Cr)/3., 2. ) ) )
	1406	! -- charm --
	1407	!IF( Cr > 0. ) THEN ; zpsi = Cr * (3.Cr + 1.) / ( (Cr + 1.) (Cr + 1.) )
	1408	!ELSE ; zpsi = 0.
[10413]	1409	!ENDIF
[10425]	1410	! -- van albada 1 --
	1411	!!zpsi = (CrCr + Cr) / (CrCr +1)
	1412	! -- smart --
	1413	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, 4. ) ) )
	1414	! -- umist --
	1415	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
[10413]	1416
[10425]	1417	! high order flux corrected by the limiter
	1418	pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFLzpsi ) Rj * 0.5
[10413]	1419
[10425]	1420	ENDIF
	1421	END DO
[10413]	1422	END DO
	1423	END DO
[10425]	1424	CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.) ! lateral boundary cond.
[10413]	1425	!
	1426	END SUBROUTINE limiter_x
	1427
[10446]	1428
	1429	SUBROUTINE limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )
[10413]	1430	!!---------------------------------------------------------------------
	1431	!! * ROUTINE limiter_y *
	1432	!!
	1433	!! ** Purpose : compute flux limiter
	1434	!!----------------------------------------------------------------------
[10446]	1435	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
	1436	REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pv ! ice i-velocity => u*e2
	1437	REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pt ! ice tracer
	1438	REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pfv_ups ! upstream flux
	1439	REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfv_ho ! high order flux
[10413]	1440	!
	1441	REAL(wp) :: Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr
[10425]	1442	INTEGER :: ji, jj, jl ! dummy loop indices
	1443	REAL(wp), DIMENSION (jpi,jpj,jpl) :: zslpy ! tracer slopes
[10413]	1444	!!----------------------------------------------------------------------
	1445	!
[10425]	1446	DO jl = 1, jpl
	1447	DO jj = 2, jpjm1
	1448	DO ji = fs_2, fs_jpim1 ! vector opt.
	1449	zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1)
	1450	END DO
[10413]	1451	END DO
	1452	END DO
[10425]	1453	CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.) ! lateral boundary cond.
[10413]	1454
[10425]	1455	DO jl = 1, jpl
	1456	DO jj = 2, jpjm1
	1457	DO ji = fs_2, fs_jpim1 ! vector opt.
	1458	vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj)
[10413]	1459
[10425]	1460	Rjm = zslpy(ji,jj-1,jl)
	1461	Rj = zslpy(ji,jj ,jl)
	1462	Rjp = zslpy(ji,jj+1,jl)
[10413]	1463
[11081]	1464	IF( np_limiter == 3 ) THEN
[10413]	1465
[10425]	1466	IF( pv(ji,jj) > 0. ) THEN ; Rr = Rjm
	1467	ELSE ; Rr = Rjp
	1468	ENDIF
[10413]	1469
[10425]	1470	zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)
	1471	IF( Rj > 0. ) THEN
	1472	zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pv(ji,jj)), &
	1473	& MIN( 2. * Rr * 0.5 * ABS(pv(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
	1474	ELSE
	1475	zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)), &
	1476	& MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
	1477	ENDIF
	1478	pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter
[10413]	1479
[11081]	1480	ELSEIF( np_limiter == 2 ) THEN
[10413]	1481
[10425]	1482	IF( Rj /= 0. ) THEN
	1483	IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj
	1484	ELSE ; Cr = Rjp / Rj
	1485	ENDIF
	1486	ELSE
	1487	Cr = 0.
	1488	ENDIF
[10413]	1489
[10425]	1490	! -- superbee --
	1491	zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
	1492	! -- van albada 2 --
	1493	!!zpsi = 2.Cr / (CrCr+1.)
	1494	! -- sweby (with beta=1) --
	1495	!!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
	1496	! -- van Leer --
	1497	!!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
	1498	! -- ospre --
	1499	!!zpsi = 1.5 * ( CrCr + Cr ) / ( CrCr + Cr + 1. )
	1500	! -- koren --
	1501	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( (1.+2Cr)/3., 2. ) ) )
	1502	! -- charm --
	1503	!IF( Cr > 0. ) THEN ; zpsi = Cr * (3.Cr + 1.) / ( (Cr + 1.) (Cr + 1.) )
	1504	!ELSE ; zpsi = 0.
	1505	!ENDIF
	1506	! -- van albada 1 --
	1507	!!zpsi = (CrCr + Cr) / (CrCr +1)
	1508	! -- smart --
	1509	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, 4. ) ) )
	1510	! -- umist --
	1511	!!zpsi = MAX( 0., MIN( 2.Cr, MIN( 0.25+0.75Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
[10413]	1512
[10425]	1513	! high order flux corrected by the limiter
	1514	pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFLzpsi ) Rj * 0.5
	1515
	1516	ENDIF
	1517	END DO
[10413]	1518	END DO
	1519	END DO
[10425]	1520	CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.) ! lateral boundary cond.
[10413]	1521	!
	1522	END SUBROUTINE limiter_y
	1523
[11081]	1524
	1525	SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
	1526	!!-------------------------------------------------------------------
	1527	!! * ROUTINE Hbig *
	1528	!!
	1529	!! ** Purpose : Thickness correction in case advection scheme creates
[12379]	1530	!! abnormally thick ice or snow
[11081]	1531	!!
	1532	!! ** Method : 1- check whether ice thickness is larger than the surrounding 9-points
	1533	!! (before advection) and reduce it by adapting ice concentration
	1534	!! 2- check whether snow thickness is larger than the surrounding 9-points
	1535	!! (before advection) and reduce it by sending the excess in the ocean
	1536	!! 3- check whether snow load deplets the snow-ice interface below sea level$
	1537	!! and reduce it by sending the excess in the ocean
	1538	!! 4- correct pond fraction to avoid a_ip > a_i
	1539	!!
	1540	!! ** input : Max thickness of the surrounding 9-points
	1541	!!-------------------------------------------------------------------
	1542	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
	1543	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max ! max ice thick from surrounding 9-pts
	1544	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip
	1545	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s
	1546	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i
	1547	!
	1548	INTEGER :: ji, jj, jk, jl ! dummy loop indices
	1549	REAL(wp) :: z1_dt, zhip, zhi, zhs, zvs_excess, zfra
	1550	REAL(wp), DIMENSION(jpi,jpj) :: zswitch
	1551	!!-------------------------------------------------------------------
	1552	!
	1553	z1_dt = 1._wp / pdt
	1554	!
	1555	DO jl = 1, jpl
	1556
	1557	DO jj = 1, jpj
	1558	DO ji = 1, jpi
	1559	IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
	1560	!
	1561	! ! -- check h_ip -- !
	1562	! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
	1563	IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
	1564	zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) )
	1565	IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN
	1566	pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl)
	1567	ENDIF
	1568	ENDIF
	1569	!
	1570	! ! -- check h_i -- !
	1571	! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
	1572	zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl)
	1573	IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
	1574	pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) ) !-- bound h_i to hi_max (99 m)
	1575	ENDIF
	1576	!
	1577	! ! -- check h_s -- !
	1578	! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
	1579	zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl)
	1580	IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
	1581	zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
	1582	!
	1583	wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt
	1584	hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
	1585	!
	1586	pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
	1587	pv_s(ji,jj,jl) = pa_i(ji,jj,jl) * phs_max(ji,jj,jl)
	1588	ENDIF
	1589	!
	1590	! ! -- check snow load -- !
	1591	! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean
	1592	! this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin)
	1593	! this imposed mini can artificially make the snow very thick (if concentration decreases drastically)
	1594	zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
	1595	IF( zvs_excess > 0._wp ) THEN
	1596	zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 )
	1597	wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt
	1598	hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
	1599	!
	1600	pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
	1601	pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess
	1602	ENDIF
	1603
	1604	ENDIF
	1605	END DO
	1606	END DO
	1607	END DO
	1608	! !-- correct pond fraction to avoid a_ip > a_i
	1609	WHERE( pa_ip(:,:,:) > pa_i(:,:,:) ) pa_ip(:,:,:) = pa_i(:,:,:)
	1610	!
	1611	!
	1612	END SUBROUTINE Hbig
	1613
[8586]	1614	#else
	1615	!!----------------------------------------------------------------------
[9570]	1616	!! Default option Dummy module NO SI3 sea-ice model
[8586]	1617	!!----------------------------------------------------------------------
	1618	#endif
	1619
	1620	!!======================================================================
	1621	END MODULE icedyn_adv_umx

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

Download in other formats: