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/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13923_Tiling_Cleanup_MPI3_LoopFusion/src/ICE/icedyn_adv_umx.F90 @ 13946

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

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

Download in other formats: