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

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

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

Download in other formats: