New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

icedyn_adv_umx.F90 in NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE – NEMO

Context Navigation

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

Last change on this file since 14288 was 13877, checked in by mocavero, 4 years ago
Align branch with the trunk@13747
Property svn:keywords set to `Id`
File size: 89.8 KB

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

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

Download in other formats: