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

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

icedyn_adv_umx.F90 in NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_GPU/src/ICE – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_GPU/src/ICE/icedyn_adv_umx.F90 @ 13153

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

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

Download in other formats: