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

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

icedyn_adv_umx.F90 in NEMO/trunk/src/ICE – NEMO

Context Navigation

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

Last change on this file since 15037 was 15037, checked in by smasson, 3 years ago
trunk: small optimization following [15033], #2699
Property svn:keywords set to `Id`
File size: 89.0 KB

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

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

Download in other formats: