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/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icedyn_adv_umx.F90 @ 12443

Last change on this file since 12443 was 12443, checked in by davestorkey, 4 years ago
2020/KERNEL-03_Storkey_Coward_RK3_stage2: More variable renaming: atfp -> rn_atfp (use namelist parameter everywhere) rdtbt -> rDt_e nn_baro -> nn_e rn_scal_load -> rn_load rau0 -> rho0
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	!
[12443]	1507	zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rho0-rhoi) * r1_rhos )
[12377]	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.

Download in other formats: