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/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icedyn_adv_umx.F90 @ 13886

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

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

Download in other formats: