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

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

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

Download in other formats: