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/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/ICE/icedyn_adv_umx.F90 @ 14021

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

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

Download in other formats: