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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icedyn_adv_umx.F90 in NEMO/trunk/src/ICE – NEMO

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

Last change on this file since 10418 was 10418, checked in by clem, 5 years ago

cosmetic changes

  • Property svn:keywords set to Id
File size: 86.0 KB
Line 
1MODULE icedyn_adv_umx
2   !!==============================================================================
3   !!                       ***  MODULE  icedyn_adv_umx  ***
4   !! sea-ice : advection using the ULTIMATE-MACHO scheme
5   !!==============================================================================
6   !! History :  3.6  !  2014-11  (C. Rousset, G. Madec)  Original code
7   !!            4.0  !  2018     (many people)           SI3 [aka Sea Ice cube]
8   !!----------------------------------------------------------------------
9#if defined key_si3
10   !!----------------------------------------------------------------------
11   !!   'key_si3'                                       SI3 sea-ice model
12   !!----------------------------------------------------------------------
13   !!   ice_dyn_adv_umx   : update the tracer trend with the 3D advection trends using a TVD scheme
14   !!   ultimate_x(_y)    : compute a tracer value at velocity points using ULTIMATE scheme at various orders
15   !!   macho             : ???
16   !!   nonosc            : compute monotonic tracer fluxes by a non-oscillatory algorithm
17   !!----------------------------------------------------------------------
18   USE phycst         ! physical constant
19   USE dom_oce        ! ocean domain
20   USE sbc_oce , ONLY : nn_fsbc   ! update frequency of surface boundary condition
21   USE ice            ! sea-ice variables
22   USE icevar         ! sea-ice: operations
23   !
24   USE in_out_manager ! I/O manager
25   USE lib_mpp        ! MPP library
26   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
27   USE lbclnk         ! lateral boundary conditions (or mpp links)
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   ice_dyn_adv_umx   ! called by icedyn_adv.F90
33     
34   REAL(wp) ::   z1_6   = 1._wp /   6._wp   ! =1/6
35   REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120
36
37   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: amaxu, amaxv
38   
39   ! advect H all the way (and get V=H*A at the end)
40   LOGICAL :: ll_thickness = .FALSE.
41   
42   ! look for 9 points around in nonosc limiter 
43   LOGICAL :: ll_9points = .FALSE.  ! false=better h?
44
45   ! use HgradU in nonosc limiter 
46   LOGICAL :: ll_HgradU = .TRUE.   ! no effect?
47
48   ! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme
49   LOGICAL :: ll_neg = .TRUE.      ! keep TRUE
50   
51   ! limit the fluxes
52   LOGICAL :: ll_zeroup1 = .FALSE. ! false ok if Hbig otherwise needed for 2D sinon on a des valeurs de H trop fortes !!
53   LOGICAL :: ll_zeroup2 = .FALSE. ! false ok for 1D, 2D, 3D
54   LOGICAL :: ll_zeroup4 = .FALSE. ! false ok for 1D, 2D, 3D
55   LOGICAL :: ll_zeroup5 = .FALSE. ! false ok for 1D, 2D
56
57   ! fluxes that are limited are u*H, then (u*H)*(ua/u) is used for V (only for nonosc)
58   LOGICAL :: ll_clem   = .TRUE.   ! simpler than rachid and works
59
60   ! First advect H as H*=Hdiv(u), then use H* for H(n+1)=H(n)-div(uH*)
61   LOGICAL :: ll_gurvan = .FALSE.  ! must be false for 1D case !!
62
63   ! First guess as div(uH) (-true-) or Hdiv(u)+ugradH (-false-)
64   LOGICAL :: ll_1stguess_clem = .FALSE. ! better negative values but less good h
65
66   ! advect (or not) open water. If not, retrieve it from advection of A
67   LOGICAL :: ll_ADVopw = .FALSE.  !
68   
69   ! alternate directions for upstream
70   LOGICAL :: ll_upsxy = .TRUE.
71
72   ! alternate directions for high order
73   LOGICAL :: ll_hoxy = .TRUE.
74   
75   ! prelimiter: use it to avoid overshoot in H
76   LOGICAL :: ll_prelimiter_zalesak = .TRUE.  ! from: Zalesak(1979) eq. 14 => true is better for 1D but false is better in 3D (for h and negative values) => pb in x-y?
77   LOGICAL :: ll_prelimiter_devore  = .FALSE.  ! from: Devore eq. 11 (worth than zalesak)
78
79   ! iterate on the limiter (only nonosc)
80   LOGICAL :: ll_limiter_it2 = .FALSE.
81   
82
83   !! * Substitutions
84#  include "vectopt_loop_substitute.h90"
85   !!----------------------------------------------------------------------
86   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
87   !! $Id$
88   !! Software governed by the CeCILL licence     (./LICENSE)
89   !!----------------------------------------------------------------------
90CONTAINS
91
92   SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice,  &
93      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
94      !!----------------------------------------------------------------------
95      !!                  ***  ROUTINE ice_dyn_adv_umx  ***
96      !!
97      !! **  Purpose :   Compute the now trend due to total advection of
98      !!                 tracers and add it to the general trend of tracer equations
99      !!                 using an "Ultimate-Macho" scheme
100      !!
101      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
102      !!----------------------------------------------------------------------
103      INTEGER                     , INTENT(in   ) ::   kn_umx     ! order of the scheme (1-5=UM or 20=CEN2)
104      INTEGER                     , INTENT(in   ) ::   kt         ! time step
105      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity
106      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity
107      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
108      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
109      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
110      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psv_i      ! salt content
111      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
112      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
113      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
114      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
115      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
116      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
117      !
118      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
119      INTEGER  ::   icycle                  ! number of sub-timestep for the advection
120      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers
121      REAL(wp) ::   zcfl , zdt
122      REAL(wp), DIMENSION(jpi,jpj) ::   zudy, zvdx, zcu_box, zcv_box, zua_ho, zva_ho
123      REAL(wp), DIMENSION(jpi,jpj) ::   zhvar
124      REAL(wp), DIMENSION(jpi,jpj) ::   zati1, zati2, z1_ai, z1_aip
125      !!----------------------------------------------------------------------
126      !
127      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
128      !
129      !
130      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !       
131      zcfl  =            MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
132      zcfl  = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
133      IF( lk_mpp )   CALL mpp_max( zcfl )
134
135      IF( zcfl > 0.5 ) THEN   ;   icycle = 2 
136      ELSE                    ;   icycle = 1 
137      ENDIF
138     
139      zdt = rdt_ice / REAL(icycle)
140
141      ! --- transport --- !
142      zudy(:,:) = pu_ice(:,:) * e2u(:,:)
143      zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
144
145      ! --- define velocity for advection: u*grad(H) --- !
146      DO jj = 2, jpjm1
147         DO ji = fs_2, fs_jpim1
148            IF    ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp
149            ELSEIF( pu_ice(ji,jj)                   >  0._wp ) THEN   ;   zcu_box(ji,jj) = pu_ice(ji-1,jj)
150            ELSE                                                      ;   zcu_box(ji,jj) = pu_ice(ji  ,jj)
151            ENDIF
152
153            IF    ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp
154            ELSEIF( pv_ice(ji,jj)                   >  0._wp ) THEN   ;   zcv_box(ji,jj) = pv_ice(ji,jj-1)
155            ELSE                                                      ;   zcv_box(ji,jj) = pv_ice(ji,jj  )
156            ENDIF
157         END DO
158      END DO
159
160      IF( ll_zeroup2 ) THEN
161         IF(.NOT. ALLOCATED(amaxu))       ALLOCATE(amaxu (jpi,jpj))
162         IF(.NOT. ALLOCATED(amaxv))       ALLOCATE(amaxv (jpi,jpj))
163      ENDIF
164      !---------------!
165      !== advection ==!
166      !---------------!
167      DO jt = 1, icycle
168
169         IF( ll_ADVopw ) THEN
170            zamsk = 1._wp
171            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) )   ! Open water area
172            zamsk = 0._wp
173         ELSE
174            zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
175         ENDIF
176         
177         DO jl = 1, jpl
178            !
179            WHERE( pa_i(:,:,jl) >= epsi20 )   ;   z1_ai(:,:) = 1._wp / pa_i(:,:,jl)
180            ELSEWHERE                         ;   z1_ai(:,:) = 0.
181            END WHERE
182            !
183            WHERE( pa_ip(:,:,jl) >= epsi20 )  ;   z1_aip(:,:) = 1._wp / pa_ip(:,:,jl)
184            ELSEWHERE                         ;   z1_aip(:,:) = 0.
185            END WHERE
186            !
187            IF( ll_zeroup2 ) THEN
188               DO jj = 2, jpjm1
189                  DO ji = fs_2, fs_jpim1
190                     amaxu(ji,jj)=MAX( pa_i(ji,jj,jl), pa_i(ji,jj-1,jl), pa_i(ji,jj+1,jl), &
191                        &                              pa_i(ji+1,jj,jl), pa_i(ji+1,jj-1,jl), pa_i(ji+1,jj+1,jl) )
192                     amaxv(ji,jj)=MAX( pa_i(ji,jj,jl), pa_i(ji-1,jj,jl), pa_i(ji+1,jj,jl), &
193                        &                              pa_i(ji,jj+1,jl), pa_i(ji-1,jj+1,jl), pa_i(ji+1,jj+1,jl) )
194                 END DO
195               END DO
196               CALL lbc_lnk_multi(amaxu, 'T', 1., amaxv, 'T', 1.)
197            ENDIF
198            !
199            zamsk = 1._wp
200            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl), pa_i(:,:,jl), &        ! Ice area
201               &          zua_ho, zva_ho )
202            zamsk = 0._wp
203            !
204            IF( ll_thickness ) THEN
205               zua_ho(:,:) = zudy(:,:)
206               zva_ho(:,:) = zvdx(:,:)
207            ENDIF
208            !
209            zhvar(:,:) = pv_i(:,:,jl) * z1_ai(:,:)
210            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
211            IF( ll_thickness )   pv_i(:,:,jl) = zhvar(:,:) * pa_i(:,:,jl)
212            !
213            zhvar(:,:) = pv_s(:,:,jl) * z1_ai(:,:)
214            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
215            IF( ll_thickness )   pv_s(:,:,jl) = zhvar(:,:) * pa_i(:,:,jl)
216            !
217            zhvar(:,:) = psv_i(:,:,jl) * z1_ai(:,:)
218            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
219            !
220            zhvar(:,:) = poa_i(:,:,jl) * z1_ai(:,:)
221            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
222            !
223            DO jk = 1, nlay_i
224               zhvar(:,:) = pe_i(:,:,jk,jl) * z1_ai(:,:)
225               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
226            END DO
227            !
228            DO jk = 1, nlay_s
229               zhvar(:,:) = pe_s(:,:,jk,jl) * z1_ai(:,:)
230               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
231            END DO
232            !
233            IF ( ln_pnd_H12 ) THEN
234               !
235               zamsk = 1._wp
236               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl), pa_ip(:,:,jl), &   ! mp fraction
237                  &          zua_ho, zva_ho )
238               zamsk = 0._wp
239
240               zhvar(:,:) = pv_ip(:,:,jl) * z1_ai(:,:)
241               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
242            ENDIF
243            !
244            !
245         END DO
246         !
247         IF( .NOT. ll_ADVopw ) THEN
248            zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
249            DO jj = 2, jpjm1
250               DO ji = fs_2, fs_jpim1
251                  pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                                                  ! Open water area
252                     &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt
253               END DO
254            END DO
255            CALL lbc_lnk( pato_i(:,:), 'T',  1. )
256         ENDIF
257         !
258      END DO
259      !
260   END SUBROUTINE ice_dyn_adv_umx
261
262   
263   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho )
264      !!----------------------------------------------------------------------
265      !!                  ***  ROUTINE adv_umx  ***
266      !!
267      !! **  Purpose :   Compute the now trend due to total advection of
268      !!       tracers and add it to the general trend of tracer equations
269      !!
270      !! **  Method  :   TVD scheme, i.e. 2nd order centered scheme with
271      !!       corrected flux (monotonic correction)
272      !!       note: - this advection scheme needs a leap-frog time scheme
273      !!
274      !! ** Action : - pt  the after advective tracer
275      !!----------------------------------------------------------------------
276      REAL(wp)                    , INTENT(in   )           ::   pamsk          ! advection of concentration (1) or other tracers (0)
277      INTEGER                     , INTENT(in   )           ::   kn_umx         ! order of the scheme (1-5=UM or 20=CEN2)
278      INTEGER                     , INTENT(in   )           ::   jt             ! number of sub-iteration
279      INTEGER                     , INTENT(in   )           ::   kt             ! number of iteration
280      REAL(wp)                    , INTENT(in   )           ::   pdt            ! tracer time-step
281      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pu   , pv      ! 2 ice velocity components => u*e2
282      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   puc  , pvc     ! 2 ice velocity components => u*e2 or u*a*e2u
283      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pubox, pvbox   ! upstream velocity
284      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   pt             ! tracer field
285      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   ptc            ! tracer content field
286      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho ! high order u*a fluxes
287      !
288      INTEGER  ::   ji, jj           ! dummy loop indices 
289      REAL(wp) ::   ztra             ! local scalar
290      INTEGER  ::   kn_limiter = 1   ! 1=nonosc ; 2=superbee ; 3=h3(rachid)
291      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ho , zfv_ho , zt_u, zt_v, zpt
292      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfv_ups, zt_ups   ! only for nonosc
293      !!----------------------------------------------------------------------
294      !
295      !  upstream (_ups) advection with initial mass fluxes
296      ! ---------------------------------------------------
297
298      IF( ll_gurvan .AND. pamsk==0. ) THEN
299         DO jj = 2, jpjm1
300            DO ji = fs_2, fs_jpim1
301               pt(ji,jj) = ( pt (ji,jj) + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) &
302                  &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
303            END DO
304         END DO
305         CALL lbc_lnk( pt, 'T', 1. )
306      ENDIF
307
308     
309      IF( .NOT. ll_upsxy ) THEN
310
311         ! fluxes in both x-y directions
312         DO jj = 1, jpjm1
313            DO ji = 1, fs_jpim1
314               IF( ll_clem ) THEN
315                  zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj)
316                  zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1)
317               ELSE
318                  zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj)
319                  zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1)
320               ENDIF
321            END DO
322         END DO
323
324      ELSE
325         !
326         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
327            ! flux in x-direction
328            DO jj = 1, jpjm1
329               DO ji = 1, fs_jpim1
330                  IF( ll_clem ) THEN
331                     zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj)
332                  ELSE
333                     zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj)
334                  ENDIF
335               END DO
336            END DO
337           
338            ! first guess of tracer content from u-flux
339            DO jj = 2, jpjm1
340               DO ji = fs_2, fs_jpim1   ! vector opt.
341                  IF( ll_clem ) THEN
342                     IF( ll_gurvan ) THEN
343                        zpt(ji,jj) = ( pt(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
344                     ELSE
345                        zpt(ji,jj) = ( pt(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
346                           &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
347                           &         ) * tmask(ji,jj,1)
348                     ENDIF
349                  ELSE
350                     zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) &
351                        &         * tmask(ji,jj,1)
352                  ENDIF
353!!                  IF( ji==26 .AND. jj==86) THEN
354!!                     WRITE(numout,*) '************************'
355!!                     WRITE(numout,*) 'zpt upstream',zpt(ji,jj)
356!!                  ENDIF
357               END DO
358            END DO
359            CALL lbc_lnk( zpt, 'T', 1. )
360            !
361            ! flux in y-direction
362            DO jj = 1, jpjm1
363               DO ji = 1, fs_jpim1
364                  IF( ll_clem ) THEN
365                     zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1)
366                  ELSE
367                     zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * zpt(ji,jj+1)
368                  ENDIF
369               END DO
370            END DO
371
372         !
373         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
374            ! flux in y-direction
375            DO jj = 1, jpjm1
376               DO ji = 1, fs_jpim1
377                  IF( ll_clem ) THEN
378                     zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1)
379                  ELSE
380                     zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1)
381                  ENDIF
382               END DO
383            END DO
384
385            ! first guess of tracer content from v-flux
386            DO jj = 2, jpjm1
387               DO ji = fs_2, fs_jpim1   ! vector opt.
388                  IF( ll_clem ) THEN
389                     IF( ll_gurvan ) THEN
390                        zpt(ji,jj) = ( pt(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
391                     ELSE
392                        zpt(ji,jj) = ( pt(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) &
393                        &                        + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) &
394                        &         * tmask(ji,jj,1)
395                     ENDIF
396                  ELSE
397                     zpt(ji,jj) = ( ptc(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) &
398                        &         * tmask(ji,jj,1)
399                  ENDIF
400!!                  IF( ji==26 .AND. jj==86) THEN
401!!                     WRITE(numout,*) '************************'
402!!                     WRITE(numout,*) 'zpt upstream',zpt(ji,jj)
403!!                  ENDIF
404                END DO
405            END DO
406            CALL lbc_lnk( zpt, 'T', 1. )
407            !
408            ! flux in x-direction
409            DO jj = 1, jpjm1
410               DO ji = 1, fs_jpim1
411                  IF( ll_clem ) THEN
412                     zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj)
413                  ELSE
414                     zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * zpt(ji+1,jj)
415                  ENDIF
416               END DO
417            END DO
418            !
419         ENDIF
420         
421      ENDIF
422
423      IF( ll_clem .AND. kn_limiter /= 1 )  &
424         & CALL ctl_stop( 'STOP', 'icedyn_adv_umx: ll_clem incompatible with limiters other than nonosc' )
425
426      IF( ll_zeroup2 ) THEN
427         DO jj = 1, jpjm1
428            DO ji = 1, fs_jpim1   ! vector opt.
429               IF( amaxu(ji,jj) == 0._wp )   zfu_ups(ji,jj) = 0._wp
430               IF( amaxv(ji,jj) == 0._wp )   zfv_ups(ji,jj) = 0._wp
431            END DO
432         END DO
433      ENDIF
434
435      ! guess after content field with upstream scheme
436      DO jj = 2, jpjm1
437         DO ji = fs_2, fs_jpim1
438            ztra          = - (   zfu_ups(ji,jj) - zfu_ups(ji-1,jj  )   &
439               &                + zfv_ups(ji,jj) - zfv_ups(ji  ,jj-1) ) * r1_e1e2t(ji,jj)
440            IF( ll_clem ) THEN
441               IF( ll_gurvan ) THEN
442                  zt_ups(ji,jj) = ( pt (ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
443               ELSE
444                  zt_ups(ji,jj) = ( pt (ji,jj) + pdt * ztra + ( pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) )   &
445                     &                                      +   pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) ) &
446                     &                                        * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1)
447               ENDIF
448            ELSE
449               zt_ups(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
450            ENDIF
451!!            IF( ji==26 .AND. jj==86) THEN
452!!               WRITE(numout,*) '**************************'
453!!               WRITE(numout,*) 'zt upstream',zt_ups(ji,jj)
454!!            ENDIF
455         END DO
456      END DO
457      CALL lbc_lnk( zt_ups, 'T', 1. )
458
459      ! High order (_ho) fluxes
460      ! -----------------------
461      SELECT CASE( kn_umx )
462         !
463      CASE ( 20 )                          !== centered second order ==!
464         !
465         CALL cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho,  &
466            &       zt_ups, zfu_ups, zfv_ups )
467         !
468      CASE ( 1:5 )                         !== 1st to 5th order ULTIMATE-MACHO scheme ==!
469         !
470         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,  &
471            &        zt_ups, zfu_ups, zfv_ups )
472         !
473      END SELECT
474
475      IF( ll_thickness ) THEN
476         ! final trend with corrected fluxes
477         ! ------------------------------------
478         DO jj = 2, jpjm1
479            DO ji = fs_2, fs_jpim1
480               IF( ll_gurvan ) THEN
481                  ztra       = - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj) 
482               ELSE
483                  ztra       = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) )  & 
484                     &           + ( pt(ji,jj) * ( pu(ji,jj) - pu(ji-1,jj) ) * (1.-pamsk) ) &
485                     &           + ( pt(ji,jj) * ( pv(ji,jj) - pv(ji,jj-1) ) * (1.-pamsk) ) ) * r1_e1e2t(ji,jj)
486               ENDIF
487               pt(ji,jj) = ( pt(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
488
489               IF( pt(ji,jj) < -epsi20 ) THEN
490                  WRITE(numout,*) 'T<0 ',pt(ji,jj)
491               ENDIF
492               
493               IF( pt(ji,jj) < 0._wp .AND. pt(ji,jj) >= -epsi20 )   pt(ji,jj) = 0._wp
494               
495!!               IF( ji==26 .AND. jj==86) THEN
496!!                  WRITE(numout,*) 'zt high order',pt(ji,jj)
497!!               ENDIF
498            END DO
499         END DO
500         CALL lbc_lnk( pt, 'T',  1. )
501      ENDIF
502   
503      ! Rachid trick
504      ! ------------
505      IF( ll_clem ) THEN
506         IF( pamsk == 0. ) THEN
507            DO jj = 1, jpjm1
508               DO ji = 1, fs_jpim1
509                  IF( ABS( puc(ji,jj) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN
510                     zfu_ho (ji,jj) = zfu_ho (ji,jj) * puc(ji,jj) / pu(ji,jj)
511                     zfu_ups(ji,jj) = zfu_ups(ji,jj) * puc(ji,jj) / pu(ji,jj)
512                  ELSE
513                     zfu_ho (ji,jj) = 0._wp
514                     zfu_ups(ji,jj) = 0._wp
515                  ENDIF
516                  !
517                  IF( ABS( pvc(ji,jj) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN
518                     zfv_ho (ji,jj) = zfv_ho (ji,jj) * pvc(ji,jj) / pv(ji,jj)
519                     zfv_ups(ji,jj) = zfv_ups(ji,jj) * pvc(ji,jj) / pv(ji,jj)
520                  ELSE
521                     zfv_ho (ji,jj) = 0._wp 
522                     zfv_ups(ji,jj) = 0._wp 
523                  ENDIF
524               ENDDO
525            ENDDO
526         ENDIF
527      ENDIF
528
529      IF( ll_zeroup5 ) THEN
530         DO jj = 2, jpjm1
531            DO ji = 2, fs_jpim1   ! vector opt.
532               zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
533                  &                      - ( zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  ) * tmask(ji,jj,1)
534               IF( zpt(ji,jj) < 0. ) THEN
535                  zfu_ho(ji,jj)   = zfu_ups(ji,jj)
536                  zfu_ho(ji-1,jj) = zfu_ups(ji-1,jj)
537                  zfv_ho(ji,jj)   = zfv_ups(ji,jj)
538                  zfv_ho(ji,jj-1) = zfv_ups(ji,jj-1)
539               ENDIF
540            END DO
541         END DO
542         CALL lbc_lnk_multi( zfu_ho, 'U',  -1., zfv_ho, 'V',  -1. )
543      ENDIF
544
545      ! output high order fluxes u*a
546      ! ----------------------------
547      IF( PRESENT( pua_ho ) ) THEN
548         DO jj = 1, jpjm1
549            DO ji = 1, fs_jpim1
550               pua_ho(ji,jj) = zfu_ho(ji,jj)
551               pva_ho(ji,jj) = zfv_ho(ji,jj)
552            END DO
553         END DO
554      ENDIF
555
556
557      IF( .NOT.ll_thickness ) THEN
558         ! final trend with corrected fluxes
559         ! ------------------------------------
560         DO jj = 2, jpjm1
561            DO ji = fs_2, fs_jpim1 
562               ztra = - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj) * pdt 
563
564               ptc(ji,jj) = ( ptc(ji,jj) + ztra ) * tmask(ji,jj,1)
565                             
566!!               IF( ji==26 .AND. jj==86) THEN
567!!                  WRITE(numout,*) 'ztc high order',ptc(ji,jj)
568!!               ENDIF
569               
570            END DO
571         END DO
572         CALL lbc_lnk( ptc, 'T',  1. )
573      ENDIF
574     
575      !
576   END SUBROUTINE adv_umx
577
578   SUBROUTINE cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, &
579      &             pt_ups, pfu_ups, pfv_ups )
580      !!---------------------------------------------------------------------
581      !!                    ***  ROUTINE macho  ***
582      !!     
583      !! **  Purpose :   compute 
584      !!
585      !! **  Method  :   ... ???
586      !!                 TIM = transient interpolation Modeling
587      !!
588      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
589      !!----------------------------------------------------------------------
590      REAL(wp)                    , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
591      INTEGER                     , INTENT(in   ) ::   kn_limiter       ! limiter
592      INTEGER                     , INTENT(in   ) ::   jt               ! number of sub-iteration
593      INTEGER                     , INTENT(in   ) ::   kt               ! number of iteration
594      REAL(wp)                    , INTENT(in   ) ::   pdt              ! tracer time-step
595      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt               ! tracer fields
596      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
597      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components
598      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc              ! tracer content at before time step
599      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes
600      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content
601      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes
602      !
603      INTEGER  ::   ji, jj    ! dummy loop indices
604      LOGICAL  ::   ll_xy = .TRUE.
605      REAL(wp), DIMENSION(jpi,jpj) ::   zzt
606      !!----------------------------------------------------------------------
607      !
608      IF( .NOT.ll_xy ) THEN   !-- no alternate directions --!
609         !
610         DO jj = 1, jpjm1
611            DO ji = 1, fs_jpim1
612               IF( ll_clem ) THEN
613                  pfu_ho(ji,jj) = 0.5 * pu(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) )
614                  pfv_ho(ji,jj) = 0.5 * pv(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) )
615               ELSE
616                  pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) )
617                  pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) )
618               ENDIF
619            END DO
620         END DO
621         IF    ( kn_limiter == 1 ) THEN
622            IF( ll_clem ) THEN
623               CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
624            ELSE
625               CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
626            ENDIF
627         ELSEIF( kn_limiter == 2 ) THEN
628            CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
629            CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
630         ELSEIF( kn_limiter == 3 ) THEN
631            CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
632            CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
633         ENDIF
634         !
635      ELSE                    !-- alternate directions --!
636         !
637         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
638            !
639            ! flux in x-direction
640            DO jj = 1, jpjm1
641               DO ji = 1, fs_jpim1
642                  IF( ll_clem ) THEN
643                     pfu_ho(ji,jj) = 0.5 * pu(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) )
644                  ELSE
645                     pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) )
646                  ENDIF
647               END DO
648            END DO
649            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
650            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
651
652            ! first guess of tracer content from u-flux
653            DO jj = 2, jpjm1
654               DO ji = fs_2, fs_jpim1   ! vector opt.
655                  IF( ll_clem ) THEN
656                     zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)        &
657                           &                  + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) &
658                           &         * tmask(ji,jj,1)
659                  ELSE                     
660                     zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
661                  ENDIF
662               END DO
663            END DO
664            CALL lbc_lnk( zzt, 'T', 1. )
665
666            ! flux in y-direction
667            DO jj = 1, jpjm1
668               DO ji = 1, fs_jpim1
669                  IF( ll_clem ) THEN
670                     pfv_ho(ji,jj) = 0.5 * pv(ji,jj) * ( zzt(ji,jj) + zzt(ji,jj+1) )
671                  ELSE                     
672                     pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( zzt(ji,jj) + zzt(ji,jj+1) )
673                  ENDIF
674               END DO
675            END DO
676            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
677            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
678
679         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
680            !
681            ! flux in y-direction
682            DO jj = 1, jpjm1
683               DO ji = 1, fs_jpim1
684                  IF( ll_clem ) THEN
685                     pfv_ho(ji,jj) = 0.5 * pv(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) )
686                  ELSE                     
687                     pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) )
688                  ENDIF
689               END DO
690            END DO
691            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
692            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
693            !
694            ! first guess of tracer content from v-flux
695            DO jj = 2, jpjm1
696               DO ji = fs_2, fs_jpim1   ! vector opt.
697                  IF( ll_clem ) THEN
698                     zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) &
699                        &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) &
700                        &         * tmask(ji,jj,1)
701                  ELSE
702                     zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
703                  ENDIF
704               END DO
705            END DO
706            CALL lbc_lnk( zzt, 'T', 1. )
707            !
708            ! flux in x-direction
709            DO jj = 1, jpjm1
710               DO ji = 1, fs_jpim1
711                  IF( ll_clem ) THEN
712                     pfu_ho(ji,jj) = 0.5 * pu(ji,jj) * ( zzt(ji,jj) + zzt(ji+1,jj) )
713                  ELSE
714                     pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( zzt(ji,jj) + zzt(ji+1,jj) )
715                  ENDIF
716               END DO
717            END DO
718            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
719            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
720
721         ENDIF
722         IF( ll_clem ) THEN
723            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 )
724         ELSE
725            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 )
726         ENDIF
727         
728      ENDIF
729   
730   END SUBROUTINE cen2
731
732   
733   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, &
734      &              pt_ups, pfu_ups, pfv_ups )
735      !!---------------------------------------------------------------------
736      !!                    ***  ROUTINE macho  ***
737      !!     
738      !! **  Purpose :   compute 
739      !!
740      !! **  Method  :   ... ???
741      !!                 TIM = transient interpolation Modeling
742      !!
743      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
744      !!----------------------------------------------------------------------
745      REAL(wp)                    , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
746      INTEGER                     , INTENT(in   ) ::   kn_limiter       ! limiter
747      INTEGER                     , INTENT(in   ) ::   kn_umx           ! order of the scheme (1-5=UM or 20=CEN2)
748      INTEGER                     , INTENT(in   ) ::   jt               ! number of sub-iteration
749      INTEGER                     , INTENT(in   ) ::   kt               ! number of iteration
750      REAL(wp)                    , INTENT(in   ) ::   pdt              ! tracer time-step
751      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt               ! tracer fields
752      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
753      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components
754      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox     ! upstream velocity
755      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc              ! tracer content at before time step
756      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u, pt_v       ! tracer at u- and v-points
757      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes
758      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content
759      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes
760      !
761      INTEGER  ::   ji, jj    ! dummy loop indices
762      REAL(wp) ::   ztra
763      REAL(wp), DIMENSION(jpi,jpj) ::   zzt, zzfu_ho, zzfv_ho
764      !!----------------------------------------------------------------------
765      !
766      IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
767         !
768         !                                                        !--  ultimate interpolation of pt at u-point  --!
769         CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
770         !                                                        !--  limiter in x --!
771         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
772         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
773         !                                                        !--  advective form update in zzt  --!
774
775         IF( ll_1stguess_clem ) THEN
776
777            ! first guess of tracer content from u-flux
778            DO jj = 2, jpjm1
779               DO ji = fs_2, fs_jpim1   ! vector opt.
780                  IF( ll_clem ) THEN
781                     IF( ll_gurvan ) THEN
782                        zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
783                     ELSE
784                        zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
785                           &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
786                           &         ) * tmask(ji,jj,1)
787                     ENDIF
788                  ELSE
789                     zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
790                  ENDIF
791               END DO
792            END DO
793            CALL lbc_lnk( zzt, 'T', 1. )
794
795         ELSE
796
797            DO jj = 2, jpjm1
798               DO ji = fs_2, fs_jpim1   ! vector opt.
799                  IF( ll_gurvan ) THEN
800                     zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  &
801                        &                   - pt   (ji,jj) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj)
802                  ELSE
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) * pamsk
805                  ENDIF
806                  zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
807               END DO
808            END DO
809            CALL lbc_lnk( zzt, 'T', 1. )
810         ENDIF
811         !
812         !                                                        !--  ultimate interpolation of pt at v-point  --!
813         IF( ll_hoxy ) THEN
814            CALL ultimate_y( kn_umx, pdt, zzt, pv, pvc, pt_v, pfv_ho )
815         ELSE
816            CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
817         ENDIF
818         !                                                        !--  limiter in y --!
819         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
820         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
821         !         
822         !
823      ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
824         !
825         !                                                        !--  ultimate interpolation of pt at v-point  --!
826         CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
827         !                                                        !--  limiter in y --!
828         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho )
829         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups )
830         !                                                        !--  advective form update in zzt  --!
831         IF( ll_1stguess_clem ) THEN
832           
833            ! first guess of tracer content from v-flux
834            DO jj = 2, jpjm1
835               DO ji = fs_2, fs_jpim1   ! vector opt.
836                  IF( ll_clem ) THEN
837                     IF( ll_gurvan ) THEN
838                        zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
839                     ELSE
840                        zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) &
841                           &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
842                           &         ) * tmask(ji,jj,1)
843                     ENDIF
844                  ELSE
845                     zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) &
846                        &         * tmask(ji,jj,1)
847                  ENDIF
848                END DO
849            END DO
850            CALL lbc_lnk( zzt, 'T', 1. )
851           
852         ELSE
853           
854            DO jj = 2, jpjm1
855               DO ji = fs_2, fs_jpim1
856                  IF( ll_gurvan ) THEN
857                     zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  &
858                        &                   - pt   (ji,jj) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj)
859                  ELSE
860                     zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  &
861                        &                   - pt   (ji,jj) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk
862                  ENDIF
863                  zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
864               END DO
865            END DO
866            CALL lbc_lnk( zzt, 'T', 1. )
867         ENDIF
868         !
869         !                                                        !--  ultimate interpolation of pt at u-point  --!
870         IF( ll_hoxy ) THEN
871            CALL ultimate_x( kn_umx, pdt, zzt, pu, puc, pt_u, pfu_ho )
872         ELSE
873            CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
874         ENDIF
875         !                                                        !--  limiter in x --!
876         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho )
877         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups )
878         !
879         !
880      ENDIF
881
882     
883      IF( kn_limiter == 1 ) THEN
884         IF( .NOT. ll_limiter_it2 ) THEN
885            IF( ll_clem ) THEN
886               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
887            ELSE
888               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
889            ENDIF
890         ELSE
891            zzfu_ho(:,:) = pfu_ho(:,:)
892            zzfv_ho(:,:) = pfv_ho(:,:)
893            ! 1st iteration of nonosc (limit the flux with the upstream solution)
894            IF( ll_clem ) THEN
895               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho )
896            ELSE
897               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho )
898            ENDIF
899            ! guess after content field with high order
900            DO jj = 2, jpjm1
901               DO ji = fs_2, fs_jpim1
902                  ztra       = - ( zzfu_ho(ji,jj) - zzfu_ho(ji-1,jj) + zzfv_ho(ji,jj) - zzfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj)
903                  zzt(ji,jj) =   ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
904               END DO
905            END DO
906            CALL lbc_lnk( zzt, 'T', 1. )
907            ! 2nd iteration of nonosc (limit the flux with the limited high order solution)
908            IF( ll_clem ) THEN
909               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho )
910            ELSE
911               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho )
912            ENDIF
913         ENDIF
914      ENDIF
915      !
916   END SUBROUTINE macho
917
918
919   SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho )
920      !!---------------------------------------------------------------------
921      !!                    ***  ROUTINE ultimate_x  ***
922      !!     
923      !! **  Purpose :   compute 
924      !!
925      !! **  Method  :   ... ???
926      !!                 TIM = transient interpolation Modeling
927      !!
928      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
929      !!----------------------------------------------------------------------
930      INTEGER                     , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
931      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
932      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu        ! ice i-velocity component
933      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc       ! ice i-velocity * A component
934      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields
935      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u      ! tracer at u-point
936      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfu_ho    ! high order flux
937      !
938      INTEGER  ::   ji, jj             ! dummy loop indices
939      REAL(wp) ::   zcu, zdx2, zdx4    !   -      -
940      REAL(wp), DIMENSION(jpi,jpj) ::   ztu1, ztu2, ztu3, ztu4
941      !!----------------------------------------------------------------------
942      !
943      !                                                     !--  Laplacian in i-direction  --!
944      DO jj = 2, jpjm1         ! First derivative (gradient)
945         DO ji = 1, fs_jpim1
946            ztu1(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
947         END DO
948         !                     ! Second derivative (Laplacian)
949         DO ji = fs_2, fs_jpim1
950            ztu2(ji,jj) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj)
951         END DO
952      END DO
953      CALL lbc_lnk( ztu2, 'T', 1. )
954      !
955      !                                                     !--  BiLaplacian in i-direction  --!
956      DO jj = 2, jpjm1         ! Third derivative
957         DO ji = 1, fs_jpim1
958            ztu3(ji,jj) = ( ztu2(ji+1,jj) - ztu2(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
959         END DO
960         !                     ! Fourth derivative
961         DO ji = fs_2, fs_jpim1
962            ztu4(ji,jj) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj)
963         END DO
964      END DO
965      CALL lbc_lnk( ztu4, 'T', 1. )
966      !
967      !
968      SELECT CASE (kn_umx )
969      !
970      CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21)
971         !       
972         DO jj = 1, jpjm1
973            DO ji = 1, fs_jpim1   ! vector opt.
974               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                              pt(ji+1,jj) + pt(ji,jj)   &
975                  &                                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) )
976            END DO
977         END DO
978         !
979      CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23)
980         !
981         DO jj = 1, jpjm1
982            DO ji = 1, fs_jpim1   ! vector opt.
983               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
984               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                                   pt(ji+1,jj) + pt(ji,jj)   &
985                  &                                               -              zcu   * ( pt(ji+1,jj) - pt(ji,jj) ) ) 
986            END DO
987         END DO
988         
989      CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24)
990         !
991         DO jj = 1, jpjm1
992            DO ji = 1, fs_jpim1   ! vector opt.
993               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
994               zdx2 = e1u(ji,jj) * e1u(ji,jj)
995!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
996               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (         (                         pt  (ji+1,jj) + pt  (ji,jj)        &
997                  &                                               -              zcu   * ( pt  (ji+1,jj) - pt  (ji,jj) )  )   &
998                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                         ztu2(ji+1,jj) + ztu2(ji,jj)        &
999                  &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj) - ztu2(ji,jj) )  )   )
1000            END DO
1001         END DO
1002         !
1003      CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27)
1004         !
1005         DO jj = 1, jpjm1
1006            DO ji = 1, fs_jpim1   ! vector opt.
1007               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
1008               zdx2 = e1u(ji,jj) * e1u(ji,jj)
1009!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
1010               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (         (                   pt  (ji+1,jj) + pt  (ji,jj)        &
1011                  &                                               -          zcu * ( pt  (ji+1,jj) - pt  (ji,jj) )  )   &
1012                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                   ztu2(ji+1,jj) + ztu2(ji,jj)        &
1013                  &                                               - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) )  )   )
1014            END DO
1015         END DO
1016         !
1017      CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29)
1018         !
1019         DO jj = 1, jpjm1
1020            DO ji = 1, fs_jpim1   ! vector opt.
1021               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
1022               zdx2 = e1u(ji,jj) * e1u(ji,jj)
1023!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
1024               zdx4 = zdx2 * zdx2
1025               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (               (                   pt  (ji+1,jj) + pt  (ji,jj)       &
1026                  &                                                     -          zcu * ( pt  (ji+1,jj) - pt  (ji,jj) ) )   &
1027                  &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *     (                   ztu2(ji+1,jj) + ztu2(ji,jj)       &
1028                  &                                                     - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) )   &
1029                  &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj) + ztu4(ji,jj)       &
1030                  &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj) - ztu4(ji,jj) ) ) )
1031            END DO
1032         END DO
1033         !
1034      END SELECT
1035      !                                                     !-- High order flux in i-direction  --!
1036      IF( ll_neg ) THEN
1037         DO jj = 1, jpjm1
1038            DO ji = 1, fs_jpim1
1039               IF( pt_u(ji,jj) < 0._wp ) THEN
1040                  pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                              pt(ji+1,jj) + pt(ji,jj)   &
1041                     &                                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) )
1042               ENDIF
1043            END DO
1044         END DO
1045      ENDIF
1046
1047      DO jj = 1, jpjm1
1048         DO ji = 1, fs_jpim1   ! vector opt.
1049            IF( ll_clem ) THEN
1050               pfu_ho(ji,jj) = pu(ji,jj) * pt_u(ji,jj)
1051            ELSE
1052               pfu_ho(ji,jj) = puc(ji,jj) * pt_u(ji,jj)
1053            ENDIF
1054         END DO
1055      END DO
1056      !
1057   END SUBROUTINE ultimate_x
1058   
1059 
1060   SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho )
1061      !!---------------------------------------------------------------------
1062      !!                    ***  ROUTINE ultimate_y  ***
1063      !!     
1064      !! **  Purpose :   compute 
1065      !!
1066      !! **  Method  :   ... ???
1067      !!                 TIM = transient interpolation Modeling
1068      !!
1069      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
1070      !!----------------------------------------------------------------------
1071      INTEGER                     , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2)
1072      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
1073      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pv        ! ice j-velocity component
1074      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvc       ! ice j-velocity*A component
1075      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields
1076      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_v      ! tracer at v-point
1077      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pfv_ho    ! high order flux
1078      !
1079      INTEGER  ::   ji, jj       ! dummy loop indices
1080      REAL(wp) ::   zcv, zdy2, zdy4    !   -      -
1081      REAL(wp), DIMENSION(jpi,jpj) ::   ztv1, ztv2, ztv3, ztv4
1082      !!----------------------------------------------------------------------
1083      !
1084      !                                                     !--  Laplacian in j-direction  --!
1085      DO jj = 1, jpjm1         ! First derivative (gradient)
1086         DO ji = fs_2, fs_jpim1
1087            ztv1(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
1088         END DO
1089      END DO
1090      DO jj = 2, jpjm1         ! Second derivative (Laplacian)
1091         DO ji = fs_2, fs_jpim1
1092            ztv2(ji,jj) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj)
1093         END DO
1094      END DO
1095      CALL lbc_lnk( ztv2, 'T', 1. )
1096      !
1097      !                                                     !--  BiLaplacian in j-direction  --!
1098      DO jj = 1, jpjm1         ! First derivative
1099         DO ji = fs_2, fs_jpim1
1100            ztv3(ji,jj) = ( ztv2(ji,jj+1) - ztv2(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
1101         END DO
1102      END DO
1103      DO jj = 2, jpjm1         ! Second derivative
1104         DO ji = fs_2, fs_jpim1
1105            ztv4(ji,jj) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj)
1106         END DO
1107      END DO
1108      CALL lbc_lnk( ztv4, 'T', 1. )
1109      !
1110      !
1111      SELECT CASE (kn_umx )
1112      !
1113      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21)
1114         DO jj = 1, jpjm1
1115            DO ji = 1, fs_jpim1
1116               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                             ( pt(ji,jj+1) + pt(ji,jj) )  &
1117                  &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) )
1118            END DO
1119         END DO
1120         !
1121      CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23)
1122         DO jj = 1, jpjm1
1123            DO ji = 1, fs_jpim1
1124               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1125               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (        ( pt(ji,jj+1) + pt(ji,jj) )  &
1126                  &                                     - zcv * ( pt(ji,jj+1) - pt(ji,jj) ) )
1127            END DO
1128         END DO
1129         CALL lbc_lnk( pt_v, 'V',  1. )
1130         !
1131      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24)
1132         DO jj = 1, jpjm1
1133            DO ji = 1, fs_jpim1
1134               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1135               zdy2 = e2v(ji,jj) * e2v(ji,jj)
1136!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
1137               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                 ( pt  (ji,jj+1) + pt  (ji,jj)       &
1138                  &                                     -                        zcv   * ( pt  (ji,jj+1) - pt  (ji,jj) ) )   &
1139                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1) + ztv2(ji,jj)       &
1140                  &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) )
1141            END DO
1142         END DO
1143         !
1144      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27)
1145         DO jj = 1, jpjm1
1146            DO ji = 1, fs_jpim1
1147               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1148               zdy2 = e2v(ji,jj) * e2v(ji,jj)
1149!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
1150               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                           ( pt  (ji,jj+1) + pt  (ji,jj)       &
1151                  &                                               -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) ) )   &
1152                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                   ztv2(ji,jj+1) + ztv2(ji,jj)       &
1153                  &                                               - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) )
1154            END DO
1155         END DO
1156         !
1157      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29)
1158         DO jj = 1, jpjm1
1159            DO ji = 1, fs_jpim1
1160               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
1161               zdy2 = e2v(ji,jj) * e2v(ji,jj)
1162!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
1163               zdy4 = zdy2 * zdy2
1164               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                 ( pt  (ji,jj+1) + pt  (ji,jj)      &
1165                  &                                                     -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) ) )  &
1166                  &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) *     (                   ztv2(ji,jj+1) + ztv2(ji,jj)      &
1167                  &                                                     - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) )  &
1168                  &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1) + ztv4(ji,jj)      &
1169                  &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1) - ztv4(ji,jj) ) ) )
1170            END DO
1171         END DO
1172         !
1173      END SELECT
1174      !                                                     !-- High order flux in j-direction  --!
1175      IF( ll_neg ) THEN
1176         DO jj = 1, jpjm1
1177            DO ji = 1, fs_jpim1
1178               IF( pt_v(ji,jj) < 0._wp ) THEN
1179                  pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                             ( pt(ji,jj+1) + pt(ji,jj) )  &
1180                     &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) )
1181               ENDIF
1182            END DO
1183         END DO
1184      ENDIF
1185
1186      DO jj = 1, jpjm1
1187         DO ji = 1, fs_jpim1   ! vector opt.
1188            IF( ll_clem ) THEN
1189               pfv_ho(ji,jj) = pv(ji,jj) * pt_v(ji,jj)
1190            ELSE
1191               pfv_ho(ji,jj) = pvc(ji,jj) * pt_v(ji,jj)
1192            ENDIF
1193         END DO
1194      END DO
1195      !
1196   END SUBROUTINE ultimate_y
1197     
1198
1199   SUBROUTINE nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho )
1200      !!---------------------------------------------------------------------
1201      !!                    ***  ROUTINE nonosc  ***
1202      !!     
1203       !! **  Purpose :   compute monotonic tracer fluxes from the upstream
1204      !!       scheme and the before field by a nonoscillatory algorithm
1205      !!
1206      !! **  Method  :   ... ???
1207      !!       warning : pt and pt_low must be masked, but the boundaries
1208      !!       conditions on the fluxes are not necessary zalezak (1979)
1209      !!       drange (1995) multi-dimensional forward-in-time and upstream-
1210      !!       in-space based differencing for fluid
1211      !!----------------------------------------------------------------------
1212      REAL(wp)                     , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0)
1213      REAL(wp)                     , INTENT(in   ) ::   pdt              ! tracer time-step
1214      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pu               ! ice i-velocity => u*e2
1215      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   puc              ! ice i-velocity *A => u*e2*a
1216      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pv               ! ice j-velocity => v*e1
1217      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pvc              ! ice j-velocity *A => v*e1*a
1218      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   ptc, pt, pt_low  ! before field & upstream guess of after field
1219      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   pfv_low, pfu_low ! upstream flux
1220      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux
1221      !
1222      INTEGER  ::   ji, jj    ! dummy loop indices
1223      REAL(wp) ::   zpos, zneg, zbig, zsml, z1_dt, zpos2, zneg2   ! local scalars
1224      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign, zcoef        !   -      -
1225      REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zti_low, ztj_low, zzt
1226      !!----------------------------------------------------------------------
1227      zbig = 1.e+40_wp
1228      zsml = epsi20
1229
1230      IF( ll_zeroup2 ) THEN
1231         DO jj = 1, jpjm1
1232            DO ji = 1, fs_jpim1   ! vector opt.
1233               IF( amaxu(ji,jj) == 0._wp )   pfu_ho(ji,jj) = 0._wp
1234               IF( amaxv(ji,jj) == 0._wp )   pfv_ho(ji,jj) = 0._wp
1235            END DO
1236         END DO
1237      ENDIF
1238     
1239      IF( ll_zeroup4 ) THEN
1240         DO jj = 1, jpjm1
1241            DO ji = 1, fs_jpim1   ! vector opt.
1242               IF( pfu_low(ji,jj) == 0._wp )   pfu_ho(ji,jj) = 0._wp
1243               IF( pfv_low(ji,jj) == 0._wp )   pfv_ho(ji,jj) = 0._wp
1244            END DO
1245         END DO
1246      ENDIF
1247
1248
1249      IF( ll_zeroup1 ) THEN
1250         DO jj = 2, jpjm1
1251            DO ji = fs_2, fs_jpim1
1252               IF( ll_gurvan ) THEN
1253                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
1254                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
1255               ELSE
1256                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
1257                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  &
1258                     &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1259                     &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1260                     &         ) * tmask(ji,jj,1)
1261               ENDIF
1262               IF( zzt(ji,jj) < 0._wp ) THEN
1263                  pfu_ho(ji,jj)   = pfu_low(ji,jj)
1264                  pfv_ho(ji,jj)   = pfv_low(ji,jj)
1265                  WRITE(numout,*) '*** 1 negative high order zzt ***',ji,jj,zzt(ji,jj)
1266               ENDIF
1267!!               IF( ji==26 .AND. jj==86) THEN
1268!!                  WRITE(numout,*) 'zzt high order',zzt(ji,jj)
1269!!                  WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt
1270!!                  WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt
1271!!                  WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt
1272!!                  WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt
1273!!               ENDIF
1274               IF( ll_gurvan ) THEN
1275                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
1276                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
1277               ELSE
1278                  zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
1279                     &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  &
1280                     &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1281                     &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1282                     &         ) * 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) &
1297                     &         ) * tmask(ji,jj,1)
1298               ENDIF
1299               IF( zzt(ji,jj) < 0._wp ) THEN
1300                  WRITE(numout,*) '*** 3 negative high order zzt ***',ji,jj,zzt(ji,jj)
1301               ENDIF
1302            END DO
1303         END DO
1304         CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. )
1305      ENDIF
1306
1307     
1308      ! antidiffusive flux : high order minus low order
1309      ! --------------------------------------------------
1310      DO jj = 1, jpjm1
1311         DO ji = 1, fs_jpim1   ! vector opt.
1312            pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_low(ji,jj)
1313            pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_low(ji,jj)
1314          END DO
1315      END DO
1316
1317      ! extreme case where pfu_ho has to be zero
1318      ! ----------------------------------------
1319      !                                    pfu_ho
1320      !                           *         --->
1321      !                        |      |  *   |        |
1322      !                        |      |      |    *   |   
1323      !                        |      |      |        |    *
1324      !            t_low :       i-1     i       i+1       i+2   
1325      IF( ll_prelimiter_zalesak ) THEN
1326         
1327         DO jj = 2, jpjm1
1328            DO ji = fs_2, fs_jpim1 
1329               zti_low(ji,jj)= pt_low(ji+1,jj  )
1330               ztj_low(ji,jj)= pt_low(ji  ,jj+1)
1331            END DO
1332         END DO
1333         CALL lbc_lnk_multi( zti_low, 'T', 1., ztj_low, 'T', 1. )
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) ) &
1443                  &    ) * ( 1. - pamsk )
1444               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) ) &
1445                  &    ) * ( 1. - pamsk )
1446            ELSE
1447               zneg2 = 0. ; zpos2 = 0.
1448            ENDIF
1449            !
1450            !                                  ! up & down beta terms
1451            IF( (zpos+zpos2) > 0. ) THEN ; zbetup(ji,jj) = MAX( 0._wp, zup - pt_low(ji,jj) ) / (zpos+zpos2) * e1e2t(ji,jj) * z1_dt
1452            ELSE                         ; zbetup(ji,jj) = 0. ! zbig
1453            ENDIF
1454            !
1455            IF( (zneg+zneg2) > 0. ) THEN ; zbetdo(ji,jj) = MAX( 0._wp, pt_low(ji,jj) - zdo ) / (zneg+zneg2) * e1e2t(ji,jj) * z1_dt
1456            ELSE                         ; zbetdo(ji,jj) = 0. ! zbig
1457            ENDIF
1458            !
1459            ! if all the points are outside ice cover
1460            IF( zup == -zbig )   zbetup(ji,jj) = 0. ! zbig
1461            IF( zdo ==  zbig )   zbetdo(ji,jj) = 0. ! zbig           
1462            !
1463
1464!!            IF( ji==26 .AND. jj==86) THEN
1465!               WRITE(numout,*) '-----------------'
1466!               WRITE(numout,*) 'zpos',zpos,zpos2
1467!               WRITE(numout,*) 'zneg',zneg,zneg2
1468!               WRITE(numout,*) 'puc/pu',ABS(puc(ji,jj))/MAX(epsi20, ABS(pu(ji,jj)))
1469!               WRITE(numout,*) 'pvc/pv',ABS(pvc(ji,jj))/MAX(epsi20, ABS(pv(ji,jj)))
1470!               WRITE(numout,*) 'pucm1/pu',ABS(puc(ji-1,jj))/MAX(epsi20, ABS(pu(ji-1,jj)))
1471!               WRITE(numout,*) 'pvcm1/pv',ABS(pvc(ji,jj-1))/MAX(epsi20, ABS(pv(ji,jj-1)))
1472!               WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)+pfu_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt
1473!               WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)+pfv_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt
1474!               WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)+pfu_low(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt
1475!               WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)+pfv_low(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt
1476!               WRITE(numout,*) 'pfu_low',pfu_low(ji,jj) * r1_e1e2t(ji,jj) * pdt
1477!               WRITE(numout,*) 'pfv_low',pfv_low(ji,jj) * r1_e1e2t(ji,jj) * pdt
1478!               WRITE(numout,*) 'pfu_lowm1',pfu_low(ji-1,jj) * r1_e1e2t(ji,jj) * pdt
1479!               WRITE(numout,*) 'pfv_lowm1',pfv_low(ji,jj-1) * r1_e1e2t(ji,jj) * pdt
1480!               
1481!               WRITE(numout,*) 'pt',pt(ji,jj)
1482!               WRITE(numout,*) 'ptim1',pt(ji-1,jj)
1483!               WRITE(numout,*) 'ptjm1',pt(ji,jj-1)
1484!               WRITE(numout,*) 'pt_low',pt_low(ji,jj)
1485!               WRITE(numout,*) 'zbetup',zbetup(ji,jj)
1486!               WRITE(numout,*) 'zbetdo',zbetdo(ji,jj)
1487!               WRITE(numout,*) 'zup',zup
1488!               WRITE(numout,*) 'zdo',zdo
1489!            ENDIF
1490            !
1491         END DO
1492      END DO
1493      CALL lbc_lnk_multi( zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
1494
1495     
1496      ! monotonic flux in the y direction
1497      ! ---------------------------------
1498      DO jj = 1, jpjm1
1499         DO ji = 1, fs_jpim1   ! vector opt.
1500            zau = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji+1,jj) )
1501            zbu = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji+1,jj) )
1502            zcu = 0.5  + SIGN( 0.5 , pfu_ho(ji,jj) )
1503            !
1504            zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu )
1505           
1506            pfu_ho(ji,jj) = pfu_ho(ji,jj) * zcoef + pfu_low(ji,jj)
1507
1508!!            IF( ji==26 .AND. jj==86) THEN
1509!!               WRITE(numout,*) 'coefU',zcoef
1510!!               WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt
1511!!               WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt
1512!!            ENDIF
1513
1514         END DO
1515      END DO
1516
1517      DO jj = 1, jpjm1
1518         DO ji = 1, fs_jpim1   ! vector opt.
1519            zav = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji,jj+1) )
1520            zbv = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji,jj+1) )
1521            zcv = 0.5  + SIGN( 0.5 , pfv_ho(ji,jj) )
1522            !
1523            zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv )
1524           
1525            pfv_ho(ji,jj) = pfv_ho(ji,jj) * zcoef + pfv_low(ji,jj)
1526
1527!!            IF( ji==26 .AND. jj==86) THEN
1528!!               WRITE(numout,*) 'coefV',zcoef
1529!!               WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt
1530!!               WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt
1531!!            ENDIF
1532         END DO
1533      END DO
1534
1535      ! clem test
1536      DO jj = 2, jpjm1
1537         DO ji = 2, fs_jpim1   ! vector opt.
1538            IF( ll_gurvan ) THEN
1539               zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
1540                  &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
1541            ELSE
1542               zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj)  &
1543                  &                     - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj)  &
1544                  &                     + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1545                  &                     + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &
1546                  &         ) * 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.