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_r13296_HPC-07_mocavero_mpi3/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_adv_umx.F90 @ 13630

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

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

Download in other formats: