Context Navigation

icedyn_adv_umx.F90 @ 12406

Last change on this file since 12406 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Property svn:keywords set to Id

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

Download in other formats: