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/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_adv_umx.F90 @ 12744

Last change on this file since 12744 was 12744, checked in by clem, 4 years ago
make sure all pond lids are set to 0 when not using this option
Property svn:keywords set to `Id`
File size: 85.7 KB

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

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

Download in other formats: