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_r13898_Tiling_Cleanup_MPI3/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13898_Tiling_Cleanup_MPI3/src/ICE/icedyn_adv_umx.F90 @ 13906

Last change on this file since 13906 was 13906, checked in by mocavero, 4 years ago
Merge with dev_r13296_HPC-07_mocavero_mpi3
Property svn:keywords set to `Id`
File size: 90.0 KB

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

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

Download in other formats: