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.
Changeset 10419 for NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_adv_umx.F90 – NEMO

Ignore:
Timestamp:
2018-12-19T20:46:30+01:00 (5 years ago)
Author:
smasson
Message:

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10418, see #2133

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_adv_umx.F90

    r10397 r10419  
    1414   !!   ultimate_x(_y)    : compute a tracer value at velocity points using ULTIMATE scheme at various orders 
    1515   !!   macho             : ??? 
    16    !!   nonosc_2d         : compute monotonic tracer fluxes by a non-oscillatory algorithm  
     16   !!   nonosc            : compute monotonic tracer fluxes by a non-oscillatory algorithm  
    1717   !!---------------------------------------------------------------------- 
    1818   USE phycst         ! physical constant 
    … …  
    2020   USE sbc_oce , ONLY : nn_fsbc   ! update frequency of surface boundary condition 
    2121   USE ice            ! sea-ice variables 
     22   USE icevar         ! sea-ice: operations 
    2223   ! 
    2324   USE in_out_manager ! I/O manager 
    … …  
    3435   REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120 
    3536 
     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 
    3683   !! * Substitutions 
    3784#  include "vectopt_loop_substitute.h90" 
    … …  
    3986   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    4087   !! $Id$ 
    41    !! Software governed by the CeCILL license (see ./LICENSE) 
     88   !! Software governed by the CeCILL licence     (./LICENSE) 
    4289   !!---------------------------------------------------------------------- 
    4390CONTAINS 
    4491 
    45    SUBROUTINE ice_dyn_adv_umx( k_order, kt, pu_ice, pv_ice,  & 
    46       &                    pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     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 ) 
    4794      !!---------------------------------------------------------------------- 
    4895      !!                  ***  ROUTINE ice_dyn_adv_umx  *** 
    … …  
    54101      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    55102      !!---------------------------------------------------------------------- 
    56       INTEGER                     , INTENT(in   ) ::   k_order    ! order of the scheme (1-5 or 20) 
     103      INTEGER                     , INTENT(in   ) ::   kn_umx     ! order of the scheme (1-5=UM or 20=CEN2) 
    57104      INTEGER                     , INTENT(in   ) ::   kt         ! time step 
    58105      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity 
    … …  
    70117      ! 
    71118      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices 
    72       INTEGER  ::   initad                  ! number of sub-timestep for the advection 
    73       INTEGER  ::   ipl                     ! third dimention of tracer array 
    74  
    75       REAL(wp) ::   zusnit, zdt 
    76       REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! send zcflnow and receive zcflprv 
    77       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zudy, zvdx, zcu_box, zcv_box 
    78       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zpato 
     119      INTEGER  ::   icycle                  ! number of sub-timestep for the advection 
     120      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers 
     121      REAL(wp) ::   zdt 
     122      REAL(wp), DIMENSION(1)       ::   zcflprv, zcflnow   ! send zcflnow and receive zcflprv 
     123      REAL(wp), DIMENSION(jpi,jpj) ::   zudy, zvdx 
     124      REAL(wp), DIMENSION(jpi,jpj) ::   zati1, zati2 
     125 
     126 
     127 
     128      REAL(wp), DIMENSION(jpi,jpj)     ::   zcu_box, zcv_box 
     129      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho 
     130      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai, z1_aip 
     131      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhvar 
     132 
    79133      !!---------------------------------------------------------------------- 
    80134      ! 
    81135      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' 
    82136      ! 
    83       ALLOCATE( zudy(jpi,jpj) , zvdx(jpi,jpj) , zcu_box(jpi,jpj) , zcv_box(jpi,jpj) ) 
    84       ALLOCATE( zpato(jpi,jpj,1) ) 
    85137      ! 
    86138      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- ! 
    … …  
    93145      CALL mpp_delay_max( 'icedyn_adv_umx', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 ) 
    94146 
    95       IF( zcflprv(1) > .5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp 
    96       ELSE                         ;   initad = 1   ;   zusnit = 1.0_wp 
     147      IF( zcflprv(1) > .5 ) THEN   ;   icycle = 2 
     148      ELSE                         ;   icycle = 1 
    97149      ENDIF 
    98  
    99       zdt = rdt_ice / REAL(initad) 
     150       
     151      zdt = rdt_ice / REAL(icycle) 
    100152 
    101153      ! --- transport --- ! 
    … …  
    118170      END DO 
    119171 
    120       zpato (:,:,1) = pato_i(:,:) 
    121  
     172      IF( ll_zeroup2 ) THEN 
     173         IF(.NOT. ALLOCATED(amaxu))       ALLOCATE(amaxu (jpi,jpj,jpl)) 
     174         IF(.NOT. ALLOCATED(amaxv))       ALLOCATE(amaxv (jpi,jpj,jpl)) 
     175      ENDIF 
    122176      !---------------! 
    123177      !== advection ==! 
    124178      !---------------! 
    125       DO jt = 1, initad 
     179      DO jt = 1, icycle 
     180 
     181!!$         IF( ll_ADVopw ) THEN 
     182!!$            zamsk = 1._wp 
     183!!$            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) )   ! Open water area  
     184!!$            zamsk = 0._wp 
     185!!$         ELSE 
     186            zati1(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
     187!!$         ENDIF 
    126188          
    127          l_full_nf_update = .FALSE.   ! false: disable full North fold update (performances) 
    128          CALL adv_umx( k_order, kt,   1, zdt, zudy, zvdx, zcu_box, zcv_box, zpato(:,:,1) )        ! Open water area  
    129          CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,:) )         ! Ice area 
    130          CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,:) )         ! Ice  volume 
    131          CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, psv_i(:,:,:) )        ! Salt content 
    132          CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i(:,:,:) )        ! Age content 
     189         WHERE( pa_i(:,:,:) >= epsi20 )   ;   z1_ai(:,:,:) = 1._wp / pa_i(:,:,:) 
     190         ELSEWHERE                        ;   z1_ai(:,:,:) = 0. 
     191         END WHERE 
     192            ! 
     193         WHERE( pa_ip(:,:,:) >= epsi20 )  ;   z1_aip(:,:,:) = 1._wp / pa_ip(:,:,:) 
     194         ELSEWHERE                        ;   z1_aip(:,:,:) = 0. 
     195         END WHERE 
     196         ! 
     197         IF( ll_zeroup2 ) THEN 
     198            DO jl = 1, jpl 
     199               DO jj = 2, jpjm1 
     200                  DO ji = fs_2, fs_jpim1 
     201                     amaxu(ji,jj,jl)=MAX( pa_i(ji,jj,jl), pa_i(ji,jj-1,jl), pa_i(ji,jj+1,jl), & 
     202                        &                                 pa_i(ji+1,jj,jl), pa_i(ji+1,jj-1,jl), pa_i(ji+1,jj+1,jl) ) 
     203                     amaxv(ji,jj,jl)=MAX( pa_i(ji,jj,jl), pa_i(ji-1,jj,jl), pa_i(ji+1,jj,jl), & 
     204                        &                                 pa_i(ji,jj+1,jl), pa_i(ji-1,jj+1,jl), pa_i(ji+1,jj+1,jl) ) 
     205                  END DO 
     206               END DO 
     207            END DO 
     208            CALL lbc_lnk_multi('icedyn_adv_umx', amaxu, 'T', 1., amaxv, 'T', 1.) 
     209         ENDIF 
     210         ! 
     211         DO jl = 1, jpl 
     212            zua_ho(:,:,jl) = zudy(:,:) 
     213            zva_ho(:,:,jl) = zvdx(:,:) 
     214         END DO 
     215          
     216         zamsk = 1._wp 
     217         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_i, pa_i, zua_ho, zva_ho ) ! Ice area 
     218         zamsk = 0._wp 
     219         ! 
     220         IF( ll_thickness ) THEN 
     221            DO jl = 1, jpl 
     222               zua_ho(:,:,jl) = zudy(:,:) 
     223               zva_ho(:,:,jl) = zvdx(:,:) 
     224            END DO 
     225         ENDIF 
     226            ! 
     227         zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 
     228         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_i )    ! Ice volume 
     229         IF( ll_thickness )   pv_i(:,:,:) = zhvar(:,:,:) * pa_i(:,:,:) 
     230         ! 
     231         zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 
     232         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_s )    ! Snw volume 
     233         IF( ll_thickness )   pv_s(:,:,:) = zhvar(:,:,:) * pa_i(:,:,:) 
     234         ! 
     235         zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:) 
     236         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, psv_i )    ! Salt content 
     237         ! 
     238         zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:) 
     239         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, poa_i )    ! Age content 
     240         ! 
    133241         DO jk = 1, nlay_i 
    134             CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,:) )   ! Ice  heat content 
    135          END DO 
    136          CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,:) )         ! Snow volume 
     242            zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:) 
     243            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_i(:,:,jk,:) ) ! Ice heat content 
     244         END DO 
     245         ! 
    137246         DO jk = 1, nlay_s 
    138             CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,jk,:) )   ! Snow heat content 
    139          END DO 
     247            zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:) 
     248            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_s(:,:,jk,:) ) ! Snw heat content 
     249         END DO 
     250            ! 
    140251         IF ( ln_pnd_H12 ) THEN 
    141             CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,:) )     ! Melt pond fraction 
    142             CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,:) )     ! Melt pond volume 
     252            ! 
     253            DO jl = 1, jpl 
     254               zua_ho(:,:,jl) = zudy(:,:) 
     255               zva_ho(:,:,jl) = zvdx(:,:) 
     256            END DO 
     257             
     258            zamsk = 1._wp 
     259            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_ip, pa_ip, zua_ho, zva_ho ) ! mp fraction 
     260            zamsk = 0._wp 
     261            ! 
     262            zhvar(:,:,:) = pv_ip(:,:,:) * z1_ai(:,:,:) 
     263            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_ip ) ! mp volume 
    143264         ENDIF 
    144  
    145          l_full_nf_update = jt == initad   ! true: enable full North fold update (performances) when exiting the loop 
    146          CALL lbc_lnk( 'icedyn_adv_umx', zpato, 'T',  1. ) 
    147          CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T',  1. ) 
    148          CALL lbc_lnk( 'icedyn_adv_umx', pe_s, 'T',  1. ) 
    149          IF ( ln_pnd_H12 ) THEN 
    150             CALL lbc_lnk_multi( 'icedyn_adv_umx',  pa_i, 'T',  1., pv_i, 'T',  1., psv_i, 'T',  1., & 
    151                                                 & poa_i, 'T',  1., pv_s, 'T',  1., pa_ip, 'T',  1., & 
    152                                                 & pv_ip, 'T',  1. ) 
    153          ELSE 
    154             CALL lbc_lnk_multi( 'icedyn_adv_umx',  pa_i, 'T',  1., pv_i, 'T',  1., psv_i, 'T',  1., & 
    155                                                 & poa_i, 'T',  1., pv_s, 'T',  1. ) 
    156          ENDIF 
    157           
     265         ! 
     266         ! 
     267!!$         IF( .NOT. ll_ADVopw ) THEN 
     268            zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
     269            DO jj = 2, jpjm1 
     270               DO ji = fs_2, fs_jpim1 
     271                  pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                                                  ! Open water area 
     272                     &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt 
     273               END DO 
     274            END DO 
     275            CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T',  1. ) 
     276!!$         ENDIF 
     277         ! 
    158278      END DO 
    159279      ! 
    160       pato_i(:,:) = zpato (:,:,1) 
    161       ! 
    162       DEALLOCATE( zudy, zvdx, zcu_box, zcv_box, zpato ) 
    163       ! 
    164280   END SUBROUTINE ice_dyn_adv_umx 
    165281 
    166282    
    167    SUBROUTINE adv_umx( k_order, kt, ipl, pdt, puc, pvc, pubox, pvbox, ptc ) 
     283   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho ) 
    168284      !!---------------------------------------------------------------------- 
    169285      !!                  ***  ROUTINE adv_umx  *** 
    … …  
    178294      !! ** Action : - pt  the after advective tracer 
    179295      !!---------------------------------------------------------------------- 
    180       INTEGER                         , INTENT(in   ) ::   k_order        ! order of the ULTIMATE scheme 
    181       INTEGER                         , INTENT(in   ) ::   kt             ! number of iteration 
    182       INTEGER                         , INTENT(in   ) ::   ipl            ! third dimension of tracer array 
    183       REAL(wp)                        , INTENT(in   ) ::   pdt            ! tracer time-step 
    184       REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   puc  , pvc     ! 2 ice velocity components => u*e2 
    185       REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   pubox, pvbox   ! upstream velocity 
    186       REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(inout) ::   ptc            ! tracer content field 
     296      REAL(wp)                    , INTENT(in   )           ::   pamsk          ! advection of concentration (1) or other tracers (0) 
     297      INTEGER                     , INTENT(in   )           ::   kn_umx         ! order of the scheme (1-5=UM or 20=CEN2) 
     298      INTEGER                     , INTENT(in   )           ::   jt             ! number of sub-iteration 
     299      INTEGER                     , INTENT(in   )           ::   kt             ! number of iteration 
     300      REAL(wp)                    , INTENT(in   )           ::   pdt            ! tracer time-step 
     301      REAL(wp), DIMENSION(:,:  ), INTENT(in   )           ::   pu   , pv      ! 2 ice velocity components => u*e2 
     302      REAL(wp), DIMENSION(:,:,:), INTENT(in   )           ::   puc  , pvc     ! 2 ice velocity components => u*e2 or u*a*e2u 
     303      REAL(wp), DIMENSION(:,:  ), INTENT(in   )           ::   pubox, pvbox   ! upstream velocity 
     304      REAL(wp), DIMENSION(:,:,:), INTENT(inout)           ::   pt             ! tracer field 
     305      REAL(wp), DIMENSION(:,:,:), INTENT(inout)           ::   ptc            ! tracer content field 
     306      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho ! high order u*a fluxes 
    187307      ! 
    188308      INTEGER  ::   ji, jj, jl       ! dummy loop indices   
    189  
    190309      REAL(wp) ::   ztra             ! local scalar 
    191       REAL(wp), DIMENSION(jpi,jpj,ipl) ::   zfu_ups, zfu_ho, zt_u, zt_ups 
    192       REAL(wp), DIMENSION(jpi,jpj,ipl) ::   zfv_ups, zfv_ho, zt_v, ztrd 
    193  
    194       DO jl = 1, ipl 
    195       !!---------------------------------------------------------------------- 
    196       ! 
    197       !  upstream advection with initial mass fluxes & intermediate update 
    198       ! -------------------------------------------------------------------- 
    199       DO jj = 1, jpjm1         ! upstream tracer flux in the i and j direction 
    200          DO ji = 1, fs_jpim1   ! vector opt. 
    201             zfu_ups(ji,jj,jl) = MAX( puc(ji,jj), 0._wp ) * ptc(ji,jj,jl) + MIN( puc(ji,jj), 0._wp ) * ptc(ji+1,jj,jl) 
    202             zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj), 0._wp ) * ptc(ji,jj,jl) + MIN( pvc(ji,jj), 0._wp ) * ptc(ji,jj+1,jl) 
     310      INTEGER  ::   kn_limiter = 1   ! 1=nonosc ; 2=superbee ; 3=h3(rachid) 
     311      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zfu_ho , zfv_ho , zt_u, zt_v, zpt 
     312      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zfu_ups, zfv_ups, zt_ups   ! only for nonosc  
     313      !!---------------------------------------------------------------------- 
     314      ! 
     315      !  upstream (_ups) advection with initial mass fluxes 
     316      ! --------------------------------------------------- 
     317 
     318      IF( ll_gurvan .AND. pamsk==0. ) THEN 
     319         DO jl = 1, jpl 
     320            DO jj = 2, jpjm1 
     321               DO ji = fs_2, fs_jpim1 
     322                  pt(ji,jj,jl) = ( pt (ji,jj,jl) + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj)     & 
     323                     &                           + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) )   & 
     324                     &           * tmask(ji,jj,1) 
     325               END DO 
     326            END DO 
     327         END DO 
     328         CALL lbc_lnk( 'icedyn_adv_umx', pt, 'T', 1. ) 
     329      ENDIF 
     330 
     331       
     332      IF( .NOT. ll_upsxy ) THEN 
     333 
     334         ! fluxes in both x-y directions 
     335         DO jl = 1, jpl 
     336            DO jj = 1, jpjm1 
     337               DO ji = 1, fs_jpim1 
     338                  IF( ll_clem ) THEN 
     339                     zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 
     340                     zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 
     341                  ELSE 
     342                     zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * pt(ji+1,jj,jl) 
     343                     zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj+1,jl) 
     344                  ENDIF 
     345               END DO 
     346            END DO 
     347         END DO 
     348 
     349      ELSE 
     350         ! 
     351         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     352            ! flux in x-direction 
     353            DO jl = 1, jpl 
     354               DO jj = 1, jpjm1 
     355                  DO ji = 1, fs_jpim1 
     356                     IF( ll_clem ) THEN 
     357                        zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 
     358                     ELSE 
     359                        zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * pt(ji+1,jj,jl) 
     360                     ENDIF 
     361                  END DO 
     362               END DO 
     363            END DO 
     364             
     365            ! first guess of tracer content from u-flux 
     366            DO jl = 1, jpl 
     367               DO jj = 2, jpjm1 
     368                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     369                     IF( ll_clem ) THEN 
     370                        IF( ll_gurvan ) THEN 
     371                           zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     372                        ELSE 
     373                           zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     374                              &            + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     375                              &            ) * tmask(ji,jj,1) 
     376                        ENDIF 
     377                     ELSE 
     378                        zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) & 
     379                           &         * tmask(ji,jj,1) 
     380                     ENDIF 
     381                     !!                  IF( ji==26 .AND. jj==86) THEN 
     382                     !!                     WRITE(numout,*) '************************' 
     383                     !!                     WRITE(numout,*) 'zpt upstream',zpt(ji,jj) 
     384                     !!                  ENDIF 
     385                  END DO 
     386               END DO 
     387            END DO 
     388            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
     389            ! 
     390            ! flux in y-direction 
     391            DO jl = 1, jpl 
     392               DO jj = 1, jpjm1 
     393                  DO ji = 1, fs_jpim1 
     394                     IF( ll_clem ) THEN 
     395                        zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl) 
     396                     ELSE 
     397                        zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * zpt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * zpt(ji,jj+1,jl) 
     398                     ENDIF 
     399                  END DO 
     400               END DO 
     401            END DO 
     402         ! 
     403         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
     404            ! flux in y-direction 
     405            DO jl = 1, jpl 
     406               DO jj = 1, jpjm1 
     407                  DO ji = 1, fs_jpim1 
     408                     IF( ll_clem ) THEN 
     409                        zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 
     410                     ELSE 
     411                        zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj+1,jl) 
     412                     ENDIF 
     413                  END DO 
     414               END DO 
     415            END DO 
     416 
     417            ! first guess of tracer content from v-flux 
     418            DO jl = 1, jpl 
     419               DO jj = 2, jpjm1 
     420                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     421                     IF( ll_clem ) THEN 
     422                        IF( ll_gurvan ) THEN 
     423                           zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     424                        ELSE 
     425                           zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 
     426                              &            + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 
     427                              &            * tmask(ji,jj,1) 
     428                        ENDIF 
     429                     ELSE 
     430                        zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) & 
     431                           &            * tmask(ji,jj,1) 
     432                     ENDIF 
     433                     !!                  IF( ji==26 .AND. jj==86) THEN 
     434                     !!                     WRITE(numout,*) '************************' 
     435                     !!                     WRITE(numout,*) 'zpt upstream',zpt(ji,jj) 
     436                     !!                  ENDIF 
     437                  END DO 
     438               END DO 
     439            END DO 
     440            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
     441            ! 
     442            ! flux in x-direction 
     443            DO jl = 1, jpl 
     444               DO jj = 1, jpjm1 
     445                  DO ji = 1, fs_jpim1 
     446                     IF( ll_clem ) THEN 
     447                        zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl) 
     448                     ELSE 
     449                        zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * zpt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * zpt(ji+1,jj,jl) 
     450                     ENDIF 
     451                  END DO 
     452               END DO 
     453            END DO 
     454            ! 
     455         ENDIF 
     456          
     457      ENDIF 
     458 
     459      IF( ll_clem .AND. kn_limiter /= 1 )  & 
     460         & CALL ctl_stop( 'STOP', 'icedyn_adv_umx: ll_clem incompatible with limiters other than nonosc' ) 
     461 
     462      IF( ll_zeroup2 ) THEN 
     463         DO jl = 1, jpl 
     464            DO jj = 1, jpjm1 
     465               DO ji = 1, fs_jpim1   ! vector opt. 
     466                  IF( amaxu(ji,jj,jl) == 0._wp )   zfu_ups(ji,jj,jl) = 0._wp 
     467                  IF( amaxv(ji,jj,jl) == 0._wp )   zfv_ups(ji,jj,jl) = 0._wp 
     468               END DO 
     469            END DO 
     470         END DO 
     471      ENDIF 
     472 
     473      ! guess after content field with upstream scheme 
     474      DO jl = 1, jpl 
     475         DO jj = 2, jpjm1 
     476            DO ji = fs_2, fs_jpim1 
     477               ztra          = - (   zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj  ,jl)   & 
     478                  &                + zfv_ups(ji,jj,jl) - zfv_ups(ji  ,jj-1,jl) ) * r1_e1e2t(ji,jj) 
     479               IF( ll_clem ) THEN 
     480                  IF( ll_gurvan ) THEN 
     481                     zt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1) 
     482                  ELSE 
     483                     zt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + pdt * ztra + ( pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) )   & 
     484                        &                                            +   pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) ) & 
     485                        &                                              * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 
     486                  ENDIF 
     487               ELSE 
     488                  zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1) 
     489               ENDIF 
     490               !!            IF( ji==26 .AND. jj==86) THEN 
     491               !!               WRITE(numout,*) '**************************' 
     492               !!               WRITE(numout,*) 'zt upstream',zt_ups(ji,jj) 
     493               !!            ENDIF 
     494            END DO 
    203495         END DO 
    204496      END DO 
    205        
    206          DO jj = 2, jpjm1            ! total intermediate advective trends 
    207             DO ji = fs_2, fs_jpim1   ! vector opt. 
    208                ztra = - (   zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj  ,jl)   & 
    209                   &       + zfv_ups(ji,jj,jl) - zfv_ups(ji  ,jj-1,jl)   ) * r1_e1e2t(ji,jj) 
    210             ! 
    211                ztrd(ji,jj,jl) =                         ztra                         ! upstream trend [ -div(uh) or -div(uhT) ]   
    212                zt_ups (ji,jj,jl) = ( ptc(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1)   ! guess after content field with monotonic scheme 
    213             END DO 
    214          END DO 
    215       END DO 
    216        
     497      CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. ) 
     498 
    217499      ! High order (_ho) fluxes  
    218500      ! ----------------------- 
    219       SELECT CASE( k_order ) 
    220       CASE ( 20 )                          ! centered second order 
    221          DO jl = 1, ipl 
     501      SELECT CASE( kn_umx ) 
     502         ! 
     503      CASE ( 20 )                          !== centered second order ==! 
     504         ! 
     505         CALL cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho,  & 
     506            &       zt_ups, zfu_ups, zfv_ups ) 
     507         ! 
     508      CASE ( 1:5 )                         !== 1st to 5th order ULTIMATE-MACHO scheme ==! 
     509         ! 
     510         CALL macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, zt_u, zt_v, zfu_ho, zfv_ho,  & 
     511            &        zt_ups, zfu_ups, zfv_ups ) 
     512         ! 
     513      END SELECT 
     514 
     515      IF( ll_thickness ) THEN 
     516         ! final trend with corrected fluxes 
     517         ! ------------------------------------ 
     518         DO jl = 1, jpl 
     519            DO jj = 2, jpjm1 
     520               DO ji = fs_2, fs_jpim1 
     521                  IF( ll_gurvan ) THEN 
     522                     ztra       = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj)  
     523                  ELSE 
     524                     ztra       = ( - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) )  &  
     525                        &           + ( pt(ji,jj,jl) * ( pu(ji,jj) - pu(ji-1,jj) ) * (1.-pamsk) ) & 
     526                        &           + ( pt(ji,jj,jl) * ( pv(ji,jj) - pv(ji,jj-1) ) * (1.-pamsk) ) ) * r1_e1e2t(ji,jj) 
     527                  ENDIF 
     528                  pt(ji,jj,jl) = ( pt(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1) 
     529                   
     530                  IF( pt(ji,jj,jl) < -epsi20 ) THEN 
     531                     WRITE(numout,*) 'T<0 ',pt(ji,jj,jl) 
     532                  ENDIF 
     533                   
     534                  IF( pt(ji,jj,jl) < 0._wp .AND. pt(ji,jj,jl) >= -epsi20 )   pt(ji,jj,jl) = 0._wp 
     535                   
     536                  !!               IF( ji==26 .AND. jj==86) THEN 
     537                  !!                  WRITE(numout,*) 'zt high order',pt(ji,jj) 
     538                  !!               ENDIF 
     539               END DO 
     540            END DO 
     541         END DO 
     542         CALL lbc_lnk( 'icedyn_adv_umx', pt, 'T',  1. ) 
     543      ENDIF 
     544    
     545      ! Rachid trick 
     546      ! ------------ 
     547      IF( ll_clem ) THEN 
     548         IF( pamsk == 0. ) THEN 
     549            DO jl = 1, jpl 
     550               DO jj = 1, jpjm1 
     551                  DO ji = 1, fs_jpim1 
     552                     IF( ABS( puc(ji,jj,jl) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN 
     553                        zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj) 
     554                        zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj) 
     555                     ELSE 
     556                        zfu_ho (ji,jj,jl) = 0._wp 
     557                        zfu_ups(ji,jj,jl) = 0._wp 
     558                     ENDIF 
     559                     ! 
     560                     IF( ABS( pvc(ji,jj,jl) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN 
     561                        zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj) 
     562                        zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj) 
     563                     ELSE 
     564                        zfv_ho (ji,jj,jl) = 0._wp   
     565                        zfv_ups(ji,jj,jl) = 0._wp   
     566                     ENDIF 
     567                  END DO 
     568               END DO 
     569            END DO 
     570         ENDIF 
     571      ENDIF 
     572 
     573      IF( ll_zeroup5 ) THEN 
     574         DO jl = 1, jpl 
     575            DO jj = 2, jpjm1 
     576               DO ji = 2, fs_jpim1   ! vector opt. 
     577                  zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     578                     &                            - ( zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  ) * tmask(ji,jj,1) 
     579                  IF( zpt(ji,jj,jl) < 0. ) THEN 
     580                     zfu_ho(ji  ,jj,jl) = zfu_ups(ji  ,jj,jl) 
     581                     zfu_ho(ji-1,jj,jl) = zfu_ups(ji-1,jj,jl) 
     582                     zfv_ho(ji  ,jj,jl) = zfv_ups(ji  ,jj,jl) 
     583                     zfv_ho(ji,jj-1,jl) = zfv_ups(ji,jj-1,jl) 
     584                  ENDIF 
     585               END DO 
     586            END DO 
     587         END DO 
     588         CALL lbc_lnk_multi( 'icedyn_adv_umx', zfu_ho, 'U',  -1., zfv_ho, 'V',  -1. ) 
     589      ENDIF 
     590 
     591      ! output high order fluxes u*a 
     592      ! ---------------------------- 
     593      IF( PRESENT( pua_ho ) ) THEN 
     594         DO jl = 1, jpl 
    222595            DO jj = 1, jpjm1 
    223                DO ji = 1, fs_jpim1   ! vector opt. 
    224                   zfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj) * ( ptc(ji,jj,jl) + ptc(ji+1,jj,jl) ) 
    225                   zfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj) * ( ptc(ji,jj,jl) + ptc(ji,jj+1,jl) ) 
    226                END DO 
    227             END DO 
    228          END DO 
    229          ! 
    230       CASE ( 1:5 )                      ! 1st to 5th order ULTIMATE-MACHO scheme 
    231          CALL macho( k_order, kt, ipl, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v ) 
    232          ! 
    233          DO jl = 1, ipl 
     596               DO ji = 1, fs_jpim1 
     597                  pua_ho(ji,jj,jl) = zfu_ho(ji,jj,jl) 
     598                  pva_ho(ji,jj,jl) = zfv_ho(ji,jj,jl) 
     599               END DO 
     600            END DO 
     601         END DO 
     602      ENDIF 
     603 
     604 
     605      IF( .NOT.ll_thickness ) THEN 
     606         ! final trend with corrected fluxes 
     607         ! ------------------------------------ 
     608         DO jl = 1, jpl 
    234609            DO jj = 2, jpjm1 
    235                DO ji = 1, fs_jpim1   ! vector opt. 
    236                   zfu_ho(ji,jj,jl) = puc(ji,jj) * zt_u(ji,jj,jl) 
    237                END DO 
    238             END DO 
    239          END DO 
    240          DO jl = 1, ipl 
     610               DO ji = fs_2, fs_jpim1  
     611                  ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj) * pdt   
     612                   
     613                  ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra ) * tmask(ji,jj,1) 
     614                   
     615                  !!               IF( ji==26 .AND. jj==86) THEN 
     616                  !!                  WRITE(numout,*) 'ztc high order',ptc(ji,jj) 
     617                  !!               ENDIF 
     618                   
     619               END DO 
     620            END DO 
     621         END DO 
     622         CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T',  1. ) 
     623      ENDIF 
     624       
     625      ! 
     626   END SUBROUTINE adv_umx 
     627 
     628   SUBROUTINE cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, & 
     629      &             pt_ups, pfu_ups, pfv_ups ) 
     630      !!--------------------------------------------------------------------- 
     631      !!                    ***  ROUTINE macho  *** 
     632      !!      
     633      !! **  Purpose :   compute   
     634      !! 
     635      !! **  Method  :   ... ??? 
     636      !!                 TIM = transient interpolation Modeling  
     637      !! 
     638      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
     639      !!---------------------------------------------------------------------- 
     640      REAL(wp)                    , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0) 
     641      INTEGER                     , INTENT(in   ) ::   kn_limiter       ! limiter 
     642      INTEGER                     , INTENT(in   ) ::   jt               ! number of sub-iteration 
     643      INTEGER                     , INTENT(in   ) ::   kt               ! number of iteration 
     644      REAL(wp)                    , INTENT(in   ) ::   pdt              ! tracer time-step 
     645      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt               ! tracer fields 
     646      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components 
     647      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components 
     648      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptc              ! tracer content at before time step  
     649      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes  
     650      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content  
     651      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes  
     652      ! 
     653      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
     654      LOGICAL  ::   ll_xy = .TRUE. 
     655      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zzt 
     656      !!---------------------------------------------------------------------- 
     657      ! 
     658      IF( .NOT.ll_xy ) THEN   !-- no alternate directions --! 
     659         ! 
     660         DO jl = 1, jpl 
    241661            DO jj = 1, jpjm1 
    242                DO ji = fs_2, fs_jpim1   ! vector opt. 
    243                   zfv_ho(ji,jj,jl) = pvc(ji,jj) * zt_v(ji,jj,jl) 
    244                END DO 
    245             END DO 
    246          END DO 
    247          ! 
    248       END SELECT 
     662               DO ji = 1, fs_jpim1 
     663                  IF( ll_clem ) THEN 
     664                     pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 
     665                     pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 
     666                  ELSE 
     667                     pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 
     668                     pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 
     669                  ENDIF 
     670               END DO 
     671            END DO 
     672         END DO 
     673         IF    ( kn_limiter == 1 ) THEN 
     674            IF( ll_clem ) THEN 
     675               CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     676            ELSE 
     677               CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     678            ENDIF 
     679         ELSEIF( kn_limiter == 2 ) THEN 
     680            CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 
     681            CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 
     682         ELSEIF( kn_limiter == 3 ) THEN 
     683            CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
     684            CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     685         ENDIF 
     686         ! 
     687      ELSE                    !-- alternate directions --! 
     688         ! 
     689         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     690            ! 
     691            ! flux in x-direction 
     692            DO jl = 1, jpl 
     693               DO jj = 1, jpjm1 
     694                  DO ji = 1, fs_jpim1 
     695                     IF( ll_clem ) THEN 
     696                        pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 
     697                     ELSE 
     698                        pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 
     699                     ENDIF 
     700                  END DO 
     701               END DO 
     702            END DO 
     703            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 
     704            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
     705 
     706            ! first guess of tracer content from u-flux 
     707            DO jl = 1, jpl 
     708               DO jj = 2, jpjm1 
     709                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     710                     IF( ll_clem ) THEN 
     711                        zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)        & 
     712                           &            + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 
     713                           &         * tmask(ji,jj,1) 
     714                     ELSE                      
     715                        zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     716                     ENDIF 
     717                  END DO 
     718               END DO 
     719            END DO 
     720            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
     721 
     722            ! flux in y-direction 
     723            DO jl = 1, jpl 
     724               DO jj = 1, jpjm1 
     725                  DO ji = 1, fs_jpim1 
     726                     IF( ll_clem ) THEN 
     727                        pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji,jj+1,jl) ) 
     728                     ELSE                      
     729                        pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( zzt(ji,jj,jl) + zzt(ji,jj+1,jl) ) 
     730                     ENDIF 
     731                  END DO 
     732               END DO 
     733            END DO 
     734            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 
     735            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     736 
     737         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
     738            ! 
     739            ! flux in y-direction 
     740            DO jl = 1, jpl 
     741               DO jj = 1, jpjm1 
     742                  DO ji = 1, fs_jpim1 
     743                     IF( ll_clem ) THEN 
     744                        pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 
     745                     ELSE                      
     746                        pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 
     747                     ENDIF 
     748                  END DO 
     749               END DO 
     750            END DO 
     751            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 
     752            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     753            ! 
     754            ! first guess of tracer content from v-flux 
     755            DO jl = 1, jpl 
     756               DO jj = 2, jpjm1 
     757                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     758                     IF( ll_clem ) THEN 
     759                        zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 
     760                           &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 
     761                           &         * tmask(ji,jj,1) 
     762                     ELSE 
     763                        zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     764                     ENDIF 
     765                  END DO 
     766               END DO 
     767            END DO 
     768            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
     769            ! 
     770            ! flux in x-direction 
     771            DO jl = 1, jpl 
     772               DO jj = 1, jpjm1 
     773                  DO ji = 1, fs_jpim1 
     774                     IF( ll_clem ) THEN 
     775                        pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji+1,jj,jl) ) 
     776                     ELSE 
     777                        pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( zzt(ji,jj,jl) + zzt(ji+1,jj,jl) ) 
     778                     ENDIF 
     779                  END DO 
     780               END DO 
     781            END DO 
     782            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 
     783            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
     784 
     785         ENDIF 
     786         IF( ll_clem ) THEN 
     787            IF( kn_limiter == 1 )   CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     788         ELSE 
     789            IF( kn_limiter == 1 )   CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     790         ENDIF 
    249791          
    250       ! antidiffusive flux : high order minus low order 
    251       ! -------------------------------------------------- 
    252       DO jl = 1, ipl 
    253          DO jj = 2, jpjm1 
    254             DO ji = 1, fs_jpim1   ! vector opt. 
    255                zfu_ho(ji,jj,jl) = zfu_ho(ji,jj,jl) - zfu_ups(ji,jj,jl) 
    256             END DO 
    257          END DO 
    258       END DO 
    259       DO jl = 1, ipl 
    260          DO jj = 1, jpjm1 
    261             DO ji = fs_2, fs_jpim1   ! vector opt. 
    262                zfv_ho(ji,jj,jl) = zfv_ho(ji,jj,jl) - zfv_ups(ji,jj,jl) 
    263             END DO 
    264          END DO 
    265       END DO 
    266  
    267       CALL lbc_lnk("icedyn_adv_umx",zt_ups, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign) 
    268        
    269       ! monotonicity algorithm 
    270       ! ------------------------- 
    271       CALL nonosc_2d( ipl, ptc, zfu_ho, zfv_ho, zt_ups, pdt ) 
    272        
    273       ! final trend with corrected fluxes 
    274       ! ------------------------------------ 
    275       DO jl = 1, ipl 
    276          DO jj = 2, jpjm1 
    277             DO ji = fs_2, fs_jpim1   ! vector opt.   
    278                ztra       = ztrd(ji,jj,jl)  - (  zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj  ,jl)   & 
    279                   &                         +    zfv_ho(ji,jj,jl) - zfv_ho(ji  ,jj-1,jl) ) * r1_e1e2t(ji,jj)   
    280                ptc(ji,jj,jl) = ptc(ji,jj,jl) + pdt * ztra * tmask(ji,jj,1) 
    281             END DO 
    282          END DO 
    283       END DO 
    284       ! 
    285    END SUBROUTINE adv_umx 
    286  
    287    SUBROUTINE macho( k_order, kt, ipl, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v ) 
     792      ENDIF 
     793    
     794   END SUBROUTINE cen2 
     795 
     796    
     797   SUBROUTINE macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, pt_u, pt_v, pfu_ho, pfv_ho, & 
     798      &              pt_ups, pfu_ups, pfv_ups ) 
     799      !!--------------------------------------------------------------------- 
     800      !!                    ***  ROUTINE macho  *** 
     801      !!      
     802      !! **  Purpose :   compute   
     803      !! 
     804      !! **  Method  :   ... ??? 
     805      !!                 TIM = transient interpolation Modeling  
     806      !! 
     807      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
     808      !!---------------------------------------------------------------------- 
     809      REAL(wp)                    , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0) 
     810      INTEGER                     , INTENT(in   ) ::   kn_limiter       ! limiter 
     811      INTEGER                     , INTENT(in   ) ::   kn_umx           ! order of the scheme (1-5=UM or 20=CEN2) 
     812      INTEGER                     , INTENT(in   ) ::   jt               ! number of sub-iteration 
     813      INTEGER                     , INTENT(in   ) ::   kt               ! number of iteration 
     814      REAL(wp)                    , INTENT(in   ) ::   pdt              ! tracer time-step 
     815      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt               ! tracer fields 
     816      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pu, pv           ! 2 ice velocity components 
     817      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components 
     818      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pubox, pvbox     ! upstream velocity 
     819      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptc              ! tracer content at before time step  
     820      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_u, pt_v       ! tracer at u- and v-points  
     821      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho, pfv_ho   ! high order fluxes  
     822      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content  
     823      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes  
     824      ! 
     825      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
     826      REAL(wp) ::   ztra 
     827      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zzt, zzfu_ho, zzfv_ho 
     828      !!---------------------------------------------------------------------- 
     829      ! 
     830      IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     831         ! 
     832         !                                                        !--  ultimate interpolation of pt at u-point  --! 
     833         CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) 
     834         !                                                        !--  limiter in x --! 
     835         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 
     836         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
     837         !                                                        !--  advective form update in zzt  --! 
     838 
     839         IF( ll_1stguess_clem ) THEN 
     840 
     841            ! first guess of tracer content from u-flux 
     842            DO jl = 1, jpl 
     843               DO jj = 2, jpjm1 
     844                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     845                     IF( ll_clem ) THEN 
     846                        IF( ll_gurvan ) THEN 
     847                           zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     848                        ELSE 
     849                           zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     850                              &            + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     851                              &         ) * tmask(ji,jj,1) 
     852                        ENDIF 
     853                     ELSE 
     854                        zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     855                     ENDIF 
     856                  END DO 
     857               END DO 
     858            END DO 
     859            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
     860 
     861         ELSE 
     862 
     863            DO jl = 1, jpl 
     864               DO jj = 2, jpjm1 
     865                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     866                     IF( ll_gurvan ) THEN 
     867                        zzt(ji,jj,jl) = pt(ji,jj,jl) - pubox(ji,jj   ) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj)  & 
     868                           &                         - pt   (ji,jj,jl) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj) 
     869                     ELSE 
     870                        zzt(ji,jj,jl) = pt(ji,jj,jl) - pubox(ji,jj   ) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj)  & 
     871                           &                         - pt   (ji,jj,jl) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk 
     872                     ENDIF 
     873                     zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 
     874                  END DO 
     875               END DO 
     876            END DO 
     877            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
     878         ENDIF 
     879         ! 
     880         !                                                        !--  ultimate interpolation of pt at v-point  --! 
     881         IF( ll_hoxy ) THEN 
     882            CALL ultimate_y( kn_umx, pdt, zzt, pv, pvc, pt_v, pfv_ho ) 
     883         ELSE 
     884            CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) 
     885         ENDIF 
     886         !                                                        !--  limiter in y --! 
     887         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 
     888         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     889         !          
     890         ! 
     891      ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
     892         ! 
     893         !                                                        !--  ultimate interpolation of pt at v-point  --! 
     894         CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) 
     895         !                                                        !--  limiter in y --! 
     896         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 
     897         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     898         !                                                        !--  advective form update in zzt  --! 
     899         IF( ll_1stguess_clem ) THEN 
     900             
     901            ! first guess of tracer content from v-flux  
     902            DO jl = 1, jpl 
     903               DO jj = 2, jpjm1 
     904                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     905                     IF( ll_clem ) THEN 
     906                        IF( ll_gurvan ) THEN 
     907                           zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     908                        ELSE 
     909                           zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 
     910                              &            + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     911                              &         ) * tmask(ji,jj,1) 
     912                        ENDIF 
     913                     ELSE 
     914                        zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) & 
     915                           &         * tmask(ji,jj,1) 
     916                     ENDIF 
     917                  END DO 
     918               END DO 
     919            END DO 
     920            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
     921             
     922         ELSE 
     923             
     924            DO jl = 1, jpl 
     925               DO jj = 2, jpjm1 
     926                  DO ji = fs_2, fs_jpim1 
     927                     IF( ll_gurvan ) THEN 
     928                        zzt(ji,jj,jl) = pt(ji,jj,jl) - pvbox(ji,jj   ) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj)  & 
     929                           &                         - pt   (ji,jj,jl) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) 
     930                     ELSE 
     931                        zzt(ji,jj,jl) = pt(ji,jj,jl) - pvbox(ji,jj   ) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj)  & 
     932                           &                         - pt   (ji,jj,jl) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk 
     933                     ENDIF 
     934                     zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 
     935                  END DO 
     936               END DO 
     937            END DO 
     938            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
     939         ENDIF 
     940         ! 
     941         !                                                        !--  ultimate interpolation of pt at u-point  --! 
     942         IF( ll_hoxy ) THEN 
     943            CALL ultimate_x( kn_umx, pdt, zzt, pu, puc, pt_u, pfu_ho ) 
     944         ELSE 
     945            CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) 
     946         ENDIF 
     947         !                                                        !--  limiter in x --! 
     948         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 
     949         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
     950         ! 
     951         ! 
     952      ENDIF 
     953 
     954      
     955      IF( kn_limiter == 1 ) THEN 
     956         IF( .NOT. ll_limiter_it2 ) THEN 
     957            IF( ll_clem ) THEN 
     958               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     959            ELSE 
     960               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     961            ENDIF 
     962         ELSE 
     963            zzfu_ho(:,:,:) = pfu_ho(:,:,:) 
     964            zzfv_ho(:,:,:) = pfv_ho(:,:,:) 
     965            ! 1st iteration of nonosc (limit the flux with the upstream solution) 
     966            IF( ll_clem ) THEN 
     967               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho ) 
     968            ELSE 
     969               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho ) 
     970            ENDIF 
     971            ! guess after content field with high order 
     972            DO jl = 1, jpl 
     973               DO jj = 2, jpjm1 
     974                  DO ji = fs_2, fs_jpim1 
     975                     ztra          = - ( zzfu_ho(ji,jj,jl) - zzfu_ho(ji-1,jj,jl) + zzfv_ho(ji,jj,jl) - zzfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj) 
     976                     zzt(ji,jj,jl) =   ( ptc(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1) 
     977                  END DO 
     978               END DO 
     979            END DO 
     980            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
     981            ! 2nd iteration of nonosc (limit the flux with the limited high order solution) 
     982            IF( ll_clem ) THEN 
     983               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho ) 
     984            ELSE 
     985               CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho ) 
     986            ENDIF 
     987         ENDIF 
     988      ENDIF 
     989      ! 
     990   END SUBROUTINE macho 
     991 
     992 
     993   SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) 
    288994      !!--------------------------------------------------------------------- 
    289995      !!                    ***  ROUTINE ultimate_x  *** 
    … …  
    2961002      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    2971003      !!---------------------------------------------------------------------- 
    298       INTEGER                         , INTENT(in   ) ::   k_order    ! order of the ULTIMATE scheme 
    299       INTEGER                         , INTENT(in   ) ::   kt         ! number of iteration 
    300       INTEGER                         , INTENT(in   ) ::   ipl        ! third dimension of tracer array 
    301       REAL(wp)                        , INTENT(in   ) ::   pdt        ! tracer time-step 
    302       REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(in   ) ::   ptc        ! tracer fields 
    303       REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   puc, pvc   ! 2 ice velocity components 
    304       REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   pubox, pvbox   ! upstream velocity 
    305       REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(  out) ::   pt_u, pt_v ! tracer at u- and v-points  
    306       ! 
    307       INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    308       REAL(wp) ::   zc_box    !   -      - 
    309       REAL(wp), DIMENSION(jpi,jpj,ipl) :: zzt 
    310       !!---------------------------------------------------------------------- 
    311       ! 
    312       IF( MOD( (kt - 1) / nn_fsbc , 2 ) == 0 ) THEN         !==  odd ice time step:  adv_x then adv_y  ==! 
    313          ! 
    314          !                                                           !--  ultimate interpolation of pt at u-point  --! 
    315          CALL ultimate_x( k_order, ipl, pdt, ptc, puc, pt_u ) 
    316          ! 
    317          !                                                           !--  advective form update in zzt  --! 
    318          DO jl = 1, ipl 
    319             DO jj = 2, jpjm1 
    320                DO ji = fs_2, fs_jpim1   ! vector opt. 
    321                   zzt(ji,jj,jl) = ptc(ji,jj,jl) - pubox(ji,jj) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj)  & 
    322                    &            - ptc(ji,jj,jl) * pdt * ( puc (ji,jj) - puc (ji-1,jj) ) * r1_e1e2t(ji,jj) 
    323                   zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 
    324                END DO 
    325             END DO 
    326          END DO 
    327          CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
    328          ! 
    329          !                                                           !--  ultimate interpolation of pt at v-point  --! 
    330          CALL ultimate_y( k_order, ipl, pdt, zzt, pvc, pt_v ) 
    331          ! 
    332       ELSE                                                  !==  even ice time step:  adv_y then adv_x  ==! 
    333          ! 
    334          !                                                           !--  ultimate interpolation of pt at v-point  --! 
    335          CALL ultimate_y( k_order, ipl, pdt, ptc, pvc, pt_v ) 
    336          ! 
    337          !                                                           !--  advective form update in zzt  --! 
    338          DO jl = 1, ipl 
    339             DO jj = 2, jpjm1 
    340                DO ji = fs_2, fs_jpim1 
    341                   zzt(ji,jj,jl) = ptc(ji,jj,jl) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj)  & 
    342                      &                    - ptc  (ji,jj,jl) * pdt * ( pvc (ji,jj) - pvc (ji,jj-1) ) * r1_e1e2t(ji,jj) 
    343                   zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 
    344                END DO 
    345             END DO 
    346          END DO 
    347          CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
    348          ! 
    349          !                                                           !--  ultimate interpolation of pt at u-point  --! 
    350          CALL ultimate_x( k_order, ipl, pdt, zzt, puc, pt_u ) 
    351          !       
    352       ENDIF       
    353       ! 
    354    END SUBROUTINE macho 
    355  
    356  
    357    SUBROUTINE ultimate_x( k_order, ipl, pdt, pt, puc, pt_u ) 
    358       !!--------------------------------------------------------------------- 
    359       !!                    ***  ROUTINE ultimate_x  *** 
    360       !!      
    361       !! **  Purpose :   compute   
    362       !! 
    363       !! **  Method  :   ... ??? 
    364       !!                 TIM = transient interpolation Modeling  
    365       !! 
    366       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    367       !!---------------------------------------------------------------------- 
    368       INTEGER                         , INTENT(in   ) ::   k_order   ! ocean time-step index 
    369       INTEGER                         , INTENT(in   ) ::   ipl        ! third dimension of tracer array 
    370       REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step 
    371       REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   puc       ! ice i-velocity component 
    372       REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(in   ) ::   pt        ! tracer fields 
    373       REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(  out) ::   pt_u      ! tracer at u-point  
    374       ! 
    375       INTEGER  ::   ji, jj, jl       ! dummy loop indices 
     1004      INTEGER                     , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2) 
     1005      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step 
     1006      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pu        ! ice i-velocity component 
     1007      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   puc       ! ice i-velocity * A component 
     1008      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt        ! tracer fields 
     1009      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_u      ! tracer at u-point  
     1010      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho    ! high order flux  
     1011      ! 
     1012      INTEGER  ::   ji, jj, jl             ! dummy loop indices 
    3761013      REAL(wp) ::   zcu, zdx2, zdx4    !   -      - 
    377       REAL(wp), DIMENSION(jpi,jpj) :: ztu1, ztu3 
    378       REAL(wp), DIMENSION(jpi,jpj,ipl) :: ztu2, ztu4 
     1014      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztu1, ztu2, ztu3, ztu4 
    3791015      !!---------------------------------------------------------------------- 
    3801016      ! 
    3811017      !                                                     !--  Laplacian in i-direction  --! 
    382       DO jl = 1, ipl 
     1018      DO jl = 1, jpl 
    3831019         DO jj = 2, jpjm1         ! First derivative (gradient) 
    3841020            DO ji = 1, fs_jpim1 
    385                ztu1(ji,jj) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 
     1021               ztu1(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 
    3861022            END DO 
    3871023            !                     ! Second derivative (Laplacian) 
    3881024            DO ji = fs_2, fs_jpim1 
    389                ztu2(ji,jj,jl) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj) 
     1025               ztu2(ji,jj,jl) = ( ztu1(ji,jj,jl) - ztu1(ji-1,jj,jl) ) * r1_e1t(ji,jj) 
    3901026            END DO 
    3911027         END DO 
    … …  
    3941030      ! 
    3951031      !                                                     !--  BiLaplacian in i-direction  --! 
    396       DO jl = 1, ipl 
     1032      DO jl = 1, jpl 
    3971033         DO jj = 2, jpjm1         ! Third derivative 
    3981034            DO ji = 1, fs_jpim1 
    399                ztu3(ji,jj) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 
     1035               ztu3(ji,jj,jl) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 
    4001036            END DO 
    4011037            !                     ! Fourth derivative 
    4021038            DO ji = fs_2, fs_jpim1 
    403                ztu4(ji,jj,jl) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj) 
     1039               ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj) 
    4041040            END DO 
    4051041         END DO 
    … …  
    4081044      ! 
    4091045      ! 
    410       SELECT CASE (k_order ) 
     1046      SELECT CASE (kn_umx ) 
    4111047      ! 
    4121048      CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21) 
    4131049         !         
    414          DO jl = 1, ipl 
    415             DO jj = 2, jpjm1 
     1050         DO jl = 1, jpl 
     1051            DO jj = 1, jpjm1 
    4161052               DO ji = 1, fs_jpim1   ! vector opt. 
    417                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                               pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    418                      &                                       - SIGN( 1._wp, puc(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
     1053                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                           pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
     1054                     &                                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
    4191055               END DO 
    4201056            END DO 
    … …  
    4231059      CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23) 
    4241060         ! 
    425          DO jl = 1, ipl 
    426             DO jj = 2, jpjm1 
     1061         DO jl = 1, jpl 
     1062            DO jj = 1, jpjm1 
    4271063               DO ji = 1, fs_jpim1   ! vector opt. 
    428                   zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    429                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                               pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    430                      &                                              -              zcu   * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )  
     1064                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     1065                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
     1066                     &                                               -              zcu   * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )  
    4311067               END DO 
    4321068            END DO 
    … …  
    4351071      CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24) 
    4361072         ! 
    437          DO jl = 1, ipl 
    438             DO jj = 2, jpjm1 
     1073         DO jl = 1, jpl 
     1074            DO jj = 1, jpjm1 
    4391075               DO ji = 1, fs_jpim1   ! vector opt. 
    440                   zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     1076                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    4411077                  zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    4421078!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    443                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                     pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        & 
    444                      &                                              -              zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   & 
    445                      &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                        ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        & 
    446                      &                                              - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )   ) 
     1079                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        & 
     1080                     &                                               -              zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   & 
     1081                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                         ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        & 
     1082                     &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )   ) 
    4471083               END DO 
    4481084            END DO 
    … …  
    4511087      CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27) 
    4521088         ! 
    453          DO jl = 1, ipl 
    454             DO jj = 2, jpjm1 
     1089         DO jl = 1, jpl 
     1090            DO jj = 1, jpjm1 
    4551091               DO ji = 1, fs_jpim1   ! vector opt. 
    456                   zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     1092                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    4571093                  zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    4581094!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    459                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                     pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        & 
    460                      &                                                    -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   & 
    461                      &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                        ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        & 
    462                      &                                                    - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )   ) 
     1095                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        & 
     1096                     &                                               -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   & 
     1097                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                   ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        & 
     1098                     &                                               - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )   ) 
    4631099               END DO 
    4641100            END DO 
    … …  
    4671103      CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29) 
    4681104         ! 
    469          DO jl = 1, ipl 
    470             DO jj = 2, jpjm1 
     1105         DO jl = 1, jpl 
     1106            DO jj = 1, jpjm1 
    4711107               DO ji = 1, fs_jpim1   ! vector opt. 
    472                   zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     1108                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    4731109                  zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    474 !!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
     1110                  !!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    4751111                  zdx4 = zdx2 * zdx2 
    476                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (           (                   pt  (ji+1,jj,jl) + pt  (ji,jj,jl)       & 
    477                      &                                                    -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) )   & 
    478                      &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *    (                   ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)       & 
    479                      &                                                    - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) )   & 
    480                      &       + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl)       & 
    481                      &                                              - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) ) 
     1112                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (            (                   pt  (ji+1,jj,jl) + pt  (ji,jj,jl)       & 
     1113                     &                                                     -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) )   & 
     1114                     &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *     (                   ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)       & 
     1115                     &                                                     - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) )   & 
     1116                     &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl)       & 
     1117                     &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) ) 
    4821118               END DO 
    4831119            END DO 
    … …  
    4851121         ! 
    4861122      END SELECT 
     1123      !                                                     !-- High order flux in i-direction  --! 
     1124      IF( ll_neg ) THEN 
     1125         DO jl = 1, jpl 
     1126            DO jj = 1, jpjm1 
     1127               DO ji = 1, fs_jpim1 
     1128                  IF( pt_u(ji,jj,jl) < 0._wp ) THEN 
     1129                     pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                           pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
     1130                        &                                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
     1131                  ENDIF 
     1132               END DO 
     1133            END DO 
     1134         END DO 
     1135      ENDIF 
     1136 
     1137      DO jl = 1, jpl 
     1138         DO jj = 1, jpjm1 
     1139            DO ji = 1, fs_jpim1   ! vector opt. 
     1140               IF( ll_clem ) THEN 
     1141                  pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 
     1142               ELSE 
     1143                  pfu_ho(ji,jj,jl) = puc(ji,jj,jl) * pt_u(ji,jj,jl) 
     1144               ENDIF 
     1145            END DO 
     1146         END DO 
     1147      END DO 
    4871148      ! 
    4881149   END SUBROUTINE ultimate_x 
    4891150    
    4901151  
    491    SUBROUTINE ultimate_y( k_order, ipl, pdt, pt, pvc, pt_v ) 
     1152   SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) 
    4921153      !!--------------------------------------------------------------------- 
    4931154      !!                    ***  ROUTINE ultimate_y  *** 
    … …  
    5001161      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    5011162      !!---------------------------------------------------------------------- 
    502       INTEGER                         , INTENT(in   ) ::   k_order   ! ocean time-step index 
    503       INTEGER                         , INTENT(in   ) ::   ipl       ! third dimension of tracer array 
    504       REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step 
    505       REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   pvc       ! ice j-velocity component 
    506       REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(in   ) ::   pt        ! tracer fields 
    507       REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(  out) ::   pt_v      ! tracer at v-point  
     1163      INTEGER                     , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2) 
     1164      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step 
     1165      REAL(wp), DIMENSION(:,:  ), INTENT(in   ) ::   pv        ! ice j-velocity component 
     1166      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pvc       ! ice j-velocity*A component 
     1167      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt        ! tracer fields 
     1168      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_v      ! tracer at v-point  
     1169      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfv_ho    ! high order flux  
    5081170      ! 
    5091171      INTEGER  ::   ji, jj, jl       ! dummy loop indices 
    5101172      REAL(wp) ::   zcv, zdy2, zdy4    !   -      - 
    511       REAL(wp), DIMENSION(jpi,jpj) :: ztv1, ztv3 
    512       REAL(wp), DIMENSION(jpi,jpj,ipl) :: ztv2, ztv4 
     1173      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztv1, ztv2, ztv3, ztv4 
    5131174      !!---------------------------------------------------------------------- 
    5141175      ! 
    5151176      !                                                     !--  Laplacian in j-direction  --! 
    516       DO jl = 1, ipl 
     1177      DO jl = 1, jpl 
    5171178         DO jj = 1, jpjm1         ! First derivative (gradient) 
    5181179            DO ji = fs_2, fs_jpim1 
    519                ztv1(ji,jj) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
     1180               ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
    5201181            END DO 
    5211182         END DO 
    5221183         DO jj = 2, jpjm1         ! Second derivative (Laplacian) 
    5231184            DO ji = fs_2, fs_jpim1 
    524                ztv2(ji,jj,jl) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj) 
     1185               ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj) 
    5251186            END DO 
    5261187         END DO 
    … …  
    5291190      ! 
    5301191      !                                                     !--  BiLaplacian in j-direction  --! 
    531       DO jl = 1, ipl 
     1192      DO jl = 1, jpl 
    5321193         DO jj = 1, jpjm1         ! First derivative 
    5331194            DO ji = fs_2, fs_jpim1 
    534             ztv3(ji,jj) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
     1195               ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
    5351196            END DO 
    5361197         END DO 
    5371198         DO jj = 2, jpjm1         ! Second derivative 
    5381199            DO ji = fs_2, fs_jpim1 
    539                ztv4(ji,jj,jl) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj) 
     1200               ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj) 
    5401201            END DO 
    5411202         END DO 
    … …  
    5441205      ! 
    5451206      ! 
    546       SELECT CASE (k_order ) 
    547       ! 
     1207      SELECT CASE (kn_umx ) 
     1208         ! 
    5481209      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21) 
    549          DO jl = 1, ipl 
     1210         DO jl = 1, jpl 
    5501211            DO jj = 1, jpjm1 
    551                DO ji = fs_2, fs_jpim1 
    552                   pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                           ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
    553                      &                                     - SIGN( 1._wp, pvc(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
     1212               DO ji = 1, fs_jpim1 
     1213                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                          ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
     1214                     &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
    5541215               END DO 
    5551216            END DO 
    … …  
    5571218         ! 
    5581219      CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23) 
    559          DO jl = 1, ipl 
     1220         DO jl = 1, jpl 
    5601221            DO jj = 1, jpjm1 
    561                DO ji = fs_2, fs_jpim1 
    562                   zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     1222               DO ji = 1, fs_jpim1 
     1223                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    5631224                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (     ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
    5641225                     &                                     - zcv * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
    … …  
    5691230         ! 
    5701231      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24) 
    571          DO jl = 1, ipl 
     1232         DO jl = 1, jpl 
    5721233            DO jj = 1, jpjm1 
    573                DO ji = fs_2, fs_jpim1 
    574                   zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     1234               DO ji = 1, fs_jpim1 
     1235                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    5751236                  zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    5761237!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    … …  
    5841245         ! 
    5851246      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27) 
    586          DO jl = 1, ipl 
     1247         DO jl = 1, jpl 
    5871248            DO jj = 1, jpjm1 
    588                DO ji = fs_2, fs_jpim1 
    589                   zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     1249               DO ji = 1, fs_jpim1 
     1250                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    5901251                  zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    5911252!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    … …  
    5991260         ! 
    6001261      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29) 
    601          DO jl = 1, ipl 
     1262         DO jl = 1, jpl 
    6021263            DO jj = 1, jpjm1 
    603                DO ji = fs_2, fs_jpim1 
    604                   zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     1264               DO ji = 1, fs_jpim1 
     1265                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    6051266                  zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    6061267!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    … …  
    6171278         ! 
    6181279      END SELECT 
     1280      !                                                     !-- High order flux in j-direction  --! 
     1281      IF( ll_neg ) THEN 
     1282         DO jl = 1, jpl 
     1283            DO jj = 1, jpjm1 
     1284               DO ji = 1, fs_jpim1 
     1285                  IF( pt_v(ji,jj,jl) < 0._wp ) THEN 
     1286                     pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                          ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
     1287                        &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
     1288                  ENDIF 
     1289               END DO 
     1290            END DO 
     1291         END DO 
     1292      ENDIF 
     1293 
     1294      DO jl = 1, jpl 
     1295         DO jj = 1, jpjm1 
     1296            DO ji = 1, fs_jpim1   ! vector opt. 
     1297               IF( ll_clem ) THEN 
     1298                  pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 
     1299               ELSE 
     1300                  pfv_ho(ji,jj,jl) = pvc(ji,jj,jl) * pt_v(ji,jj,jl) 
     1301               ENDIF 
     1302            END DO 
     1303         END DO 
     1304      END DO 
    6191305      ! 
    6201306   END SUBROUTINE ultimate_y 
    621     
    622    
    623    SUBROUTINE nonosc_2d( ipl, pbef, paa, pbb, paft, pdt ) 
     1307      
     1308 
     1309   SUBROUTINE nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho ) 
    6241310      !!--------------------------------------------------------------------- 
    6251311      !!                    ***  ROUTINE nonosc  *** 
    6261312      !!      
    627       !! **  Purpose :   compute monotonic tracer fluxes from the upstream  
     1313       !! **  Purpose :   compute monotonic tracer fluxes from the upstream  
    6281314      !!       scheme and the before field by a nonoscillatory algorithm  
    6291315      !! 
    6301316      !! **  Method  :   ... ??? 
    631       !!       warning : pbef and paft must be masked, but the boundaries 
     1317      !!       warning : pt and pt_low must be masked, but the boundaries 
    6321318      !!       conditions on the fluxes are not necessary zalezak (1979) 
    6331319      !!       drange (1995) multi-dimensional forward-in-time and upstream- 
    6341320      !!       in-space based differencing for fluid 
    6351321      !!---------------------------------------------------------------------- 
    636       INTEGER                          , INTENT(in   ) ::   ipl          ! third dimension of tracer array 
    637       REAL(wp)                         , INTENT(in   ) ::   pdt          ! tracer time-step 
    638       REAL(wp), DIMENSION (jpi,jpj,ipl), INTENT(in   ) ::   pbef, paft   ! before & after field 
    639       REAL(wp), DIMENSION (jpi,jpj,ipl), INTENT(inout) ::   paa, pbb     ! monotonic fluxes in the 2 directions 
     1322      REAL(wp)                     , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0) 
     1323      REAL(wp)                     , INTENT(in   ) ::   pdt              ! tracer time-step 
     1324      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pu               ! ice i-velocity => u*e2 
     1325      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   puc              ! ice i-velocity *A => u*e2*a 
     1326      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pv               ! ice j-velocity => v*e1 
     1327      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pvc              ! ice j-velocity *A => v*e1*a 
     1328      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   ptc, pt, pt_low  ! before field & upstream guess of after field 
     1329      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   pfv_low, pfu_low ! upstream flux 
     1330      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux 
    6401331      ! 
    6411332      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    642       INTEGER  ::   ikm1      ! local integer 
    643       REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zsml, z1_dt   ! local scalars 
    644       REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      - 
    645       REAL(wp), DIMENSION(jpi,jpj) :: zbup, zbdo, zmsk 
    646       REAL(wp), DIMENSION(jpi,jpj,ipl) :: zbetup, zbetdo, zdiv 
    647       !!---------------------------------------------------------------------- 
    648       ! 
     1333      REAL(wp) ::   zpos, zneg, zbig, zsml, z1_dt, zpos2, zneg2   ! local scalars 
     1334      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign, zcoef        !   -      - 
     1335      REAL(wp), DIMENSION(jpi,jpj    ) :: zbup, zbdo 
     1336      REAL(wp), DIMENSION(jpi,jpj,jpl) :: zbetup, zbetdo, zti_low, ztj_low, zzt 
     1337      !!---------------------------------------------------------------------- 
    6491338      zbig = 1.e+40_wp 
    650       zsml = 1.e-15_wp 
    651  
    652       ! test on divergence 
    653       DO jl = 1, ipl 
     1339      zsml = epsi20 
     1340       
     1341      IF( ll_zeroup2 ) THEN 
     1342         DO jl = 1, jpl 
     1343            DO jj = 1, jpjm1 
     1344               DO ji = 1, fs_jpim1   ! vector opt. 
     1345                  IF( amaxu(ji,jj,jl) == 0._wp )   pfu_ho(ji,jj,jl) = 0._wp 
     1346                  IF( amaxv(ji,jj,jl) == 0._wp )   pfv_ho(ji,jj,jl) = 0._wp 
     1347               END DO 
     1348            END DO 
     1349         END DO 
     1350      ENDIF 
     1351 
     1352      IF( ll_zeroup4 ) THEN 
     1353         DO jl = 1, jpl 
     1354            DO jj = 1, jpjm1 
     1355               DO ji = 1, fs_jpim1   ! vector opt. 
     1356                  IF( pfu_low(ji,jj,jl) == 0._wp )   pfu_ho(ji,jj,jl) = 0._wp 
     1357                  IF( pfv_low(ji,jj,jl) == 0._wp )   pfv_ho(ji,jj,jl) = 0._wp 
     1358               END DO 
     1359            END DO 
     1360         END DO 
     1361      ENDIF 
     1362 
     1363 
     1364      IF( ll_zeroup1 ) THEN 
     1365         DO jl = 1, jpl 
     1366            DO jj = 2, jpjm1 
     1367               DO ji = fs_2, fs_jpim1 
     1368                  IF( ll_gurvan ) THEN 
     1369                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1370                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1371                  ELSE 
     1372                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1373                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1374                        &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1375                        &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1376                        &         ) * tmask(ji,jj,1) 
     1377                  ENDIF 
     1378                  IF( zzt(ji,jj,jl) < 0._wp ) THEN 
     1379                     pfu_ho(ji,jj,jl)   = pfu_low(ji,jj,jl) 
     1380                     pfv_ho(ji,jj,jl)   = pfv_low(ji,jj,jl) 
     1381                     WRITE(numout,*) '*** 1 negative high order zzt ***',ji,jj,zzt(ji,jj,jl) 
     1382                  ENDIF 
     1383                  !!               IF( ji==26 .AND. jj==86) THEN 
     1384                  !!                  WRITE(numout,*) 'zzt high order',zzt(ji,jj) 
     1385                  !!                  WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt 
     1386                  !!                  WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt 
     1387                  !!                  WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj,jl)) * r1_e1e2t(ji,jj) * pdt 
     1388                  !!                  WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1,jl)) * r1_e1e2t(ji,jj) * pdt 
     1389                  !!               ENDIF 
     1390                  IF( ll_gurvan ) THEN 
     1391                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1392                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1393                  ELSE 
     1394                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1395                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1396                        &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1397                        &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1398                        &         ) * tmask(ji,jj,1) 
     1399                  ENDIF 
     1400                  IF( zzt(ji,jj,jl) < 0._wp ) THEN 
     1401                     pfu_ho(ji-1,jj,jl) = pfu_low(ji-1,jj,jl) 
     1402                     pfv_ho(ji,jj-1,jl) = pfv_low(ji,jj-1,jl) 
     1403                     WRITE(numout,*) '*** 2 negative high order zzt ***',ji,jj,zzt(ji,jj,jl) 
     1404                  ENDIF 
     1405                  IF( ll_gurvan ) THEN 
     1406                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1407                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1408                  ELSE 
     1409                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1410                        &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1411                        &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1412                        &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1413                        &         ) * tmask(ji,jj,1) 
     1414                  ENDIF 
     1415                  IF( zzt(ji,jj,jl) < 0._wp ) THEN 
     1416                     WRITE(numout,*) '*** 3 negative high order zzt ***',ji,jj,zzt(ji,jj,jl) 
     1417                  ENDIF 
     1418               END DO 
     1419            END DO 
     1420         END DO 
     1421         CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) 
     1422      ENDIF 
     1423 
     1424       
     1425      ! antidiffusive flux : high order minus low order 
     1426      ! -------------------------------------------------- 
     1427      DO jl = 1, jpl 
     1428         DO jj = 1, jpjm1 
     1429            DO ji = 1, fs_jpim1   ! vector opt. 
     1430               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_low(ji,jj,jl) 
     1431               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_low(ji,jj,jl) 
     1432            END DO 
     1433         END DO 
     1434      END DO 
     1435 
     1436      ! extreme case where pfu_ho has to be zero 
     1437      ! ---------------------------------------- 
     1438      !                                    pfu_ho 
     1439      !                           *         ---> 
     1440      !                        |      |  *   |        |  
     1441      !                        |      |      |    *   |     
     1442      !                        |      |      |        |    * 
     1443      !            t_low :       i-1     i       i+1       i+2    
     1444      IF( ll_prelimiter_zalesak ) THEN 
     1445          
     1446         DO jl = 1, jpl 
     1447            DO jj = 2, jpjm1 
     1448               DO ji = fs_2, fs_jpim1  
     1449                  zti_low(ji,jj,jl)= pt_low(ji+1,jj  ,jl) 
     1450                  ztj_low(ji,jj,jl)= pt_low(ji  ,jj+1,jl) 
     1451               END DO 
     1452            END DO 
     1453         END DO 
     1454         CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_low, 'T', 1., ztj_low, 'T', 1. ) 
     1455 
     1456         !! this does not work ?? 
     1457         !!            DO jj = 2, jpjm1 
     1458         !!               DO ji = fs_2, fs_jpim1 
     1459         !!                  IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji+1,jj  ) - pt_low (ji  ,jj) ) .AND.     & 
     1460         !!               &      SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj+1) - pt_low (ji  ,jj) )           & 
     1461         !!               &    ) THEN 
     1462         !!                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_low(ji+1,jj ) - zti_low(ji  ,jj) ) .AND.   & 
     1463         !!               &         SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_low(ji,jj+1 ) - ztj_low(ji  ,jj) )         & 
     1464         !!               &       ) THEN 
     1465         !!                        pfu_ho(ji,jj) = 0.  ;   pfv_ho(ji,jj) = 0. 
     1466         !!                     ENDIF 
     1467         !!                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj) - pt_low (ji-1,jj  ) ) .AND.  & 
     1468         !!               &         SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj) - pt_low (ji  ,jj-1) )        & 
     1469         !!               &       )  THEN 
     1470         !!                        pfu_ho(ji,jj) = 0.  ;   pfv_ho(ji,jj) = 0. 
     1471         !!                     ENDIF 
     1472         !!                  ENDIF 
     1473         !!                END DO 
     1474         !!            END DO             
     1475 
     1476         DO jl = 1, jpl 
     1477            DO jj = 2, jpjm1 
     1478               DO ji = fs_2, fs_jpim1 
     1479                  IF ( pfu_ho(ji,jj,jl) * ( pt_low(ji+1,jj,jl) - pt_low(ji,jj,jl) ) <= 0. .AND.  & 
     1480                     & pfv_ho(ji,jj,jl) * ( pt_low(ji,jj+1,jl) - pt_low(ji,jj,jl) ) <= 0. ) THEN 
     1481                     ! 
     1482                     IF(  pfu_ho(ji,jj,jl) * ( zti_low(ji+1,jj,jl) - zti_low(ji,jj,jl) ) <= 0 .AND.  & 
     1483                        & pfv_ho(ji,jj,jl) * ( ztj_low(ji,jj+1,jl) - ztj_low(ji,jj,jl) ) <= 0)  pfu_ho(ji,jj,jl)=0. ; pfv_ho(ji,jj,jl)=0. 
     1484                     ! 
     1485                     IF(  pfu_ho(ji,jj,jl) * ( pt_low(ji  ,jj,jl) - pt_low(ji-1,jj,jl) ) <= 0 .AND.  & 
     1486                        & pfv_ho(ji,jj,jl) * ( pt_low(ji  ,jj,jl) - pt_low(ji,jj-1,jl) ) <= 0)  pfu_ho(ji,jj,jl)=0. ; pfv_ho(ji,jj,jl)=0. 
     1487                     ! 
     1488                  ENDIF 
     1489               END DO 
     1490            END DO 
     1491         END DO 
     1492         CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond. 
     1493 
     1494      ELSEIF( ll_prelimiter_devore ) THEN 
     1495         DO jl = 1, jpl 
     1496            DO jj = 2, jpjm1 
     1497               DO ji = fs_2, fs_jpim1  
     1498                  zti_low(ji,jj,jl)= pt_low(ji+1,jj  ,jl) 
     1499                  ztj_low(ji,jj,jl)= pt_low(ji  ,jj+1,jl) 
     1500               END DO 
     1501            END DO 
     1502         END DO 
     1503         CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_low, 'T', 1., ztj_low, 'T', 1. ) 
     1504 
     1505         z1_dt = 1._wp / pdt 
     1506         DO jl = 1, jpl 
     1507            DO jj = 2, jpjm1 
     1508               DO ji = fs_2, fs_jpim1 
     1509                  zsign = SIGN( 1., pt_low(ji+1,jj,jl) - pt_low(ji,jj,jl) ) 
     1510                  pfu_ho(ji,jj,jl) =  zsign * MAX( 0. , MIN( ABS(pfu_ho(ji,jj,jl)) , & 
     1511                     &                          zsign * ( pt_low (ji  ,jj,jl) - pt_low (ji-1,jj,jl) ) * e1e2t(ji  ,jj) * z1_dt , & 
     1512                     &                          zsign * ( zti_low(ji+1,jj,jl) - zti_low(ji  ,jj,jl) ) * e1e2t(ji+1,jj) * z1_dt ) ) 
     1513 
     1514                  zsign = SIGN( 1., pt_low(ji,jj+1,jl) - pt_low(ji,jj,jl) ) 
     1515                  pfv_ho(ji,jj,jl) =  zsign * MAX( 0. , MIN( ABS(pfv_ho(ji,jj,jl)) , & 
     1516                     &                          zsign * ( pt_low (ji,jj  ,jl) - pt_low (ji,jj-1,jl) ) * e1e2t(ji,jj  ) * z1_dt , & 
     1517                     &                          zsign * ( ztj_low(ji,jj+1,jl) - ztj_low(ji,jj  ,jl) ) * e1e2t(ji,jj+1) * z1_dt ) ) 
     1518               END DO 
     1519            END DO 
     1520         END DO 
     1521         CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond. 
     1522 
     1523      ENDIF 
     1524 
     1525 
     1526      ! Search local extrema 
     1527      ! -------------------- 
     1528      ! max/min of pt & pt_low with large negative/positive value (-/+zbig) outside ice cover 
     1529      z1_dt = 1._wp / pdt 
     1530      DO jl = 1, jpl 
     1531          
     1532         DO jj = 1, jpj 
     1533            DO ji = 1, jpi 
     1534               IF    ( pt(ji,jj,jl) <= 0._wp .AND. pt_low(ji,jj,jl) <= 0._wp ) THEN 
     1535                  zbup(ji,jj) = -zbig 
     1536                  zbdo(ji,jj) =  zbig 
     1537               ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_low(ji,jj,jl) > 0._wp ) THEN 
     1538                  zbup(ji,jj) = pt_low(ji,jj,jl) 
     1539                  zbdo(ji,jj) = pt_low(ji,jj,jl) 
     1540               ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_low(ji,jj,jl) <= 0._wp ) THEN 
     1541                  zbup(ji,jj) = pt(ji,jj,jl) 
     1542                  zbdo(ji,jj) = pt(ji,jj,jl) 
     1543               ELSE 
     1544                  zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_low(ji,jj,jl) ) 
     1545                  zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_low(ji,jj,jl) ) 
     1546               ENDIF 
     1547            END DO 
     1548         END DO 
     1549 
    6541550         DO jj = 2, jpjm1 
    655             DO ji = fs_2, fs_jpim1   ! vector opt.   
    656                zdiv(ji,jj,jl) =  - (  paa(ji,jj,jl) - paa(ji-1,jj  ,jl)   & 
    657                   &              +    pbb(ji,jj,jl) - pbb(ji  ,jj-1,jl) )   
    658             END DO 
    659          END DO 
    660       END DO 
    661       CALL lbc_lnk( 'icedyn_adv_umx', zdiv, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign) 
    662  
    663       DO jl = 1, ipl 
    664          ! Determine ice masks for before and after tracers  
    665          WHERE( pbef(:,:,jl) == 0._wp .AND. paft(:,:,jl) == 0._wp .AND. zdiv(:,:,jl) == 0._wp )  
    666             zmsk(:,:) = 0._wp 
    667          ELSEWHERE                                                                                    
    668             zmsk(:,:) = 1._wp * tmask(:,:,1) 
    669          END WHERE 
    670  
    671          ! Search local extrema 
    672          ! -------------------- 
    673          ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 
    674 !         zbup(:,:) = MAX( pbef(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) ),   & 
    675 !            &             paft(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) )  ) 
    676 !         zbdo(:,:) = MIN( pbef(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) ),   & 
    677 !            &             paft(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) )  ) 
    678          zbup(:,:) = MAX( pbef(:,:,jl) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ),   & 
    679             &             paft(:,:,jl) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) )  ) 
    680          zbdo(:,:) = MIN( pbef(:,:,jl) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ),   & 
    681             &             paft(:,:,jl) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) )  ) 
    682  
    683       z1_dt = 1._wp / pdt 
    684       DO jj = 2, jpjm1 
    685          DO ji = fs_2, fs_jpim1   ! vector opt. 
     1551            DO ji = fs_2, fs_jpim1   ! vector opt. 
     1552               ! 
     1553               IF( .NOT. ll_9points ) THEN 
     1554                  zup  = MAX( zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ), zbup(ji  ,jj-1), zbup(ji  ,jj+1) )  ! search max/min in neighbourhood 
     1555                  zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1) ) 
     1556                  ! 
     1557               ELSE 
     1558                  zup  = MAX( zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ), zbup(ji  ,jj-1), zbup(ji  ,jj+1), &  ! search max/min in neighbourhood 
     1559                     &                     zbup(ji-1,jj-1), zbup(ji+1,jj+1), zbup(ji+1,jj-1), zbup(ji-1,jj+1)  ) 
     1560                  zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1), & 
     1561                     &                     zbdo(ji-1,jj-1), zbdo(ji+1,jj+1), zbdo(ji+1,jj-1), zbdo(ji-1,jj+1)  ) 
     1562               ENDIF 
     1563               ! 
     1564               zpos = MAX( 0., pfu_ho(ji-1,jj,jl) ) - MIN( 0., pfu_ho(ji  ,jj,jl) ) &  ! positive/negative part of the flux 
     1565                  & + MAX( 0., pfv_ho(ji,jj-1,jl) ) - MIN( 0., pfv_ho(ji,jj  ,jl) ) 
     1566               zneg = MAX( 0., pfu_ho(ji  ,jj,jl) ) - MIN( 0., pfu_ho(ji-1,jj,jl) ) & 
     1567                  & + MAX( 0., pfv_ho(ji,jj  ,jl) ) - MIN( 0., pfv_ho(ji,jj-1,jl) ) 
     1568               ! 
     1569               IF( ll_HgradU .AND. .NOT.ll_gurvan ) THEN 
     1570                  zneg2 = (   pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) & 
     1571                     &    ) * ( 1. - pamsk ) 
     1572                  zpos2 = ( - pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) - pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) & 
     1573                     &    ) * ( 1. - pamsk ) 
     1574               ELSE 
     1575                  zneg2 = 0. ; zpos2 = 0. 
     1576               ENDIF 
     1577               ! 
     1578               !                                  ! up & down beta terms 
     1579               IF( (zpos+zpos2) > 0. ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_low(ji,jj,jl) ) / (zpos+zpos2) * e1e2t(ji,jj) * z1_dt 
     1580               ELSE                         ; zbetup(ji,jj,jl) = 0. ! zbig 
     1581               ENDIF 
     1582               ! 
     1583               IF( (zneg+zneg2) > 0. ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_low(ji,jj,jl) - zdo ) / (zneg+zneg2) * e1e2t(ji,jj) * z1_dt 
     1584               ELSE                         ; zbetdo(ji,jj,jl) = 0. ! zbig 
     1585               ENDIF 
     1586               ! 
     1587               ! if all the points are outside ice cover 
     1588               IF( zup == -zbig )   zbetup(ji,jj,jl) = 0. ! zbig 
     1589               IF( zdo ==  zbig )   zbetdo(ji,jj,jl) = 0. ! zbig             
     1590               ! 
     1591 
     1592!!            IF( ji==26 .AND. jj==86) THEN 
     1593!               WRITE(numout,*) '-----------------' 
     1594!               WRITE(numout,*) 'zpos',zpos,zpos2 
     1595!               WRITE(numout,*) 'zneg',zneg,zneg2 
     1596!               WRITE(numout,*) 'puc/pu',ABS(puc(ji,jj))/MAX(epsi20, ABS(pu(ji,jj))) 
     1597!               WRITE(numout,*) 'pvc/pv',ABS(pvc(ji,jj))/MAX(epsi20, ABS(pv(ji,jj))) 
     1598!               WRITE(numout,*) 'pucm1/pu',ABS(puc(ji-1,jj))/MAX(epsi20, ABS(pu(ji-1,jj))) 
     1599!               WRITE(numout,*) 'pvcm1/pv',ABS(pvc(ji,jj-1))/MAX(epsi20, ABS(pv(ji,jj-1))) 
     1600!               WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)+pfu_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt 
     1601!               WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)+pfv_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt 
     1602!               WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)+pfu_low(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt 
     1603!               WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)+pfv_low(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt 
     1604!               WRITE(numout,*) 'pfu_low',pfu_low(ji,jj) * r1_e1e2t(ji,jj) * pdt 
     1605!               WRITE(numout,*) 'pfv_low',pfv_low(ji,jj) * r1_e1e2t(ji,jj) * pdt 
     1606!               WRITE(numout,*) 'pfu_lowm1',pfu_low(ji-1,jj) * r1_e1e2t(ji,jj) * pdt 
     1607!               WRITE(numout,*) 'pfv_lowm1',pfv_low(ji,jj-1) * r1_e1e2t(ji,jj) * pdt 
     1608!                
     1609!               WRITE(numout,*) 'pt',pt(ji,jj) 
     1610!               WRITE(numout,*) 'ptim1',pt(ji-1,jj) 
     1611!               WRITE(numout,*) 'ptjm1',pt(ji,jj-1) 
     1612!               WRITE(numout,*) 'pt_low',pt_low(ji,jj) 
     1613!               WRITE(numout,*) 'zbetup',zbetup(ji,jj) 
     1614!               WRITE(numout,*) 'zbetdo',zbetdo(ji,jj) 
     1615!               WRITE(numout,*) 'zup',zup 
     1616!               WRITE(numout,*) 'zdo',zdo 
     1617!            ENDIF 
    6861618            ! 
    687             zup  = MAX(   zbup(ji,jj), zbup(ji-1,jj  ), zbup(ji+1,jj  ),   &        ! search max/min in neighbourhood 
    688                &                       zbup(ji  ,jj-1), zbup(ji  ,jj+1)    ) 
    689             zdo  = MIN(   zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ),   & 
    690                &                       zbdo(ji  ,jj-1), zbdo(ji  ,jj+1)    ) 
    691                ! 
    692                zpos = MAX( 0., paa(ji-1,jj  ,jl) ) - MIN( 0., paa(ji  ,jj  ,jl) )   &        ! positive/negative  part of the flux 
    693                   & + MAX( 0., pbb(ji  ,jj-1,jl) ) - MIN( 0., pbb(ji  ,jj  ,jl) ) 
    694                zneg = MAX( 0., paa(ji  ,jj  ,jl) ) - MIN( 0., paa(ji-1,jj  ,jl) )   & 
    695                   & + MAX( 0., pbb(ji  ,jj  ,jl) ) - MIN( 0., pbb(ji  ,jj-1,jl) ) 
    696                ! 
    697                zbt = e1e2t(ji,jj) * z1_dt                                   ! up & down beta terms 
    698                zbetup(ji,jj,jl) = ( zup            - paft(ji,jj,jl) ) / ( zpos + zsml ) * zbt 
    699                zbetdo(ji,jj,jl) = ( paft(ji,jj,jl) - zdo         )    / ( zneg + zsml ) * zbt 
    7001619            END DO 
    7011620         END DO 
    … …  
    7031622      CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign) 
    7041623 
    705       DO jl = 1, ipl 
    706          ! monotonic flux in the i & j direction (paa & pbb) 
    707          ! ------------------------------------- 
    708          DO jj = 2, jpjm1 
     1624       
     1625      ! monotonic flux in the y direction 
     1626      ! --------------------------------- 
     1627      DO jl = 1, jpl 
     1628         DO jj = 1, jpjm1 
    7091629            DO ji = 1, fs_jpim1   ! vector opt. 
    7101630               zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 
    7111631               zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) 
    712                zcu = 0.5  + SIGN( 0.5 , paa(ji,jj,jl) ) 
     1632               zcu = 0.5  + SIGN( 0.5 , pfu_ho(ji,jj,jl) ) 
    7131633               ! 
    714                paa(ji,jj,jl) = paa(ji,jj,jl) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 
    715             END DO 
    716          END DO 
    717          ! 
     1634               zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu ) 
     1635 
     1636               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_low(ji,jj,jl) 
     1637 
     1638!!            IF( ji==26 .AND. jj==86) THEN 
     1639!!               WRITE(numout,*) 'coefU',zcoef 
     1640!!               WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt 
     1641!!               WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj,jl)) * r1_e1e2t(ji,jj) * pdt 
     1642!!            ENDIF 
     1643 
     1644            END DO 
     1645         END DO 
     1646 
    7181647         DO jj = 1, jpjm1 
    719             DO ji = fs_2, fs_jpim1   ! vector opt. 
     1648            DO ji = 1, fs_jpim1   ! vector opt. 
    7201649               zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 
    7211650               zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) 
    722                zcv = 0.5  + SIGN( 0.5 , pbb(ji,jj,jl) ) 
     1651               zcv = 0.5  + SIGN( 0.5 , pfv_ho(ji,jj,jl) ) 
    7231652               ! 
    724                pbb(ji,jj,jl) = pbb(ji,jj,jl) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 
    725             END DO 
    726          END DO 
     1653               zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv ) 
     1654 
     1655               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_low(ji,jj,jl) 
     1656 
     1657!!            IF( ji==26 .AND. jj==86) THEN 
     1658!!               WRITE(numout,*) 'coefV',zcoef 
     1659!!               WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt 
     1660!!               WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1,jl)) * r1_e1e2t(ji,jj) * pdt 
     1661!!            ENDIF 
     1662            END DO 
     1663         END DO 
     1664 
     1665         ! clem test 
     1666         DO jj = 2, jpjm1 
     1667            DO ji = 2, fs_jpim1   ! vector opt. 
     1668               IF( ll_gurvan ) THEN 
     1669                  zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1670                     &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1671               ELSE 
     1672                  zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1673                     &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1674                     &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1675                     &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1676                     &         ) * tmask(ji,jj,1) 
     1677               ENDIF 
     1678               IF( zzt(ji,jj,jl) < -epsi20 ) THEN 
     1679                  WRITE(numout,*) 'T<0 nonosc',zzt(ji,jj,jl) 
     1680               ENDIF 
     1681            END DO 
     1682         END DO 
     1683 
    7271684      END DO 
     1685 
     1686      ! 
    7281687      ! 
    7291688   END SUBROUTINE nonosc_2d 
     1689 
     1690   SUBROUTINE limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
     1691      !!--------------------------------------------------------------------- 
     1692      !!                    ***  ROUTINE limiter_x  *** 
     1693      !!      
     1694      !! **  Purpose :   compute flux limiter  
     1695      !!---------------------------------------------------------------------- 
     1696      REAL(wp)                   , INTENT(in   )           ::   pdt          ! tracer time-step 
     1697      REAL(wp), DIMENSION (:,:  ), INTENT(in   )           ::   pu           ! ice i-velocity => u*e2 
     1698      REAL(wp), DIMENSION (:,:,:), INTENT(in   )           ::   puc          ! ice i-velocity *A => u*e2*a 
     1699      REAL(wp), DIMENSION (:,:,:), INTENT(in   )           ::   pt           ! ice tracer 
     1700      REAL(wp), DIMENSION (:,:,:), INTENT(inout)           ::   pfu_ho       ! high order flux 
     1701      REAL(wp), DIMENSION (:,:,:), INTENT(in   ), OPTIONAL ::   pfu_ups      ! upstream flux 
     1702      ! 
     1703      REAL(wp) ::   Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr 
     1704      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
     1705      REAL(wp), DIMENSION (jpi,jpj,jpl) ::   zslpx       ! tracer slopes  
     1706      !!---------------------------------------------------------------------- 
     1707      ! 
     1708      DO jl = 1, jpl 
     1709         DO jj = 2, jpjm1 
     1710            DO ji = fs_2, fs_jpim1   ! vector opt. 
     1711               zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1) 
     1712            END DO 
     1713         END DO 
     1714      END DO 
     1715      CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.)   ! lateral boundary cond. 
     1716       
     1717      DO jl = 1, jpl 
     1718         DO jj = 2, jpjm1 
     1719            DO ji = fs_2, fs_jpim1   ! vector opt. 
     1720               uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj) 
     1721                
     1722               Rjm = zslpx(ji-1,jj,jl) 
     1723               Rj  = zslpx(ji  ,jj,jl) 
     1724               Rjp = zslpx(ji+1,jj,jl) 
     1725 
     1726               IF( PRESENT(pfu_ups) ) THEN 
     1727 
     1728                  IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm 
     1729                  ELSE                        ;   Rr = Rjp 
     1730                  ENDIF 
     1731 
     1732                  zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)      
     1733                  IF( Rj > 0. ) THEN 
     1734                     zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pu(ji,jj)),  & 
     1735                        &        MIN( 2. * Rr * 0.5 * ABS(pu(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) ) 
     1736                  ELSE 
     1737                     zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)),  & 
     1738                        &        MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) ) 
     1739                  ENDIF 
     1740                  pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter 
     1741 
     1742               ELSE 
     1743                  IF( Rj /= 0. ) THEN 
     1744                     IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj 
     1745                     ELSE                        ;   Cr = Rjp / Rj 
     1746                     ENDIF 
     1747                  ELSE 
     1748                     Cr = 0. 
     1749                     !IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm * 1.e20 
     1750                     !ELSE                        ;   Cr = Rjp * 1.e20 
     1751                     !ENDIF 
     1752                  ENDIF 
     1753 
     1754                  ! -- superbee -- 
     1755                  zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 
     1756                  ! -- van albada 2 -- 
     1757                  !!zpsi = 2.*Cr / (Cr*Cr+1.) 
     1758 
     1759                  ! -- sweby (with beta=1) -- 
     1760                  !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 
     1761                  ! -- van Leer -- 
     1762                  !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 
     1763                  ! -- ospre -- 
     1764                  !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 
     1765                  ! -- koren -- 
     1766                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 
     1767                  ! -- charm -- 
     1768                  !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 
     1769                  !ELSE                 ;   zpsi = 0. 
     1770                  !ENDIF 
     1771                  ! -- van albada 1 -- 
     1772                  !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 
     1773                  ! -- smart -- 
     1774                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 
     1775                  ! -- umist -- 
     1776                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 
     1777 
     1778                  ! high order flux corrected by the limiter 
     1779                  pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5 
     1780 
     1781               ENDIF 
     1782            END DO 
     1783         END DO 
     1784      END DO 
     1785      CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.)   ! lateral boundary cond. 
     1786      ! 
     1787   END SUBROUTINE limiter_x 
     1788 
     1789   SUBROUTINE limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     1790      !!--------------------------------------------------------------------- 
     1791      !!                    ***  ROUTINE limiter_y  *** 
     1792      !!      
     1793      !! **  Purpose :   compute flux limiter  
     1794      !!---------------------------------------------------------------------- 
     1795      REAL(wp)                   , INTENT(in   )           ::   pdt          ! tracer time-step 
     1796      REAL(wp), DIMENSION (:,:  ), INTENT(in   )           ::   pv           ! ice i-velocity => u*e2 
     1797      REAL(wp), DIMENSION (:,:,:), INTENT(in   )           ::   pvc          ! ice i-velocity *A => u*e2*a 
     1798      REAL(wp), DIMENSION (:,:,:), INTENT(in   )           ::   pt           ! ice tracer 
     1799      REAL(wp), DIMENSION (:,:,:), INTENT(inout)           ::   pfv_ho       ! high order flux 
     1800      REAL(wp), DIMENSION (:,:,:), INTENT(in   ), OPTIONAL ::   pfv_ups      ! upstream flux 
     1801      ! 
     1802      REAL(wp) ::   Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr 
     1803      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
     1804      REAL(wp), DIMENSION (jpi,jpj,jpl) ::   zslpy       ! tracer slopes  
     1805      !!---------------------------------------------------------------------- 
     1806      ! 
     1807      DO jl = 1, jpl 
     1808         DO jj = 2, jpjm1 
     1809            DO ji = fs_2, fs_jpim1   ! vector opt. 
     1810               zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1) 
     1811            END DO 
     1812         END DO 
     1813      END DO 
     1814      CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.)   ! lateral boundary cond. 
     1815 
     1816      DO jl = 1, jpl 
     1817         DO jj = 2, jpjm1 
     1818            DO ji = fs_2, fs_jpim1   ! vector opt. 
     1819               vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj) 
     1820 
     1821               Rjm = zslpy(ji,jj-1,jl) 
     1822               Rj  = zslpy(ji,jj  ,jl) 
     1823               Rjp = zslpy(ji,jj+1,jl) 
     1824 
     1825               IF( PRESENT(pfv_ups) ) THEN 
     1826 
     1827                  IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm 
     1828                  ELSE                        ;   Rr = Rjp 
     1829                  ENDIF 
     1830 
     1831                  zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)      
     1832                  IF( Rj > 0. ) THEN 
     1833                     zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pv(ji,jj)),  & 
     1834                        &        MIN( 2. * Rr * 0.5 * ABS(pv(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) ) 
     1835                  ELSE 
     1836                     zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)),  & 
     1837                        &        MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) ) 
     1838                  ENDIF 
     1839                  pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter 
     1840 
     1841               ELSE 
     1842 
     1843                  IF( Rj /= 0. ) THEN 
     1844                     IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj 
     1845                     ELSE                        ;   Cr = Rjp / Rj 
     1846                     ENDIF 
     1847                  ELSE 
     1848                     Cr = 0. 
     1849                     !IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm * 1.e20 
     1850                     !ELSE                        ;   Cr = Rjp * 1.e20 
     1851                     !ENDIF 
     1852                  ENDIF 
     1853 
     1854                  ! -- superbee -- 
     1855                  zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 
     1856                  ! -- van albada 2 -- 
     1857                  !!zpsi = 2.*Cr / (Cr*Cr+1.) 
     1858 
     1859                  ! -- sweby (with beta=1) -- 
     1860                  !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 
     1861                  ! -- van Leer -- 
     1862                  !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 
     1863                  ! -- ospre -- 
     1864                  !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 
     1865                  ! -- koren -- 
     1866                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 
     1867                  ! -- charm -- 
     1868                  !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 
     1869                  !ELSE                 ;   zpsi = 0. 
     1870                  !ENDIF 
     1871                  ! -- van albada 1 -- 
     1872                  !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 
     1873                  ! -- smart -- 
     1874                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 
     1875                  ! -- umist -- 
     1876                  !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 
     1877 
     1878                  ! high order flux corrected by the limiter 
     1879                  pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5 
     1880 
     1881               ENDIF 
     1882            END DO 
     1883         END DO 
     1884      END DO 
     1885      CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.)   ! lateral boundary cond. 
     1886      ! 
     1887   END SUBROUTINE limiter_y 
    7301888 
    7311889#else 
Note: See TracChangeset for help on using the changeset viewer.