Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

icedyn_adv_umx.F90 @ 15814

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

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.3_mirror/src/ICE/icedyn_adv_umx.F90 @ 15814

Download in other formats: