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

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

icedyn_adv_umx.F90 in NEMO/releases/release-4.0/src/ICE – NEMO

Context Navigation

source: NEMO/releases/release-4.0/src/ICE/icedyn_adv_umx.F90 @ 10929

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

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

Download in other formats: