Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

icedyn_adv_umx.F90

Last change on this file was 14258, checked in by andmirek, 4 years ago
Ticket #2482: few changes to improve code on GPU
File size: 86.6 KB

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

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

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

Context Navigation

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

Download in other formats: