Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

icedyn_adv_umx.F90

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

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

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

Context Navigation

source: NEMO/trunk/src/ICE/icedyn_adv_umx.F90

Download in other formats: