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

Last change on this file since 10413 was 10413, checked in by clem, 2 years ago

merge dev_r9947_SI3_advection with the trunk. All sette tests passed. There is probably a conservation issue with the new advection scheme that I should solve but the routine icedyn_adv_umx.F90 is going to change anyway in the next couple of weeks. At worst, the old routine can be plugged without harm

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