Context Navigation

icedyn_adv_umx.F90 @ 12399

Last change on this file since 12399 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Property svn:keywords set to Id

File size: 80.3 KB

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

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

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

Context Navigation

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

Download in other formats: