Context Navigation

icedyn_adv_umx.F90

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

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_fix_topo_minor/src/ICE/icedyn_adv_umx.F90

Download in other formats: