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/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/ICE/icedyn_adv_umx.F90 @ 12178

Last change on this file since 12178 was 12178, checked in by agn, 4 years ago
updated trunk to v 11653
Property svn:keywords set to `Id`
File size: 84.1 KB

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

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

Download in other formats: