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/releases/r4.0/r4.0-HEAD/src/ICE – NEMO

Context Navigation

source: NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_adv_umx.F90 @ 13566

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

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

Download in other formats: