Context Navigation

icedyn_adv_umx.F90 @ 15812

Last change on this file since 15812 was 15812, checked in by clem, 18 months ago
change option for UMx advection scheme (np_advS=2)
Property svn:keywords set to `Id`
File size: 92.6 KB

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

Download in other formats: