source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_adv_umx.F90 @ 10419

Last change on this file since 10419 was 10419, checked in by smasson, 22 months ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10418, see #2133

  • Property svn:keywords set to Id
File size: 94.2 KB
Line 
1MODULE 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 trend with the 3D advection trends using a TVD scheme
14   !!   ultimate_x(_y)    : compute a tracer value at velocity points using ULTIMATE scheme at various orders
15   !!   macho             : ???
16   !!   nonosc            : compute monotonic tracer fluxes by 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 lib_mpp        ! MPP library
26   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
27   USE lbclnk         ! lateral boundary conditions (or mpp links)
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   ice_dyn_adv_umx   ! called by icedyn_adv.F90
33     
34   REAL(wp) ::   z1_6   = 1._wp /   6._wp   ! =1/6
35   REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120
36
37   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: amaxu, amaxv
38   
39   ! advect H all the way (and get V=H*A at the end)
40   LOGICAL :: ll_thickness = .FALSE.
41   
42   ! look for 9 points around in nonosc limiter 
43   LOGICAL :: ll_9points = .FALSE.  ! false=better h?
44
45   ! use HgradU in nonosc limiter 
46   LOGICAL :: ll_HgradU = .TRUE.   ! no effect?
47
48   ! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme
49   LOGICAL :: ll_neg = .TRUE.      ! keep TRUE
50   
51   ! limit the fluxes
52   LOGICAL :: ll_zeroup1 = .FALSE. ! false ok if Hbig otherwise needed for 2D sinon on a des valeurs de H trop fortes !!
53   LOGICAL :: ll_zeroup2 = .FALSE. ! false ok for 1D, 2D, 3D
54   LOGICAL :: ll_zeroup4 = .FALSE. ! false ok for 1D, 2D, 3D
55   LOGICAL :: ll_zeroup5 = .FALSE. ! false ok for 1D, 2D
56
57   ! fluxes that are limited are u*H, then (u*H)*(ua/u) is used for V (only for nonosc)
58   LOGICAL :: ll_clem   = .TRUE.   ! simpler than rachid and works
59
60   ! First advect H as H*=Hdiv(u), then use H* for H(n+1)=H(n)-div(uH*)
61   LOGICAL :: ll_gurvan = .FALSE.  ! must be false for 1D case !!
62
63   ! First guess as div(uH) (-true-) or Hdiv(u)+ugradH (-false-)
64   LOGICAL :: ll_1stguess_clem = .FALSE. ! better negative values but less good h
65
66   ! advect (or not) open water. If not, retrieve it from advection of A
67   LOGICAL :: ll_ADVopw = .FALSE.  !
68   
69   ! alternate directions for upstream
70   LOGICAL :: ll_upsxy = .TRUE.
71
72   ! alternate directions for high order
73   LOGICAL :: ll_hoxy = .TRUE.
74   
75   ! prelimiter: use it to avoid overshoot in H
76   LOGICAL :: ll_prelimiter_zalesak = .TRUE.  ! from: Zalesak(1979) eq. 14 => true is better for 1D but false is better in 3D (for h and negative values) => pb in x-y?
77   LOGICAL :: ll_prelimiter_devore  = .FALSE.  ! from: Devore eq. 11 (worth than zalesak)
78
79   ! iterate on the limiter (only nonosc)
80   LOGICAL :: ll_limiter_it2 = .FALSE.
81   
82
83   !! * Substitutions
84#  include "vectopt_loop_substitute.h90"
85   !!----------------------------------------------------------------------
86   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
87   !! $Id$
88   !! Software governed by the CeCILL licence     (./LICENSE)
89   !!----------------------------------------------------------------------
90CONTAINS
91
92   SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice,  &
93      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
94      !!----------------------------------------------------------------------
95      !!                  ***  ROUTINE ice_dyn_adv_umx  ***
96      !!
97      !! **  Purpose :   Compute the now trend due to total advection of
98      !!                 tracers and add it to the general trend of tracer equations
99      !!                 using an "Ultimate-Macho" scheme
100      !!
101      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
102      !!----------------------------------------------------------------------
103      INTEGER                     , INTENT(in   ) ::   kn_umx     ! order of the scheme (1-5=UM or 20=CEN2)
104      INTEGER                     , INTENT(in   ) ::   kt         ! time step
105      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity
106      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity
107      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
108      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
109      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
110      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_i      ! salt content
111      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
112      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
113      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
114      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
115      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
116      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
117      !
118      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
119      INTEGER  ::   icycle                  ! number of sub-timestep for the advection
120      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers
121      REAL(wp) ::   zdt
122      REAL(wp), DIMENSION(1)       ::   zcflprv, zcflnow   ! send zcflnow and receive zcflprv
123      REAL(wp), DIMENSION(jpi,jpj) ::   zudy, zvdx
124      REAL(wp), DIMENSION(jpi,jpj) ::   zati1, zati2
125
126
127
128      REAL(wp), DIMENSION(jpi,jpj)     ::   zcu_box, zcv_box
129      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho
130      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai, z1_aip
131      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhvar
132
133      !!----------------------------------------------------------------------
134      !
135      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
136      !
137      !
138      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !
139      !     When needed, the advection split is applied at the next time-step in order to avoid blocking global comm.
140      !     ...this should not affect too much the stability... Was ok on the tests we did...
141      zcflnow(1) =                  MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
142      zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
143     
144      ! non-blocking global communication send zcflnow and receive zcflprv
145      CALL mpp_delay_max( 'icedyn_adv_umx', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 )
146
147      IF( zcflprv(1) > .5 ) THEN   ;   icycle = 2
148      ELSE                         ;   icycle = 1
149      ENDIF
150     
151      zdt = rdt_ice / REAL(icycle)
152
153      ! --- transport --- !
154      zudy(:,:) = pu_ice(:,:) * e2u(:,:)
155      zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
156
157      ! --- define velocity for advection: u*grad(H) --- !
158      DO jj = 2, jpjm1
159         DO ji = fs_2, fs_jpim1
160            IF    ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp
161            ELSEIF( pu_ice(ji,jj)                   >  0._wp ) THEN   ;   zcu_box(ji,jj) = pu_ice(ji-1,jj)
162            ELSE                                                      ;   zcu_box(ji,jj) = pu_ice(ji  ,jj)
163            ENDIF
164
165            IF    ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp
166            ELSEIF( pv_ice(ji,jj)                   >  0._wp ) THEN   ;   zcv_box(ji,jj) = pv_ice(ji,jj-1)
167            ELSE                                                      ;   zcv_box(ji,jj) = pv_ice(ji,jj  )
168            ENDIF
169         END DO
170      END DO
171
172      IF( ll_zeroup2 ) THEN
173         IF(.NOT. ALLOCATED(amaxu))       ALLOCATE(amaxu (jpi,jpj,jpl))
174         IF(.NOT. ALLOCATED(amaxv))       ALLOCATE(amaxv (jpi,jpj,jpl))
175      ENDIF
176      !---------------!
177      !== advection ==!
178      !---------------!
179      DO jt = 1, icycle
180
181!!$         IF( ll_ADVopw ) THEN
182!!$            zamsk = 1._wp
183!!$            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) )   ! Open water area
184!!$            zamsk = 0._wp
185!!$         ELSE
186            zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
187!!$         ENDIF
188         
189         WHERE( pa_i(:,:,:) >= epsi20 )   ;   z1_ai(:,:,:) = 1._wp / pa_i(:,:,:)
190         ELSEWHERE                        ;   z1_ai(:,:,:) = 0.
191         END WHERE
192            !
193         WHERE( pa_ip(:,:,:) >= epsi20 )  ;   z1_aip(:,:,:) = 1._wp / pa_ip(:,:,:)
194         ELSEWHERE                        ;   z1_aip(:,:,:) = 0.
195         END WHERE
196         !
197         IF( ll_zeroup2 ) THEN
198            DO jl = 1, jpl
199               DO jj = 2, jpjm1
200                  DO ji = fs_2, fs_jpim1
201                     amaxu(ji,jj,jl)=MAX( pa_i(ji,jj,jl), pa_i(ji,jj-1,jl), pa_i(ji,jj+1,jl), &
202                        &                                 pa_i(ji+1,jj,jl), pa_i(ji+1,jj-1,jl), pa_i(ji+1,jj+1,jl) )
203                     amaxv(ji,jj,jl)=MAX( pa_i(ji,jj,jl), pa_i(ji-1,jj,jl), pa_i(ji+1,jj,jl), &
204                        &                                 pa_i(ji,jj+1,jl), pa_i(ji-1,jj+1,jl), pa_i(ji+1,jj+1,jl) )
205                  END DO
206               END DO
207            END DO
208            CALL lbc_lnk_multi('icedyn_adv_umx', amaxu, 'T', 1., amaxv, 'T', 1.)
209         ENDIF
210         !
211         DO jl = 1, jpl
212            zua_ho(:,:,jl) = zudy(:,:)
213            zva_ho(:,:,jl) = zvdx(:,:)
214         END DO
215         
216         zamsk = 1._wp
217         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_i, pa_i, zua_ho, zva_ho ) ! Ice area
218         zamsk = 0._wp
219         !
220         IF( ll_thickness ) THEN
221            DO jl = 1, jpl
222               zua_ho(:,:,jl) = zudy(:,:)
223               zva_ho(:,:,jl) = zvdx(:,:)
224            END DO
225         ENDIF
226            !
227         zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:)
228         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_i )    ! Ice volume
229         IF( ll_thickness )   pv_i(:,:,:) = zhvar(:,:,:) * pa_i(:,:,:)
230         !
231         zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:)
232         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_s )    ! Snw volume
233         IF( ll_thickness )   pv_s(:,:,:) = zhvar(:,:,:) * pa_i(:,:,:)
234         !
235         zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:)
236         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, psv_i )    ! Salt content
237         !
238         zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:)
239         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, poa_i )    ! Age content
240         !
241         DO jk = 1, nlay_i
242            zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:)
243            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_i(:,:,jk,:) ) ! Ice heat content
244         END DO
245         !
246         DO jk = 1, nlay_s
247            zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:)
248            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_s(:,:,jk,:) ) ! Snw heat content
249         END DO
250            !
251         IF ( ln_pnd_H12 ) THEN
252            !
253            DO jl = 1, jpl
254               zua_ho(:,:,jl) = zudy(:,:)
255               zva_ho(:,:,jl) = zvdx(:,:)
256            END DO
257           
258            zamsk = 1._wp
259            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_ip, pa_ip, zua_ho, zva_ho ) ! mp fraction
260            zamsk = 0._wp
261            !
262            zhvar(:,:,:) = pv_ip(:,:,:) * z1_ai(:,:,:)
263            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_ip ) ! mp volume
264         ENDIF
265         !
266         !
267!!$         IF( .NOT. ll_ADVopw ) THEN
268            zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
269            DO jj = 2, jpjm1
270               DO ji = fs_2, fs_jpim1
271                  pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                                                  ! Open water area
272                     &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt
273               END DO
274            END DO
275            CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T',  1. )
276!!$         ENDIF
277         !
278      END DO
279      !
280   END SUBROUTINE ice_dyn_adv_umx
281
282   
283   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho )
284      !!----------------------------------------------------------------------
285      !!                  ***  ROUTINE adv_umx  ***
286      !!
287      !! **  Purpose :   Compute the now trend due to total advection of
288      !!       tracers and add it to the general trend of tracer equations
289      !!
290      !! **  Method  :   TVD scheme, i.e. 2nd order centered scheme with
291      !!       corrected flux (monotonic correction)
292      !!       note: - this advection scheme needs a leap-frog time scheme
293      !!
294      !! ** Action : - pt  the after advective tracer
295      !!----------------------------------------------------------------------
296      REAL(wp)                    , INTENT(in   )           ::   pamsk          ! advection of concentration (1) or other tracers (0)
297      INTEGER                     , INTENT(in   )           ::   kn_umx         ! order of the scheme (1-5=UM or 20=CEN2)
298      INTEGER                     , INTENT(in   )           ::   jt             ! number of sub-iteration
299      INTEGER                     , INTENT(in   )           ::   kt             ! number of iteration
300      REAL(wp)                    , INTENT(in   )           ::   pdt            ! tracer time-step
301      REAL(wp), DIMENSION(:,:  ), INTENT(in   )           ::   pu   , pv      ! 2 ice velocity components => u*e2
302      REAL(wp), DIMENSION(:,:,:), INTENT(in   )           ::   puc  , pvc     ! 2 ice velocity components => u*e2 or u*a*e2u
303      REAL(wp), DIMENSION(:,:  ), INTENT(in   )           ::   pubox, pvbox   ! upstream velocity
304      REAL(wp), DIMENSION(:,:,:), INTENT(inout)           ::   pt             ! tracer field
305      REAL(wp), DIMENSION(:,:,:), INTENT(inout)           ::   ptc            ! tracer content field
306      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho ! high order u*a fluxes
307      !
308      INTEGER  ::   ji, jj, jl       ! dummy loop indices 
309      REAL(wp) ::   ztra             ! local scalar
310      INTEGER  ::   kn_limiter = 1   ! 1=nonosc ; 2=superbee ; 3=h3(rachid)
311      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zfu_ho , zfv_ho , zt_u, zt_v, zpt
312      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zfu_ups, zfv_ups, zt_ups   ! only for nonosc
313      !!----------------------------------------------------------------------
314      !
315      !  upstream (_ups) advection with initial mass fluxes
316      ! ---------------------------------------------------
317
318      IF( ll_gurvan .AND. pamsk==0. ) THEN
319         DO jl = 1, jpl
320            DO jj = 2, jpjm1
321               DO ji = fs_2, fs_jpim1
322                  pt(ji,jj,jl) = ( pt (ji,jj,jl) + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj)     &
323                     &                           + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) )   &
324                     &           * tmask(ji,jj,1)
325               END DO
326            END DO
327         END DO
328         CALL lbc_lnk( 'icedyn_adv_umx', pt, 'T', 1. )
329      ENDIF
330
331     
332      IF( .NOT. ll_upsxy ) THEN
333
334         ! fluxes in both x-y directions
335         DO jl = 1, jpl
336            DO jj = 1, jpjm1
337               DO ji = 1, fs_jpim1
338                  IF( ll_clem ) THEN
339                     zfu_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)
340                     zfv_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)
341                  ELSE
342                     zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * pt(ji+1,jj,jl)
343                     zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj+1,jl)
344                  ENDIF
345               END DO
346            END DO
347         END DO
348
349      ELSE
350         !
351         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
352            ! flux in x-direction
353            DO jl = 1, jpl
354               DO jj = 1, jpjm1
355                  DO ji = 1, fs_jpim1
356                     IF( ll_clem ) THEN
357                        zfu_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)
358                     ELSE
359                        zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * pt(ji+1,jj,jl)
360                     ENDIF
361                  END DO
362               END DO
363            END DO
364           
365            ! first guess of tracer content from u-flux
366            DO jl = 1, jpl
367               DO jj = 2, jpjm1
368                  DO ji = fs_2, fs_jpim1   ! vector opt.
369                     IF( ll_clem ) THEN
370                        IF( ll_gurvan ) THEN
371                           zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
372                        ELSE
373                           zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
374                              &            + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
375                              &            ) * tmask(ji,jj,1)
376                        ENDIF
377                     ELSE
378                        zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) &
379                           &         * tmask(ji,jj,1)
380                     ENDIF
381                     !!                  IF( ji==26 .AND. jj==86) THEN
382                     !!                     WRITE(numout,*) '************************'
383                     !!                     WRITE(numout,*) 'zpt upstream',zpt(ji,jj)
384                     !!                  ENDIF
385                  END DO
386               END DO
387            END DO
388            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
389            !
390            ! flux in y-direction
391            DO jl = 1, jpl
392               DO jj = 1, jpjm1
393                  DO ji = 1, fs_jpim1
394                     IF( ll_clem ) THEN
395                        zfv_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)
396                     ELSE
397                        zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * zpt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * zpt(ji,jj+1,jl)
398                     ENDIF
399                  END DO
400               END DO
401            END DO
402         !
403         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
404            ! flux in y-direction
405            DO jl = 1, jpl
406               DO jj = 1, jpjm1
407                  DO ji = 1, fs_jpim1
408                     IF( ll_clem ) THEN
409                        zfv_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)
410                     ELSE
411                        zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj+1,jl)
412                     ENDIF
413                  END DO
414               END DO
415            END DO
416
417            ! first guess of tracer content from v-flux
418            DO jl = 1, jpl
419               DO jj = 2, jpjm1
420                  DO ji = fs_2, fs_jpim1   ! vector opt.
421                     IF( ll_clem ) THEN
422                        IF( ll_gurvan ) THEN
423                           zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
424                        ELSE
425                           zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &
426                              &            + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) &
427                              &            * tmask(ji,jj,1)
428                        ENDIF
429                     ELSE
430                        zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) &
431                           &            * tmask(ji,jj,1)
432                     ENDIF
433                     !!                  IF( ji==26 .AND. jj==86) THEN
434                     !!                     WRITE(numout,*) '************************'
435                     !!                     WRITE(numout,*) 'zpt upstream',zpt(ji,jj)
436                     !!                  ENDIF
437                  END DO
438               END DO
439            END DO
440            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. )
441            !
442            ! flux in x-direction
443            DO jl = 1, jpl
444               DO jj = 1, jpjm1
445                  DO ji = 1, fs_jpim1
446                     IF( ll_clem ) THEN
447                        zfu_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)
448                     ELSE
449                        zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * zpt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * zpt(ji+1,jj,jl)
450                     ENDIF
451                  END DO
452               END DO
453            END DO
454            !
455         ENDIF
456         
457      ENDIF
458
459      IF( ll_clem .AND. kn_limiter /= 1 )  &
460         & CALL ctl_stop( 'STOP', 'icedyn_adv_umx: ll_clem incompatible with limiters other than nonosc' )
461
462      IF( ll_zeroup2 ) THEN
463         DO jl = 1, jpl
464            DO jj = 1, jpjm1
465               DO ji = 1, fs_jpim1   ! vector opt.
466                  IF( amaxu(ji,jj,jl) == 0._wp )   zfu_ups(ji,jj,jl) = 0._wp
467                  IF( amaxv(ji,jj,jl) == 0._wp )   zfv_ups(ji,jj,jl) = 0._wp
468               END DO
469            END DO
470         END DO
471      ENDIF
472
473      ! guess after content field with upstream scheme
474      DO jl = 1, jpl
475         DO jj = 2, jpjm1
476            DO ji = fs_2, fs_jpim1
477               ztra          = - (   zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj  ,jl)   &
478                  &                + zfv_ups(ji,jj,jl) - zfv_ups(ji  ,jj-1,jl) ) * r1_e1e2t(ji,jj)
479               IF( ll_clem ) THEN
480                  IF( ll_gurvan ) THEN
481                     zt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1)
482                  ELSE
483                     zt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + pdt * ztra + ( pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) )   &
484                        &                                            +   pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) ) &
485                        &                                              * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
486                  ENDIF
487               ELSE
488                  zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1)
489               ENDIF
490               !!            IF( ji==26 .AND. jj==86) THEN
491               !!               WRITE(numout,*) '**************************'
492               !!               WRITE(numout,*) 'zt upstream',zt_ups(ji,jj)
493               !!            ENDIF
494            END DO
495         END DO
496      END DO
497      CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. )
498
499      ! High order (_ho) fluxes
500      ! -----------------------
501      SELECT CASE( kn_umx )
502         !
503      CASE ( 20 )                          !== centered second order ==!
504         !
505         CALL cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho,  &
506            &       zt_ups, zfu_ups, zfv_ups )
507         !
508      CASE ( 1:5 )                         !== 1st to 5th order ULTIMATE-MACHO scheme ==!
509         !
510         CALL macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, zt_u, zt_v, zfu_ho, zfv_ho,  &
511            &        zt_ups, zfu_ups, zfv_ups )
512         !
513      END SELECT
514
515      IF( ll_thickness ) THEN
516         ! final trend with corrected fluxes
517         ! ------------------------------------
518         DO jl = 1, jpl
519            DO jj = 2, jpjm1
520               DO ji = fs_2, fs_jpim1
521                  IF( ll_gurvan ) THEN
522                     ztra       = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj) 
523                  ELSE
524                     ztra       = ( - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) )  & 
525                        &           + ( pt(ji,jj,jl) * ( pu(ji,jj) - pu(ji-1,jj) ) * (1.-pamsk) ) &
526                        &           + ( pt(ji,jj,jl) * ( pv(ji,jj) - pv(ji,jj-1) ) * (1.-pamsk) ) ) * r1_e1e2t(ji,jj)
527                  ENDIF
528                  pt(ji,jj,jl) = ( pt(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1)
529                 
530                  IF( pt(ji,jj,jl) < -epsi20 ) THEN
531                     WRITE(numout,*) 'T<0 ',pt(ji,jj,jl)
532                  ENDIF
533                 
534                  IF( pt(ji,jj,jl) < 0._wp .AND. pt(ji,jj,jl) >= -epsi20 )   pt(ji,jj,jl) = 0._wp
535                 
536                  !!               IF( ji==26 .AND. jj==86) THEN
537                  !!                  WRITE(numout,*) 'zt high order',pt(ji,jj)
538                  !!               ENDIF
539               END DO
540            END DO
541         END DO
542         CALL lbc_lnk( 'icedyn_adv_umx', pt, 'T',  1. )
543      ENDIF
544   
545      ! Rachid trick
546      ! ------------
547      IF( ll_clem ) THEN
548         IF( pamsk == 0. ) THEN
549            DO jl = 1, jpl
550               DO jj = 1, jpjm1
551                  DO ji = 1, fs_jpim1
552                     IF( ABS( puc(ji,jj,jl) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN
553                        zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj)
554                        zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj)
555                     ELSE
556                        zfu_ho (ji,jj,jl) = 0._wp
557                        zfu_ups(ji,jj,jl) = 0._wp
558                     ENDIF
559                     !
560                     IF( ABS( pvc(ji,jj,jl) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN
561                        zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj)
562                        zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj)
563                     ELSE
564                        zfv_ho (ji,jj,jl) = 0._wp 
565                        zfv_ups(ji,jj,jl) = 0._wp 
566                     ENDIF
567                  END DO
568               END DO
569            END DO
570         ENDIF
571      ENDIF
572
573      IF( ll_zeroup5 ) THEN
574         DO jl = 1, jpl
575            DO jj = 2, jpjm1
576               DO ji = 2, fs_jpim1   ! vector opt.
577                  zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
578                     &                            - ( zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  ) * tmask(ji,jj,1)
579                  IF( zpt(ji,jj,jl) < 0. ) THEN
580                     zfu_ho(ji  ,jj,jl) = zfu_ups(ji  ,jj,jl)
581                     zfu_ho(ji-1,jj,jl) = zfu_ups(ji-1,jj,jl)
582                     zfv_ho(ji  ,jj,jl) = zfv_ups(ji  ,jj,jl)
583                     zfv_ho(ji,jj-1,jl) = zfv_ups(ji,jj-1,jl)
584                  ENDIF
585               END DO
586            END DO
587         END DO
588         CALL lbc_lnk_multi( 'icedyn_adv_umx', zfu_ho, 'U',  -1., zfv_ho, 'V',  -1. )
589      ENDIF
590
591      ! output high order fluxes u*a
592      ! ----------------------------
593      IF( PRESENT( pua_ho ) ) THEN
594         DO jl = 1, jpl
595            DO jj = 1, jpjm1
596               DO ji = 1, fs_jpim1
597                  pua_ho(ji,jj,jl) = zfu_ho(ji,jj,jl)
598                  pva_ho(ji,jj,jl) = zfv_ho(ji,jj,jl)
599               END DO
600            END DO
601         END DO
602      ENDIF
603
604
605      IF( .NOT.ll_thickness ) THEN
606         ! final trend with corrected fluxes
607         ! ------------------------------------
608         DO jl = 1, jpl
609            DO jj = 2, jpjm1
610               DO ji = fs_2, fs_jpim1 
611                  ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj) * pdt 
612                 
613                  ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra ) * tmask(ji,jj,1)
614                 
615                  !!               IF( ji==26 .AND. jj==86) THEN
616                  !!                  WRITE(numout,*) 'ztc high order',ptc(ji,jj)
617                  !!               ENDIF
618                 
619               END DO
620            END DO
621         END DO
622         CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T',  1. )
623      ENDIF
624     
625      !
626   END SUBROUTINE adv_umx
627
628   SUBROUTINE cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, &
629      &             pt_ups, pfu_ups, pfv_ups )
630      !!---------------------------------------------------------------------
631      !!                    ***  ROUTINE macho  ***
632      !!     
633      !! **  Purpose :   compute 
634      !!
635      !! **  Method  :   ... ???
636      !!                 TIM = transient interpolation Modeling
637      !!
638      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
639      !!----------------------------------------------------------------------
640      REAL(wp)                    , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
641      INTEGER                     , INTENT(in   ) ::   kn_limiter       ! limiter
642      INTEGER                     , INTENT(in   ) ::   jt               ! number of sub-iteration
643      INTEGER                     , INTENT(in   ) ::   kt               ! number of iteration
644      REAL(wp)                    , INTENT(in   ) ::   pdt              ! tracer time-step
645      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt               ! tracer fields
646      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
647      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components
648      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptc              ! tracer content at before time step
649      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes
650      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content
651      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes
652      !
653      INTEGER  ::   ji, jj, jl    ! dummy loop indices
654      LOGICAL  ::   ll_xy = .TRUE.
655      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zzt
656      !!----------------------------------------------------------------------
657      !
658      IF( .NOT.ll_xy ) THEN   !-- no alternate directions --!
659         !
660         DO jl = 1, jpl
661            DO jj = 1, jpjm1
662               DO ji = 1, fs_jpim1
663                  IF( ll_clem ) THEN
664                     pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
665                     pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
666                  ELSE
667                     pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
668                     pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
669                  ENDIF
670               END DO
671            END DO
672         END DO
673         IF    ( kn_limiter == 1 ) THEN
674            IF( ll_clem ) THEN
675               CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
676            ELSE
677               CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
678            ENDIF
679         ELSEIF( kn_limiter == 2 ) THEN
680            CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
681            CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
682         ELSEIF( kn_limiter == 3 ) THEN
683            CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
684            CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
685         ENDIF
686         !
687      ELSE                    !-- alternate directions --!
688         !
689         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
690            !
691            ! flux in x-direction
692            DO jl = 1, jpl
693               DO jj = 1, jpjm1
694                  DO ji = 1, fs_jpim1
695                     IF( ll_clem ) THEN
696                        pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
697                     ELSE
698                        pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
699                     ENDIF
700                  END DO
701               END DO
702            END DO
703            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
704            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
705
706            ! first guess of tracer content from u-flux
707            DO jl = 1, jpl
708               DO jj = 2, jpjm1
709                  DO ji = fs_2, fs_jpim1   ! vector opt.
710                     IF( ll_clem ) THEN
711                        zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)        &
712                           &            + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) &
713                           &         * tmask(ji,jj,1)
714                     ELSE                     
715                        zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
716                     ENDIF
717                  END DO
718               END DO
719            END DO
720            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
721
722            ! flux in y-direction
723            DO jl = 1, jpl
724               DO jj = 1, jpjm1
725                  DO ji = 1, fs_jpim1
726                     IF( ll_clem ) THEN
727                        pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji,jj+1,jl) )
728                     ELSE                     
729                        pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( zzt(ji,jj,jl) + zzt(ji,jj+1,jl) )
730                     ENDIF
731                  END DO
732               END DO
733            END DO
734            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
735            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
736
737         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
738            !
739            ! flux in y-direction
740            DO jl = 1, jpl
741               DO jj = 1, jpjm1
742                  DO ji = 1, fs_jpim1
743                     IF( ll_clem ) THEN
744                        pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
745                     ELSE                     
746                        pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
747                     ENDIF
748                  END DO
749               END DO
750            END DO
751            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
752            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
753            !
754            ! first guess of tracer content from v-flux
755            DO jl = 1, jpl
756               DO jj = 2, jpjm1
757                  DO ji = fs_2, fs_jpim1   ! vector opt.
758                     IF( ll_clem ) THEN
759                        zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &
760                           &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) &
761                           &         * tmask(ji,jj,1)
762                     ELSE
763                        zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
764                     ENDIF
765                  END DO
766               END DO
767            END DO
768            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
769            !
770            ! flux in x-direction
771            DO jl = 1, jpl
772               DO jj = 1, jpjm1
773                  DO ji = 1, fs_jpim1
774                     IF( ll_clem ) THEN
775                        pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji+1,jj,jl) )
776                     ELSE
777                        pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( zzt(ji,jj,jl) + zzt(ji+1,jj,jl) )
778                     ENDIF
779                  END DO
780               END DO
781            END DO
782            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
783            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
784
785         ENDIF
786         IF( ll_clem ) THEN
787            IF( kn_limiter == 1 )   CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
788         ELSE
789            IF( kn_limiter == 1 )   CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
790         ENDIF
791         
792      ENDIF
793   
794   END SUBROUTINE cen2
795
796   
797   SUBROUTINE macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, pt_u, pt_v, pfu_ho, pfv_ho, &
798      &              pt_ups, pfu_ups, pfv_ups )
799      !!---------------------------------------------------------------------
800      !!                    ***  ROUTINE macho  ***
801      !!     
802      !! **  Purpose :   compute 
803      !!
804      !! **  Method  :   ... ???
805      !!                 TIM = transient interpolation Modeling
806      !!
807      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
808      !!----------------------------------------------------------------------
809      REAL(wp)                    , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
810      INTEGER                     , INTENT(in   ) ::   kn_limiter       ! limiter
811      INTEGER                     , INTENT(in   ) ::   kn_umx           ! order of the scheme (1-5=UM or 20=CEN2)
812      INTEGER                     , INTENT(in   ) ::   jt               ! number of sub-iteration
813      INTEGER                     , INTENT(in   ) ::   kt               ! number of iteration
814      REAL(wp)                    , INTENT(in   ) ::   pdt              ! tracer time-step
815      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt               ! tracer fields
816      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
817      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components
818      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pubox, pvbox     ! upstream velocity
819      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptc              ! tracer content at before time step
820      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_u, pt_v       ! tracer at u- and v-points
821      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes
822      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content
823      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes
824      !
825      INTEGER  ::   ji, jj, jl    ! dummy loop indices
826      REAL(wp) ::   ztra
827      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zzt, zzfu_ho, zzfv_ho
828      !!----------------------------------------------------------------------
829      !
830      IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
831         !
832         !                                                        !--  ultimate interpolation of pt at u-point  --!
833         CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
834         !                                                        !--  limiter in x --!
835         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
836         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
837         !                                                        !--  advective form update in zzt  --!
838
839         IF( ll_1stguess_clem ) THEN
840
841            ! first guess of tracer content from u-flux
842            DO jl = 1, jpl
843               DO jj = 2, jpjm1
844                  DO ji = fs_2, fs_jpim1   ! vector opt.
845                     IF( ll_clem ) THEN
846                        IF( ll_gurvan ) THEN
847                           zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
848                        ELSE
849                           zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
850                              &            + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
851                              &         ) * tmask(ji,jj,1)
852                        ENDIF
853                     ELSE
854                        zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
855                     ENDIF
856                  END DO
857               END DO
858            END DO
859            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
860
861         ELSE
862
863            DO jl = 1, jpl
864               DO jj = 2, jpjm1
865                  DO ji = fs_2, fs_jpim1   ! vector opt.
866                     IF( ll_gurvan ) THEN
867                        zzt(ji,jj,jl) = pt(ji,jj,jl) - pubox(ji,jj   ) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj)  &
868                           &                         - pt   (ji,jj,jl) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj)
869                     ELSE
870                        zzt(ji,jj,jl) = pt(ji,jj,jl) - pubox(ji,jj   ) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj)  &
871                           &                         - pt   (ji,jj,jl) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk
872                     ENDIF
873                     zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1)
874                  END DO
875               END DO
876            END DO
877            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
878         ENDIF
879         !
880         !                                                        !--  ultimate interpolation of pt at v-point  --!
881         IF( ll_hoxy ) THEN
882            CALL ultimate_y( kn_umx, pdt, zzt, pv, pvc, pt_v, pfv_ho )
883         ELSE
884            CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
885         ENDIF
886         !                                                        !--  limiter in y --!
887         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
888         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
889         !         
890         !
891      ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
892         !
893         !                                                        !--  ultimate interpolation of pt at v-point  --!
894         CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
895         !                                                        !--  limiter in y --!
896         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
897         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
898         !                                                        !--  advective form update in zzt  --!
899         IF( ll_1stguess_clem ) THEN
900           
901            ! first guess of tracer content from v-flux
902            DO jl = 1, jpl
903               DO jj = 2, jpjm1
904                  DO ji = fs_2, fs_jpim1   ! vector opt.
905                     IF( ll_clem ) THEN
906                        IF( ll_gurvan ) THEN
907                           zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
908                        ELSE
909                           zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &
910                              &            + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
911                              &         ) * tmask(ji,jj,1)
912                        ENDIF
913                     ELSE
914                        zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) &
915                           &         * tmask(ji,jj,1)
916                     ENDIF
917                  END DO
918               END DO
919            END DO
920            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
921           
922         ELSE
923           
924            DO jl = 1, jpl
925               DO jj = 2, jpjm1
926                  DO ji = fs_2, fs_jpim1
927                     IF( ll_gurvan ) THEN
928                        zzt(ji,jj,jl) = pt(ji,jj,jl) - pvbox(ji,jj   ) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj)  &
929                           &                         - pt   (ji,jj,jl) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj)
930                     ELSE
931                        zzt(ji,jj,jl) = pt(ji,jj,jl) - pvbox(ji,jj   ) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj)  &
932                           &                         - pt   (ji,jj,jl) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk
933                     ENDIF
934                     zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1)
935                  END DO
936               END DO
937            END DO
938            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
939         ENDIF
940         !
941         !                                                        !--  ultimate interpolation of pt at u-point  --!
942         IF( ll_hoxy ) THEN
943            CALL ultimate_x( kn_umx, pdt, zzt, pu, puc, pt_u, pfu_ho )
944         ELSE
945            CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
946         ENDIF
947         !                                                        !--  limiter in x --!
948         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
949         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
950         !
951         !
952      ENDIF
953
954     
955      IF( kn_limiter == 1 ) THEN
956         IF( .NOT. ll_limiter_it2 ) THEN
957            IF( ll_clem ) THEN
958               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
959            ELSE
960               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
961            ENDIF
962         ELSE
963            zzfu_ho(:,:,:) = pfu_ho(:,:,:)
964            zzfv_ho(:,:,:) = pfv_ho(:,:,:)
965            ! 1st iteration of nonosc (limit the flux with the upstream solution)
966            IF( ll_clem ) THEN
967               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho )
968            ELSE
969               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho )
970            ENDIF
971            ! guess after content field with high order
972            DO jl = 1, jpl
973               DO jj = 2, jpjm1
974                  DO ji = fs_2, fs_jpim1
975                     ztra          = - ( zzfu_ho(ji,jj,jl) - zzfu_ho(ji-1,jj,jl) + zzfv_ho(ji,jj,jl) - zzfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj)
976                     zzt(ji,jj,jl) =   ( ptc(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1)
977                  END DO
978               END DO
979            END DO
980            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. )
981            ! 2nd iteration of nonosc (limit the flux with the limited high order solution)
982            IF( ll_clem ) THEN
983               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho )
984            ELSE
985               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho )
986            ENDIF
987         ENDIF
988      ENDIF
989      !
990   END SUBROUTINE macho
991
992
993   SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
994      !!---------------------------------------------------------------------
995      !!                    ***  ROUTINE ultimate_x  ***
996      !!     
997      !! **  Purpose :   compute 
998      !!
999      !! **  Method  :   ... ???
1000      !!                 TIM = transient interpolation Modeling
1001      !!
1002      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
1003      !!----------------------------------------------------------------------
1004      INTEGER                     , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
1005      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
1006      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pu        ! ice i-velocity component
1007      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   puc       ! ice i-velocity * A component
1008      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt        ! tracer fields
1009      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_u      ! tracer at u-point
1010      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho    ! high order flux
1011      !
1012      INTEGER  ::   ji, jj, jl             ! dummy loop indices
1013      REAL(wp) ::   zcu, zdx2, zdx4    !   -      -
1014      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztu1, ztu2, ztu3, ztu4
1015      !!----------------------------------------------------------------------
1016      !
1017      !                                                     !--  Laplacian in i-direction  --!
1018      DO jl = 1, jpl
1019         DO jj = 2, jpjm1         ! First derivative (gradient)
1020            DO ji = 1, fs_jpim1
1021               ztu1(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
1022            END DO
1023            !                     ! Second derivative (Laplacian)
1024            DO ji = fs_2, fs_jpim1
1025               ztu2(ji,jj,jl) = ( ztu1(ji,jj,jl) - ztu1(ji-1,jj,jl) ) * r1_e1t(ji,jj)
1026            END DO
1027         END DO
1028      END DO
1029      CALL lbc_lnk( 'icedyn_adv_umx', ztu2, 'T', 1. )
1030      !
1031      !                                                     !--  BiLaplacian in i-direction  --!
1032      DO jl = 1, jpl
1033         DO jj = 2, jpjm1         ! Third derivative
1034            DO ji = 1, fs_jpim1
1035               ztu3(ji,jj,jl) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
1036            END DO
1037            !                     ! Fourth derivative
1038            DO ji = fs_2, fs_jpim1
1039               ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj)
1040            END DO
1041         END DO
1042      END DO
1043      CALL lbc_lnk( 'icedyn_adv_umx', ztu4, 'T', 1. )
1044      !
1045      !
1046      SELECT CASE (kn_umx )
1047      !
1048      CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21)
1049         !       
1050         DO jl = 1, jpl
1051            DO jj = 1, jpjm1
1052               DO ji = 1, fs_jpim1   ! vector opt.
1053                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                           pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
1054                     &                                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
1055               END DO
1056            END DO
1057         END DO
1058         !
1059      CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23)
1060         !
1061         DO jl = 1, jpl
1062            DO jj = 1, jpjm1
1063               DO ji = 1, fs_jpim1   ! vector opt.
1064                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
1065                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
1066                     &                                               -              zcu   * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
1067               END DO
1068            END DO
1069         END DO
1070         
1071      CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24)
1072         !
1073         DO jl = 1, jpl
1074            DO jj = 1, jpjm1
1075               DO ji = 1, fs_jpim1   ! vector opt.
1076                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
1077                  zdx2 = e1u(ji,jj) * e1u(ji,jj)
1078!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
1079                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        &
1080                     &                                               -              zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   &
1081                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                         ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        &
1082                     &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )   )
1083               END DO
1084            END DO
1085         END DO
1086         !
1087      CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27)
1088         !
1089         DO jl = 1, jpl
1090            DO jj = 1, jpjm1
1091               DO ji = 1, fs_jpim1   ! vector opt.
1092                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
1093                  zdx2 = e1u(ji,jj) * e1u(ji,jj)
1094!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
1095                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        &
1096                     &                                               -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   &
1097                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                   ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        &
1098                     &                                               - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )   )
1099               END DO
1100            END DO
1101         END DO
1102         !
1103      CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29)
1104         !
1105         DO jl = 1, jpl
1106            DO jj = 1, jpjm1
1107               DO ji = 1, fs_jpim1   ! vector opt.
1108                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
1109                  zdx2 = e1u(ji,jj) * e1u(ji,jj)
1110                  !!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
1111                  zdx4 = zdx2 * zdx2
1112                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (            (                   pt  (ji+1,jj,jl) + pt  (ji,jj,jl)       &
1113                     &                                                     -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) )   &
1114                     &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *     (                   ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)       &
1115                     &                                                     - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) )   &
1116                     &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl)       &
1117                     &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) )
1118               END DO
1119            END DO
1120         END DO
1121         !
1122      END SELECT
1123      !                                                     !-- High order flux in i-direction  --!
1124      IF( ll_neg ) THEN
1125         DO jl = 1, jpl
1126            DO jj = 1, jpjm1
1127               DO ji = 1, fs_jpim1
1128                  IF( pt_u(ji,jj,jl) < 0._wp ) THEN
1129                     pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                           pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
1130                        &                                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
1131                  ENDIF
1132               END DO
1133            END DO
1134         END DO
1135      ENDIF
1136
1137      DO jl = 1, jpl
1138         DO jj = 1, jpjm1
1139            DO ji = 1, fs_jpim1   ! vector opt.
1140               IF( ll_clem ) THEN
1141                  pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl)
1142               ELSE
1143                  pfu_ho(ji,jj,jl) = puc(ji,jj,jl) * pt_u(ji,jj,jl)
1144               ENDIF
1145            END DO
1146         END DO
1147      END DO
1148      !
1149   END SUBROUTINE ultimate_x
1150   
1151 
1152   SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
1153      !!---------------------------------------------------------------------
1154      !!                    ***  ROUTINE ultimate_y  ***
1155      !!     
1156      !! **  Purpose :   compute 
1157      !!
1158      !! **  Method  :   ... ???
1159      !!                 TIM = transient interpolation Modeling
1160      !!
1161      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
1162      !!----------------------------------------------------------------------
1163      INTEGER                     , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
1164      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
1165      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pv        ! ice j-velocity component
1166      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pvc       ! ice j-velocity*A component
1167      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt        ! tracer fields
1168      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_v      ! tracer at v-point
1169      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfv_ho    ! high order flux
1170      !
1171      INTEGER  ::   ji, jj, jl       ! dummy loop indices
1172      REAL(wp) ::   zcv, zdy2, zdy4    !   -      -
1173      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztv1, ztv2, ztv3, ztv4
1174      !!----------------------------------------------------------------------
1175      !
1176      !                                                     !--  Laplacian in j-direction  --!
1177      DO jl = 1, jpl
1178         DO jj = 1, jpjm1         ! First derivative (gradient)
1179            DO ji = fs_2, fs_jpim1
1180               ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
1181            END DO
1182         END DO
1183         DO jj = 2, jpjm1         ! Second derivative (Laplacian)
1184            DO ji = fs_2, fs_jpim1
1185               ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj)
1186            END DO
1187         END DO
1188      END DO
1189      CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1. )
1190      !
1191      !                                                     !--  BiLaplacian in j-direction  --!
1192      DO jl = 1, jpl
1193         DO jj = 1, jpjm1         ! First derivative
1194            DO ji = fs_2, fs_jpim1
1195               ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
1196            END DO
1197         END DO
1198         DO jj = 2, jpjm1         ! Second derivative
1199            DO ji = fs_2, fs_jpim1
1200               ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj)
1201            END DO
1202         END DO
1203      END DO
1204      CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1. )
1205      !
1206      !
1207      SELECT CASE (kn_umx )
1208         !
1209      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21)
1210         DO jl = 1, jpl
1211            DO jj = 1, jpjm1
1212               DO ji = 1, fs_jpim1
1213                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                          ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  &
1214                     &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
1215               END DO
1216            END DO
1217         END DO
1218         !
1219      CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23)
1220         DO jl = 1, jpl
1221            DO jj = 1, jpjm1
1222               DO ji = 1, fs_jpim1
1223                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1224                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (     ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  &
1225                     &                                     - zcv * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
1226               END DO
1227            END DO
1228         END DO
1229         CALL lbc_lnk( 'icedyn_adv_umx', pt_v, 'V',  1. )
1230         !
1231      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24)
1232         DO jl = 1, jpl
1233            DO jj = 1, jpjm1
1234               DO ji = 1, fs_jpim1
1235                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1236                  zdy2 = e2v(ji,jj) * e2v(ji,jj)
1237!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
1238                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)       &
1239                     &                                     -                        zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) )   &
1240                     &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)       &
1241                     &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
1242               END DO
1243            END DO
1244         END DO
1245         !
1246      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27)
1247         DO jl = 1, jpl
1248            DO jj = 1, jpjm1
1249               DO ji = 1, fs_jpim1
1250                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1251                  zdy2 = e2v(ji,jj) * e2v(ji,jj)
1252!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
1253                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                        ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)       &
1254                     &                                               -          zcv * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) )   &
1255                     &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                   ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)       &
1256                     &                                               - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
1257               END DO
1258            END DO
1259         END DO
1260         !
1261      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29)
1262         DO jl = 1, jpl
1263            DO jj = 1, jpjm1
1264               DO ji = 1, fs_jpim1
1265                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1266                  zdy2 = e2v(ji,jj) * e2v(ji,jj)
1267!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
1268                  zdy4 = zdy2 * zdy2
1269                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)      &
1270                     &                                                     -          zcv * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) )  &
1271                     &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) *     (                   ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)      &
1272                     &                                                     - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) )  &
1273                     &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl)      &
1274                     &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) )
1275               END DO
1276            END DO
1277         END DO
1278         !
1279      END SELECT
1280      !                                                     !-- High order flux in j-direction  --!
1281      IF( ll_neg ) THEN
1282         DO jl = 1, jpl
1283            DO jj = 1, jpjm1
1284               DO ji = 1, fs_jpim1
1285                  IF( pt_v(ji,jj,jl) < 0._wp ) THEN
1286                     pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                          ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  &
1287                        &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
1288                  ENDIF
1289               END DO
1290            END DO
1291         END DO
1292      ENDIF
1293
1294      DO jl = 1, jpl
1295         DO jj = 1, jpjm1
1296            DO ji = 1, fs_jpim1   ! vector opt.
1297               IF( ll_clem ) THEN
1298                  pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl)
1299               ELSE
1300                  pfv_ho(ji,jj,jl) = pvc(ji,jj,jl) * pt_v(ji,jj,jl)
1301               ENDIF
1302            END DO
1303         END DO
1304      END DO
1305      !
1306   END SUBROUTINE ultimate_y
1307     
1308
1309   SUBROUTINE nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho )
1310      !!---------------------------------------------------------------------
1311      !!                    ***  ROUTINE nonosc  ***
1312      !!     
1313       !! **  Purpose :   compute monotonic tracer fluxes from the upstream
1314      !!       scheme and the before field by a nonoscillatory algorithm
1315      !!
1316      !! **  Method  :   ... ???
1317      !!       warning : pt and pt_low must be masked, but the boundaries
1318      !!       conditions on the fluxes are not necessary zalezak (1979)
1319      !!       drange (1995) multi-dimensional forward-in-time and upstream-
1320      !!       in-space based differencing for fluid
1321      !!----------------------------------------------------------------------
1322      REAL(wp)                     , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
1323      REAL(wp)                     , INTENT(in   ) ::   pdt              ! tracer time-step
1324      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pu               ! ice i-velocity => u*e2
1325      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   puc              ! ice i-velocity *A => u*e2*a
1326      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pv               ! ice j-velocity => v*e1
1327      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pvc              ! ice j-velocity *A => v*e1*a
1328      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   ptc, pt, pt_low  ! before field & upstream guess of after field
1329      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   pfv_low, pfu_low ! upstream flux
1330      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux
1331      !
1332      INTEGER  ::   ji, jj, jl    ! dummy loop indices
1333      REAL(wp) ::   zpos, zneg, zbig, zsml, z1_dt, zpos2, zneg2   ! local scalars
1334      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign, zcoef        !   -      -
1335      REAL(wp), DIMENSION(jpi,jpj    ) :: zbup, zbdo
1336      REAL(wp), DIMENSION(jpi,jpj,jpl) :: zbetup, zbetdo, zti_low, ztj_low, zzt
1337      !!----------------------------------------------------------------------
1338      zbig = 1.e+40_wp
1339      zsml = epsi20
1340     
1341      IF( ll_zeroup2 ) THEN
1342         DO jl = 1, jpl
1343            DO jj = 1, jpjm1
1344               DO ji = 1, fs_jpim1   ! vector opt.
1345                  IF( amaxu(ji,jj,jl) == 0._wp )   pfu_ho(ji,jj,jl) = 0._wp
1346                  IF( amaxv(ji,jj,jl) == 0._wp )   pfv_ho(ji,jj,jl) = 0._wp
1347               END DO
1348            END DO
1349         END DO
1350      ENDIF
1351
1352      IF( ll_zeroup4 ) THEN
1353         DO jl = 1, jpl
1354            DO jj = 1, jpjm1
1355               DO ji = 1, fs_jpim1   ! vector opt.
1356                  IF( pfu_low(ji,jj,jl) == 0._wp )   pfu_ho(ji,jj,jl) = 0._wp
1357                  IF( pfv_low(ji,jj,jl) == 0._wp )   pfv_ho(ji,jj,jl) = 0._wp
1358               END DO
1359            END DO
1360         END DO
1361      ENDIF
1362
1363
1364      IF( ll_zeroup1 ) THEN
1365         DO jl = 1, jpl
1366            DO jj = 2, jpjm1
1367               DO ji = fs_2, fs_jpim1
1368                  IF( ll_gurvan ) THEN
1369                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
1370                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
1371                  ELSE
1372                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
1373                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  &
1374                        &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1375                        &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1376                        &         ) * tmask(ji,jj,1)
1377                  ENDIF
1378                  IF( zzt(ji,jj,jl) < 0._wp ) THEN
1379                     pfu_ho(ji,jj,jl)   = pfu_low(ji,jj,jl)
1380                     pfv_ho(ji,jj,jl)   = pfv_low(ji,jj,jl)
1381                     WRITE(numout,*) '*** 1 negative high order zzt ***',ji,jj,zzt(ji,jj,jl)
1382                  ENDIF
1383                  !!               IF( ji==26 .AND. jj==86) THEN
1384                  !!                  WRITE(numout,*) 'zzt high order',zzt(ji,jj)
1385                  !!                  WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt
1386                  !!                  WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt
1387                  !!                  WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj,jl)) * r1_e1e2t(ji,jj) * pdt
1388                  !!                  WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1,jl)) * r1_e1e2t(ji,jj) * pdt
1389                  !!               ENDIF
1390                  IF( ll_gurvan ) THEN
1391                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
1392                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
1393                  ELSE
1394                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
1395                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  &
1396                        &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1397                        &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1398                        &         ) * tmask(ji,jj,1)
1399                  ENDIF
1400                  IF( zzt(ji,jj,jl) < 0._wp ) THEN
1401                     pfu_ho(ji-1,jj,jl) = pfu_low(ji-1,jj,jl)
1402                     pfv_ho(ji,jj-1,jl) = pfv_low(ji,jj-1,jl)
1403                     WRITE(numout,*) '*** 2 negative high order zzt ***',ji,jj,zzt(ji,jj,jl)
1404                  ENDIF
1405                  IF( ll_gurvan ) THEN
1406                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
1407                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
1408                  ELSE
1409                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
1410                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  &
1411                        &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1412                        &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1413                        &         ) * tmask(ji,jj,1)
1414                  ENDIF
1415                  IF( zzt(ji,jj,jl) < 0._wp ) THEN
1416                     WRITE(numout,*) '*** 3 negative high order zzt ***',ji,jj,zzt(ji,jj,jl)
1417                  ENDIF
1418               END DO
1419            END DO
1420         END DO
1421         CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. )
1422      ENDIF
1423
1424     
1425      ! antidiffusive flux : high order minus low order
1426      ! --------------------------------------------------
1427      DO jl = 1, jpl
1428         DO jj = 1, jpjm1
1429            DO ji = 1, fs_jpim1   ! vector opt.
1430               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_low(ji,jj,jl)
1431               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_low(ji,jj,jl)
1432            END DO
1433         END DO
1434      END DO
1435
1436      ! extreme case where pfu_ho has to be zero
1437      ! ----------------------------------------
1438      !                                    pfu_ho
1439      !                           *         --->
1440      !                        |      |  *   |        |
1441      !                        |      |      |    *   |   
1442      !                        |      |      |        |    *
1443      !            t_low :       i-1     i       i+1       i+2   
1444      IF( ll_prelimiter_zalesak ) THEN
1445         
1446         DO jl = 1, jpl
1447            DO jj = 2, jpjm1
1448               DO ji = fs_2, fs_jpim1 
1449                  zti_low(ji,jj,jl)= pt_low(ji+1,jj  ,jl)
1450                  ztj_low(ji,jj,jl)= pt_low(ji  ,jj+1,jl)
1451               END DO
1452            END DO
1453         END DO
1454         CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_low, 'T', 1., ztj_low, 'T', 1. )
1455
1456         !! this does not work ??
1457         !!            DO jj = 2, jpjm1
1458         !!               DO ji = fs_2, fs_jpim1
1459         !!                  IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji+1,jj  ) - pt_low (ji  ,jj) ) .AND.     &
1460         !!               &      SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj+1) - pt_low (ji  ,jj) )           &
1461         !!               &    ) THEN
1462         !!                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_low(ji+1,jj ) - zti_low(ji  ,jj) ) .AND.   &
1463         !!               &         SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_low(ji,jj+1 ) - ztj_low(ji  ,jj) )         &
1464         !!               &       ) THEN
1465         !!                        pfu_ho(ji,jj) = 0.  ;   pfv_ho(ji,jj) = 0.
1466         !!                     ENDIF
1467         !!                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj) - pt_low (ji-1,jj  ) ) .AND.  &
1468         !!               &         SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj) - pt_low (ji  ,jj-1) )        &
1469         !!               &       )  THEN
1470         !!                        pfu_ho(ji,jj) = 0.  ;   pfv_ho(ji,jj) = 0.
1471         !!                     ENDIF
1472         !!                  ENDIF
1473         !!                END DO
1474         !!            END DO           
1475
1476         DO jl = 1, jpl
1477            DO jj = 2, jpjm1
1478               DO ji = fs_2, fs_jpim1
1479                  IF ( pfu_ho(ji,jj,jl) * ( pt_low(ji+1,jj,jl) - pt_low(ji,jj,jl) ) <= 0. .AND.  &
1480                     & pfv_ho(ji,jj,jl) * ( pt_low(ji,jj+1,jl) - pt_low(ji,jj,jl) ) <= 0. ) THEN
1481                     !
1482                     IF(  pfu_ho(ji,jj,jl) * ( zti_low(ji+1,jj,jl) - zti_low(ji,jj,jl) ) <= 0 .AND.  &
1483                        & pfv_ho(ji,jj,jl) * ( ztj_low(ji,jj+1,jl) - ztj_low(ji,jj,jl) ) <= 0)  pfu_ho(ji,jj,jl)=0. ; pfv_ho(ji,jj,jl)=0.
1484                     !
1485                     IF(  pfu_ho(ji,jj,jl) * ( pt_low(ji  ,jj,jl) - pt_low(ji-1,jj,jl) ) <= 0 .AND.  &
1486                        & pfv_ho(ji,jj,jl) * ( pt_low(ji  ,jj,jl) - pt_low(ji,jj-1,jl) ) <= 0)  pfu_ho(ji,jj,jl)=0. ; pfv_ho(ji,jj,jl)=0.
1487                     !
1488                  ENDIF
1489               END DO
1490            END DO
1491         END DO
1492         CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond.
1493
1494      ELSEIF( ll_prelimiter_devore ) THEN
1495         DO jl = 1, jpl
1496            DO jj = 2, jpjm1
1497               DO ji = fs_2, fs_jpim1 
1498                  zti_low(ji,jj,jl)= pt_low(ji+1,jj  ,jl)
1499                  ztj_low(ji,jj,jl)= pt_low(ji  ,jj+1,jl)
1500               END DO
1501            END DO
1502         END DO
1503         CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_low, 'T', 1., ztj_low, 'T', 1. )
1504
1505         z1_dt = 1._wp / pdt
1506         DO jl = 1, jpl
1507            DO jj = 2, jpjm1
1508               DO ji = fs_2, fs_jpim1
1509                  zsign = SIGN( 1., pt_low(ji+1,jj,jl) - pt_low(ji,jj,jl) )
1510                  pfu_ho(ji,jj,jl) =  zsign * MAX( 0. , MIN( ABS(pfu_ho(ji,jj,jl)) , &
1511                     &                          zsign * ( pt_low (ji  ,jj,jl) - pt_low (ji-1,jj,jl) ) * e1e2t(ji  ,jj) * z1_dt , &
1512                     &                          zsign * ( zti_low(ji+1,jj,jl) - zti_low(ji  ,jj,jl) ) * e1e2t(ji+1,jj) * z1_dt ) )
1513
1514                  zsign = SIGN( 1., pt_low(ji,jj+1,jl) - pt_low(ji,jj,jl) )
1515                  pfv_ho(ji,jj,jl) =  zsign * MAX( 0. , MIN( ABS(pfv_ho(ji,jj,jl)) , &
1516                     &                          zsign * ( pt_low (ji,jj  ,jl) - pt_low (ji,jj-1,jl) ) * e1e2t(ji,jj  ) * z1_dt , &
1517                     &                          zsign * ( ztj_low(ji,jj+1,jl) - ztj_low(ji,jj  ,jl) ) * e1e2t(ji,jj+1) * z1_dt ) )
1518               END DO
1519            END DO
1520         END DO
1521         CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond.
1522
1523      ENDIF
1524
1525
1526      ! Search local extrema
1527      ! --------------------
1528      ! max/min of pt & pt_low with large negative/positive value (-/+zbig) outside ice cover
1529      z1_dt = 1._wp / pdt
1530      DO jl = 1, jpl
1531         
1532         DO jj = 1, jpj
1533            DO ji = 1, jpi
1534               IF    ( pt(ji,jj,jl) <= 0._wp .AND. pt_low(ji,jj,jl) <= 0._wp ) THEN
1535                  zbup(ji,jj) = -zbig
1536                  zbdo(ji,jj) =  zbig
1537               ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_low(ji,jj,jl) > 0._wp ) THEN
1538                  zbup(ji,jj) = pt_low(ji,jj,jl)
1539                  zbdo(ji,jj) = pt_low(ji,jj,jl)
1540               ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_low(ji,jj,jl) <= 0._wp ) THEN
1541                  zbup(ji,jj) = pt(ji,jj,jl)
1542                  zbdo(ji,jj) = pt(ji,jj,jl)
1543               ELSE
1544                  zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_low(ji,jj,jl) )
1545                  zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_low(ji,jj,jl) )
1546               ENDIF
1547            END DO
1548         END DO
1549
1550         DO jj = 2, jpjm1
1551            DO ji = fs_2, fs_jpim1   ! vector opt.
1552               !
1553               IF( .NOT. ll_9points ) THEN
1554                  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
1555                  zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1) )
1556                  !
1557               ELSE
1558                  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
1559                     &                     zbup(ji-1,jj-1), zbup(ji+1,jj+1), zbup(ji+1,jj-1), zbup(ji-1,jj+1)  )
1560                  zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1), &
1561                     &                     zbdo(ji-1,jj-1), zbdo(ji+1,jj+1), zbdo(ji+1,jj-1), zbdo(ji-1,jj+1)  )
1562               ENDIF
1563               !
1564               zpos = MAX( 0., pfu_ho(ji-1,jj,jl) ) - MIN( 0., pfu_ho(ji  ,jj,jl) ) &  ! positive/negative part of the flux
1565                  & + MAX( 0., pfv_ho(ji,jj-1,jl) ) - MIN( 0., pfv_ho(ji,jj  ,jl) )
1566               zneg = MAX( 0., pfu_ho(ji  ,jj,jl) ) - MIN( 0., pfu_ho(ji-1,jj,jl) ) &
1567                  & + MAX( 0., pfv_ho(ji,jj  ,jl) ) - MIN( 0., pfv_ho(ji,jj-1,jl) )
1568               !
1569               IF( ll_HgradU .AND. .NOT.ll_gurvan ) THEN
1570                  zneg2 = (   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) ) &
1571                     &    ) * ( 1. - pamsk )
1572                  zpos2 = ( - 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) ) &
1573                     &    ) * ( 1. - pamsk )
1574               ELSE
1575                  zneg2 = 0. ; zpos2 = 0.
1576               ENDIF
1577               !
1578               !                                  ! up & down beta terms
1579               IF( (zpos+zpos2) > 0. ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_low(ji,jj,jl) ) / (zpos+zpos2) * e1e2t(ji,jj) * z1_dt
1580               ELSE                         ; zbetup(ji,jj,jl) = 0. ! zbig
1581               ENDIF
1582               !
1583               IF( (zneg+zneg2) > 0. ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_low(ji,jj,jl) - zdo ) / (zneg+zneg2) * e1e2t(ji,jj) * z1_dt
1584               ELSE                         ; zbetdo(ji,jj,jl) = 0. ! zbig
1585               ENDIF
1586               !
1587               ! if all the points are outside ice cover
1588               IF( zup == -zbig )   zbetup(ji,jj,jl) = 0. ! zbig
1589               IF( zdo ==  zbig )   zbetdo(ji,jj,jl) = 0. ! zbig           
1590               !
1591
1592!!            IF( ji==26 .AND. jj==86) THEN
1593!               WRITE(numout,*) '-----------------'
1594!               WRITE(numout,*) 'zpos',zpos,zpos2
1595!               WRITE(numout,*) 'zneg',zneg,zneg2
1596!               WRITE(numout,*) 'puc/pu',ABS(puc(ji,jj))/MAX(epsi20, ABS(pu(ji,jj)))
1597!               WRITE(numout,*) 'pvc/pv',ABS(pvc(ji,jj))/MAX(epsi20, ABS(pv(ji,jj)))
1598!               WRITE(numout,*) 'pucm1/pu',ABS(puc(ji-1,jj))/MAX(epsi20, ABS(pu(ji-1,jj)))
1599!               WRITE(numout,*) 'pvcm1/pv',ABS(pvc(ji,jj-1))/MAX(epsi20, ABS(pv(ji,jj-1)))
1600!               WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)+pfu_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt
1601!               WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)+pfv_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt
1602!               WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)+pfu_low(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt
1603!               WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)+pfv_low(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt
1604!               WRITE(numout,*) 'pfu_low',pfu_low(ji,jj) * r1_e1e2t(ji,jj) * pdt
1605!               WRITE(numout,*) 'pfv_low',pfv_low(ji,jj) * r1_e1e2t(ji,jj) * pdt
1606!               WRITE(numout,*) 'pfu_lowm1',pfu_low(ji-1,jj) * r1_e1e2t(ji,jj) * pdt
1607!               WRITE(numout,*) 'pfv_lowm1',pfv_low(ji,jj-1) * r1_e1e2t(ji,jj) * pdt
1608!               
1609!               WRITE(numout,*) 'pt',pt(ji,jj)
1610!               WRITE(numout,*) 'ptim1',pt(ji-1,jj)
1611!               WRITE(numout,*) 'ptjm1',pt(ji,jj-1)
1612!               WRITE(numout,*) 'pt_low',pt_low(ji,jj)
1613!               WRITE(numout,*) 'zbetup',zbetup(ji,jj)
1614!               WRITE(numout,*) 'zbetdo',zbetdo(ji,jj)
1615!               WRITE(numout,*) 'zup',zup
1616!               WRITE(numout,*) 'zdo',zdo
1617!            ENDIF
1618            !
1619            END DO
1620         END DO
1621      END DO
1622      CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
1623
1624     
1625      ! monotonic flux in the y direction
1626      ! ---------------------------------
1627      DO jl = 1, jpl
1628         DO jj = 1, jpjm1
1629            DO ji = 1, fs_jpim1   ! vector opt.
1630               zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) )
1631               zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) )
1632               zcu = 0.5  + SIGN( 0.5 , pfu_ho(ji,jj,jl) )
1633               !
1634               zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu )
1635
1636               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_low(ji,jj,jl)
1637
1638!!            IF( ji==26 .AND. jj==86) THEN
1639!!               WRITE(numout,*) 'coefU',zcoef
1640!!               WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt
1641!!               WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj,jl)) * r1_e1e2t(ji,jj) * pdt
1642!!            ENDIF
1643
1644            END DO
1645         END DO
1646
1647         DO jj = 1, jpjm1
1648            DO ji = 1, fs_jpim1   ! vector opt.
1649               zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) )
1650               zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) )
1651               zcv = 0.5  + SIGN( 0.5 , pfv_ho(ji,jj,jl) )
1652               !
1653               zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv )
1654
1655               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_low(ji,jj,jl)
1656
1657!!            IF( ji==26 .AND. jj==86) THEN
1658!!               WRITE(numout,*) 'coefV',zcoef
1659!!               WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt
1660!!               WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1,jl)) * r1_e1e2t(ji,jj) * pdt
1661!!            ENDIF
1662            END DO
1663         END DO
1664
1665         ! clem test
1666         DO jj = 2, jpjm1
1667            DO ji = 2, fs_jpim1   ! vector opt.
1668               IF( ll_gurvan ) THEN
1669                  zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
1670                     &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
1671               ELSE
1672                  zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  &
1673                     &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  &
1674                     &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1675                     &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1676                     &         ) * tmask(ji,jj,1)
1677               ENDIF
1678               IF( zzt(ji,jj,jl) < -epsi20 ) THEN
1679                  WRITE(numout,*) 'T<0 nonosc',zzt(ji,jj,jl)
1680               ENDIF
1681            END DO
1682         END DO
1683
1684      END DO
1685
1686      !
1687      !
1688   END SUBROUTINE nonosc_2d
1689
1690   SUBROUTINE limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
1691      !!---------------------------------------------------------------------
1692      !!                    ***  ROUTINE limiter_x  ***
1693      !!     
1694      !! **  Purpose :   compute flux limiter
1695      !!----------------------------------------------------------------------
1696      REAL(wp)                   , INTENT(in   )           ::   pdt          ! tracer time-step
1697      REAL(wp), DIMENSION (:,:  ), INTENT(in   )           ::   pu           ! ice i-velocity => u*e2
1698      REAL(wp), DIMENSION (:,:,:), INTENT(in   )           ::   puc          ! ice i-velocity *A => u*e2*a
1699      REAL(wp), DIMENSION (:,:,:), INTENT(in   )           ::   pt           ! ice tracer
1700      REAL(wp), DIMENSION (:,:,:), INTENT(inout)           ::   pfu_ho       ! high order flux
1701      REAL(wp), DIMENSION (:,:,:), INTENT(in   ), OPTIONAL ::   pfu_ups      ! upstream flux
1702      !
1703      REAL(wp) ::   Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr
1704      INTEGER  ::   ji, jj, jl    ! dummy loop indices
1705      REAL(wp), DIMENSION (jpi,jpj,jpl) ::   zslpx       ! tracer slopes
1706      !!----------------------------------------------------------------------
1707      !
1708      DO jl = 1, jpl
1709         DO jj = 2, jpjm1
1710            DO ji = fs_2, fs_jpim1   ! vector opt.
1711               zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1)
1712            END DO
1713         END DO
1714      END DO
1715      CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.)   ! lateral boundary cond.
1716     
1717      DO jl = 1, jpl
1718         DO jj = 2, jpjm1
1719            DO ji = fs_2, fs_jpim1   ! vector opt.
1720               uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj)
1721               
1722               Rjm = zslpx(ji-1,jj,jl)
1723               Rj  = zslpx(ji  ,jj,jl)
1724               Rjp = zslpx(ji+1,jj,jl)
1725
1726               IF( PRESENT(pfu_ups) ) THEN
1727
1728                  IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm
1729                  ELSE                        ;   Rr = Rjp
1730                  ENDIF
1731
1732                  zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)     
1733                  IF( Rj > 0. ) THEN
1734                     zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pu(ji,jj)),  &
1735                        &        MIN( 2. * Rr * 0.5 * ABS(pu(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
1736                  ELSE
1737                     zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)),  &
1738                        &        MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) )
1739                  ENDIF
1740                  pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter
1741
1742               ELSE
1743                  IF( Rj /= 0. ) THEN
1744                     IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
1745                     ELSE                        ;   Cr = Rjp / Rj
1746                     ENDIF
1747                  ELSE
1748                     Cr = 0.
1749                     !IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm * 1.e20
1750                     !ELSE                        ;   Cr = Rjp * 1.e20
1751                     !ENDIF
1752                  ENDIF
1753
1754                  ! -- superbee --
1755                  zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
1756                  ! -- van albada 2 --
1757                  !!zpsi = 2.*Cr / (Cr*Cr+1.)
1758
1759                  ! -- sweby (with beta=1) --
1760                  !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
1761                  ! -- van Leer --
1762                  !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
1763                  ! -- ospre --
1764                  !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
1765                  ! -- koren --
1766                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
1767                  ! -- charm --
1768                  !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
1769                  !ELSE                 ;   zpsi = 0.
1770                  !ENDIF
1771                  ! -- van albada 1 --
1772                  !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
1773                  ! -- smart --
1774                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
1775                  ! -- umist --
1776                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
1777
1778                  ! high order flux corrected by the limiter
1779                  pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5
1780
1781               ENDIF
1782            END DO
1783         END DO
1784      END DO
1785      CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.)   ! lateral boundary cond.
1786      !
1787   END SUBROUTINE limiter_x
1788
1789   SUBROUTINE limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
1790      !!---------------------------------------------------------------------
1791      !!                    ***  ROUTINE limiter_y  ***
1792      !!     
1793      !! **  Purpose :   compute flux limiter
1794      !!----------------------------------------------------------------------
1795      REAL(wp)                   , INTENT(in   )           ::   pdt          ! tracer time-step
1796      REAL(wp), DIMENSION (:,:  ), INTENT(in   )           ::   pv           ! ice i-velocity => u*e2
1797      REAL(wp), DIMENSION (:,:,:), INTENT(in   )           ::   pvc          ! ice i-velocity *A => u*e2*a
1798      REAL(wp), DIMENSION (:,:,:), INTENT(in   )           ::   pt           ! ice tracer
1799      REAL(wp), DIMENSION (:,:,:), INTENT(inout)           ::   pfv_ho       ! high order flux
1800      REAL(wp), DIMENSION (:,:,:), INTENT(in   ), OPTIONAL ::   pfv_ups      ! upstream flux
1801      !
1802      REAL(wp) ::   Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr
1803      INTEGER  ::   ji, jj, jl    ! dummy loop indices
1804      REAL(wp), DIMENSION (jpi,jpj,jpl) ::   zslpy       ! tracer slopes
1805      !!----------------------------------------------------------------------
1806      !
1807      DO jl = 1, jpl
1808         DO jj = 2, jpjm1
1809            DO ji = fs_2, fs_jpim1   ! vector opt.
1810               zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1)
1811            END DO
1812         END DO
1813      END DO
1814      CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.)   ! lateral boundary cond.
1815
1816      DO jl = 1, jpl
1817         DO jj = 2, jpjm1
1818            DO ji = fs_2, fs_jpim1   ! vector opt.
1819               vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj)
1820
1821               Rjm = zslpy(ji,jj-1,jl)
1822               Rj  = zslpy(ji,jj  ,jl)
1823               Rjp = zslpy(ji,jj+1,jl)
1824
1825               IF( PRESENT(pfv_ups) ) THEN
1826
1827                  IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm
1828                  ELSE                        ;   Rr = Rjp
1829                  ENDIF
1830
1831                  zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)     
1832                  IF( Rj > 0. ) THEN
1833                     zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pv(ji,jj)),  &
1834                        &        MIN( 2. * Rr * 0.5 * ABS(pv(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
1835                  ELSE
1836                     zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)),  &
1837                        &        MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) )
1838                  ENDIF
1839                  pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter
1840
1841               ELSE
1842
1843                  IF( Rj /= 0. ) THEN
1844                     IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj
1845                     ELSE                        ;   Cr = Rjp / Rj
1846                     ENDIF
1847                  ELSE
1848                     Cr = 0.
1849                     !IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm * 1.e20
1850                     !ELSE                        ;   Cr = Rjp * 1.e20
1851                     !ENDIF
1852                  ENDIF
1853
1854                  ! -- superbee --
1855                  zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) )
1856                  ! -- van albada 2 --
1857                  !!zpsi = 2.*Cr / (Cr*Cr+1.)
1858
1859                  ! -- sweby (with beta=1) --
1860                  !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
1861                  ! -- van Leer --
1862                  !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) )
1863                  ! -- ospre --
1864                  !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. )
1865                  ! -- koren --
1866                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) )
1867                  ! -- charm --
1868                  !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) )
1869                  !ELSE                 ;   zpsi = 0.
1870                  !ENDIF
1871                  ! -- van albada 1 --
1872                  !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1)
1873                  ! -- smart --
1874                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) )
1875                  ! -- umist --
1876                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) )
1877
1878                  ! high order flux corrected by the limiter
1879                  pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5
1880
1881               ENDIF
1882            END DO
1883         END DO
1884      END DO
1885      CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.)   ! lateral boundary cond.
1886      !
1887   END SUBROUTINE limiter_y
1888
1889#else
1890   !!----------------------------------------------------------------------
1891   !!   Default option           Dummy module         NO SI3 sea-ice model
1892   !!----------------------------------------------------------------------
1893#endif
1894
1895   !!======================================================================
1896END MODULE icedyn_adv_umx
Note: See TracBrowser for help on using the repository browser.