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

Context Navigation

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

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

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

Download in other formats: