Context Navigation

icedyn_adv_umx.F90 @ 12350

Last change on this file since 12350 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.

Property svn:keywords set to Id

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

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

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

Context Navigation

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn_adv_umx.F90 @ 12350

Download in other formats: