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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

icedyn_adv_umx.F90 in NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icedyn_adv_umx.F90 @ 13553

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

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

Download in other formats: