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 @ 13149

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

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

Download in other formats: