Changeset 10439


Ignore:
Timestamp:
2018-12-21T12:18:00+01:00 (21 months ago)
Author:
clem
Message:

debug zalesak prelimiter in ice advection

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/ICE/icedyn_adv_umx.F90

    r10425 r10439  
    3434   REAL(wp) ::   z1_6   = 1._wp /   6._wp   ! =1/6 
    3535   REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120 
    36  
    37    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: amaxu, amaxv 
    3836    
    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  
    4837   ! 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.  ! 
     38   LOGICAL :: ll_neg = .TRUE. 
    6839    
    6940   ! alternate directions for upstream 
     
    7445    
    7546   ! 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     
     47   LOGICAL :: ll_prelimiter_zalesak = .TRUE.  ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D 
     48 
    8249 
    8350   !! * Substitutions 
     
    12087      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers 
    12188      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 
     89      REAL(wp), DIMENSION(1)           ::   zcflprv, zcflnow   ! send zcflnow and receive zcflprv 
     90      REAL(wp), DIMENSION(jpi,jpj)     ::   zudy, zvdx, zcu_box, zcv_box  
     91      REAL(wp), DIMENSION(jpi,jpj)     ::   zati1, zati2 
    12992      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho 
    13093      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai, z1_aip 
    13194      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhvar 
    132  
    13395      !!---------------------------------------------------------------------- 
    13496      ! 
    13597      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' 
    136       ! 
    13798      ! 
    13899      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- ! 
     
    170131      END DO 
    171132 
    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 
    176133      !---------------! 
    177134      !== advection ==! 
     
    179136      DO jt = 1, icycle 
    180137 
    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 
     138         ! record at_i before advection (for open water) 
     139         zati1(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
    188140          
     141         ! inverse of A and Ap 
    189142         WHERE( pa_i(:,:,:) >= epsi20 )   ;   z1_ai(:,:,:) = 1._wp / pa_i(:,:,:) 
    190143         ELSEWHERE                        ;   z1_ai(:,:,:) = 0. 
    191144         END WHERE 
    192             ! 
    193145         WHERE( pa_ip(:,:,:) >= epsi20 )  ;   z1_aip(:,:,:) = 1._wp / pa_ip(:,:,:) 
    194146         ELSEWHERE                        ;   z1_aip(:,:,:) = 0. 
    195147         END WHERE 
    196148         ! 
    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          ! 
     149         ! set u*a=u for advection of A only  
    211150         DO jl = 1, jpl 
    212151            zua_ho(:,:,jl) = zudy(:,:) 
     
    218157         zamsk = 0._wp 
    219158         ! 
    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             ! 
    227159         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(:,:,:) 
     160         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_i )                ! Ice volume 
    230161         ! 
    231162         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(:,:,:) 
     163         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_s )                ! Snw volume 
    234164         ! 
    235165         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 
     166         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, psv_i )               ! Salt content 
    237167         ! 
    238168         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 
     169         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, poa_i )               ! Age content 
    240170         ! 
    241171         DO jk = 1, nlay_i 
    242172            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 
     173            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 
    244174         END DO 
    245175         ! 
    246176         DO jk = 1, nlay_s 
    247177            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             ! 
     178            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 
     179         END DO 
     180         ! 
    251181         IF ( ln_pnd_H12 ) THEN 
    252             ! 
     182            ! set u*a=u for advection of Ap only  
    253183            DO jl = 1, jpl 
    254184               zua_ho(:,:,jl) = zudy(:,:) 
     
    261191            ! 
    262192            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 
     193            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_ip )                ! mp volume 
    264194         ENDIF 
    265195         ! 
    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 
     196         zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
     197         DO jj = 2, jpjm1 
     198            DO ji = fs_2, fs_jpim1 
     199               pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                                                   ! Open water area 
     200                  &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 
     201            END DO 
     202         END DO 
     203         CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T',  1. ) 
    277204         ! 
    278205      END DO 
     
    294221      !! ** Action : - pt  the after advective tracer 
    295222      !!---------------------------------------------------------------------- 
    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 
     223      REAL(wp)                        , INTENT(in   )           ::   pamsk          ! advection of concentration (1) or other tracers (0) 
     224      INTEGER                         , INTENT(in   )           ::   kn_umx         ! order of the scheme (1-5=UM or 20=CEN2) 
     225      INTEGER                         , INTENT(in   )           ::   jt             ! number of sub-iteration 
     226      INTEGER                         , INTENT(in   )           ::   kt             ! number of iteration 
     227      REAL(wp)                        , INTENT(in   )           ::   pdt            ! tracer time-step 
     228      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pu   , pv      ! 2 ice velocity components => u*e2 
     229      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   )           ::   puc  , pvc     ! 2 ice velocity components => u*e2 or u*a*e2u 
     230      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pubox, pvbox   ! upstream velocity 
     231      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   pt             ! tracer field 
     232      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   ptc            ! tracer content field 
    306233      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho ! high order u*a fluxes 
    307234      ! 
     
    310237      INTEGER  ::   kn_limiter = 1   ! 1=nonosc ; 2=superbee ; 3=h3(rachid) 
    311238      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  
     239      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zfu_ups, zfv_ups, zt_ups 
    313240      !!---------------------------------------------------------------------- 
    314241      ! 
    315242      !  upstream (_ups) advection with initial mass fluxes 
    316243      ! --------------------------------------------------- 
    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 
     244      ! 
     245      IF( .NOT. ll_upsxy ) THEN         !** no alternate directions **! 
    335246         DO jl = 1, jpl 
    336247            DO jj = 1, jpjm1 
    337248               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 
     249                  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) 
     250                  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) 
     251               END DO 
     252            END DO 
     253         END DO 
     254 
     255      ELSE                              !** alternate directions **! 
    350256         ! 
    351257         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 
     258            ! 
     259            DO jl = 1, jpl              !-- flux in x-direction 
    354260               DO jj = 1, jpjm1 
    355261                  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 
     262                     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) 
    361263                  END DO 
    362264               END DO 
    363265            END DO 
    364              
    365             ! first guess of tracer content from u-flux 
    366             DO jl = 1, jpl 
     266            ! 
     267            DO jl = 1, jpl              !-- first guess of tracer from u-flux 
    367268               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 
     269                  DO ji = fs_2, fs_jpim1 
     270                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     271                        &            + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     272                        &            ) * tmask(ji,jj,1) 
    385273                  END DO 
    386274               END DO 
     
    388276            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
    389277            ! 
    390             ! flux in y-direction 
    391             DO jl = 1, jpl 
     278            DO jl = 1, jpl              !-- flux in y-direction 
    392279               DO jj = 1, jpjm1 
    393280                  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 
     281                     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) 
    399282                  END DO 
    400283               END DO 
    401284            END DO 
    402          ! 
     285            ! 
    403286         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
    404             ! flux in y-direction 
    405             DO jl = 1, jpl 
     287            ! 
     288            DO jl = 1, jpl              !-- flux in y-direction 
    406289               DO jj = 1, jpjm1 
    407290                  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 
     291                     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) 
    413292                  END DO 
    414293               END DO 
    415294            END DO 
    416  
    417             ! first guess of tracer content from v-flux 
    418             DO jl = 1, jpl 
     295            ! 
     296            DO jl = 1, jpl              !-- first guess of tracer from v-flux 
    419297               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 
     298                  DO ji = fs_2, fs_jpim1 
     299                     zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 
     300                        &            + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 
     301                        &            * tmask(ji,jj,1) 
    437302                  END DO 
    438303               END DO 
     
    440305            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
    441306            ! 
    442             ! flux in x-direction 
    443             DO jl = 1, jpl 
     307            DO jl = 1, jpl              !-- flux in x-direction 
    444308               DO jj = 1, jpjm1 
    445309                  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 
     310                     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) 
    451311                  END DO 
    452312               END DO 
     
    456316          
    457317      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 
     318      ! 
     319      DO jl = 1, jpl                    !-- after tracer with upstream scheme 
    475320         DO jj = 2, jpjm1 
    476321            DO ji = fs_2, fs_jpim1 
    477322               ztra          = - (   zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj  ,jl)   & 
    478323                  &                + 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 
     324               zt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + pdt * ztra + ( pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) )   & 
     325                  &                                            +   pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) ) & 
     326                  &                                              * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 
    494327            END DO 
    495328         END DO 
     
    512345         ! 
    513346      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)  
     347      ! 
     348      !              --ho    --ho 
     349      ! new fluxes = u*H  *  u*a / u 
     350      ! ---------------------------- 
     351      IF( pamsk == 0. ) THEN 
     352         DO jl = 1, jpl 
     353            DO jj = 1, jpjm1 
     354               DO ji = 1, fs_jpim1 
     355                  IF( ABS( puc(ji,jj,jl) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN 
     356                     zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj) 
     357                     zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj) 
    523358                  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) 
     359                     zfu_ho (ji,jj,jl) = 0._wp 
     360                     zfu_ups(ji,jj,jl) = 0._wp 
    527361                  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) 
     362                  ! 
     363                  IF( ABS( pvc(ji,jj,jl) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN 
     364                     zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj) 
     365                     zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj) 
     366                  ELSE 
     367                     zfv_ho (ji,jj,jl) = 0._wp   
     368                     zfv_ups(ji,jj,jl) = 0._wp   
    532369                  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. ) 
     370               END DO 
     371            END DO 
     372         END DO 
    543373      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       ! ---------------------------- 
     374      ! 
     375      ! in case of advection of A: output u*a(ho) 
     376      ! ----------------------------------------- 
    593377      IF( PRESENT( pua_ho ) ) THEN 
    594378         DO jl = 1, jpl 
     
    601385         END DO 
    602386      ENDIF 
    603  
    604  
    605       IF( .NOT.ll_thickness ) THEN 
    606          ! final trend with corrected fluxes 
    607          ! ------------------------------------ 
    608          DO jl = 1, jpl 
    609             DO jj = 2, jpjm1 
    610                DO ji = fs_2, fs_jpim1  
    611                   ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj) * pdt   
    612                    
    613                   ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra ) * tmask(ji,jj,1) 
    614                    
    615                   !!               IF( ji==26 .AND. jj==86) THEN 
    616                   !!                  WRITE(numout,*) 'ztc high order',ptc(ji,jj) 
    617                   !!               ENDIF 
    618                    
    619                END DO 
    620             END DO 
    621          END DO 
    622          CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T',  1. ) 
    623       ENDIF 
    624        
     387      ! 
     388      ! final trend with corrected fluxes 
     389      ! --------------------------------- 
     390      DO jl = 1, jpl 
     391         DO jj = 2, jpjm1 
     392            DO ji = fs_2, fs_jpim1  
     393               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   
     394                
     395               ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra ) * tmask(ji,jj,1)                
     396            END DO 
     397         END DO 
     398      END DO 
     399      CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T',  1. ) 
    625400      ! 
    626401   END SUBROUTINE adv_umx 
     402 
    627403 
    628404   SUBROUTINE cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, & 
     
    638414      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    639415      !!---------------------------------------------------------------------- 
    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  
     416      REAL(wp)                        , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0) 
     417      INTEGER                         , INTENT(in   ) ::   kn_limiter       ! limiter 
     418      INTEGER                         , INTENT(in   ) ::   jt               ! number of sub-iteration 
     419      INTEGER                         , INTENT(in   ) ::   kt               ! number of iteration 
     420      REAL(wp)                        , INTENT(in   ) ::   pdt              ! tracer time-step 
     421      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt               ! tracer fields 
     422      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pu, pv           ! 2 ice velocity components 
     423      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components 
     424      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   ptc              ! tracer content at before time step  
    649425      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  
     426      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content  
     427      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes  
    652428      ! 
    653429      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    654       LOGICAL  ::   ll_xy = .TRUE. 
    655430      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zzt 
    656431      !!---------------------------------------------------------------------- 
    657432      ! 
    658       IF( .NOT.ll_xy ) THEN   !-- no alternate directions --! 
     433      IF( .NOT.ll_hoxy ) THEN           !** no alternate directions **! 
    659434         ! 
    660435         DO jl = 1, jpl 
    661436            DO jj = 1, jpjm1 
    662437               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 
     438                  pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 
     439                  pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 
    670440               END DO 
    671441            END DO 
    672442         END DO 
    673443         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 
     444            CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
    679445         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 ) 
     446            CALL limiter_x( pdt, pu, pt, pfu_ho ) 
     447            CALL limiter_y( pdt, pv, pt, pfv_ho ) 
    682448         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 ) 
     449            CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 
     450            CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 
    685451         ENDIF 
    686452         ! 
    687       ELSE                    !-- alternate directions --! 
     453      ELSE                              !** alternate directions **! 
    688454         ! 
    689455         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
    690456            ! 
    691             ! flux in x-direction 
    692             DO jl = 1, jpl 
     457            DO jl = 1, jpl              !-- flux in x-direction 
    693458               DO jj = 1, jpjm1 
    694459                  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 
     460                     pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 
    700461                  END DO 
    701462               END DO 
    702463            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 
     464            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, pt, pfu_ho ) 
     465            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 
     466 
     467            DO jl = 1, jpl              !-- first guess of tracer from u-flux 
    708468               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 
     469                  DO ji = fs_2, fs_jpim1 
     470                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)        & 
     471                        &            + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 
     472                        &         * tmask(ji,jj,1) 
    717473                  END DO 
    718474               END DO 
     
    720476            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
    721477 
    722             ! flux in y-direction 
    723             DO jl = 1, jpl 
     478            DO jl = 1, jpl              !-- flux in y-direction 
    724479               DO jj = 1, jpjm1 
    725480                  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 
     481                     pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji,jj+1,jl) ) 
    731482                  END DO 
    732483               END DO 
    733484            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 ) 
     485            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pt, pfv_ho ) 
     486            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 
    736487 
    737488         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
    738489            ! 
    739             ! flux in y-direction 
    740             DO jl = 1, jpl 
     490            DO jl = 1, jpl              !-- flux in y-direction 
    741491               DO jj = 1, jpjm1 
    742492                  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 
     493                     pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 
    748494                  END DO 
    749495               END DO 
    750496            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 ) 
     497            IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pt, pfv_ho ) 
     498            IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 
    753499            ! 
    754             ! first guess of tracer content from v-flux 
    755             DO jl = 1, jpl 
     500            DO jl = 1, jpl              !-- first guess of tracer from v-flux 
    756501               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 
     502                  DO ji = fs_2, fs_jpim1 
     503                     zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 
     504                        &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 
     505                        &         * tmask(ji,jj,1) 
    765506                  END DO 
    766507               END DO 
     
    768509            CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
    769510            ! 
    770             ! flux in x-direction 
    771             DO jl = 1, jpl 
     511            DO jl = 1, jpl              !-- flux in x-direction 
    772512               DO jj = 1, jpjm1 
    773513                  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 
     514                     pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji+1,jj,jl) ) 
    779515                  END DO 
    780516               END DO 
    781517            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 ) 
     518            IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, pt, pfu_ho ) 
     519            IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 
    784520 
    785521         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 
     522         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 ) 
    791523          
    792524      ENDIF 
     
    807539      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    808540      !!---------------------------------------------------------------------- 
    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  
     541      REAL(wp)                        , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0) 
     542      INTEGER                         , INTENT(in   ) ::   kn_limiter       ! limiter 
     543      INTEGER                         , INTENT(in   ) ::   kn_umx           ! order of the scheme (1-5=UM or 20=CEN2) 
     544      INTEGER                         , INTENT(in   ) ::   jt               ! number of sub-iteration 
     545      INTEGER                         , INTENT(in   ) ::   kt               ! number of iteration 
     546      REAL(wp)                        , INTENT(in   ) ::   pdt              ! tracer time-step 
     547      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt               ! tracer fields 
     548      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pu, pv           ! 2 ice velocity components 
     549      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   puc, pvc         ! 2 ice velocity * A components 
     550      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pubox, pvbox     ! upstream velocity 
     551      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   ptc              ! tracer content at before time step  
    820552      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_u, pt_v       ! tracer at u- and v-points  
    821553      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  
     554      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt_ups           ! upstream guess of tracer content  
     555      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout) ::   pfu_ups, pfv_ups ! upstream fluxes  
    824556      ! 
    825557      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
     
    831563         ! 
    832564         !                                                        !--  ultimate interpolation of pt at u-point  --! 
    833          CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) 
     565         CALL ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho ) 
    834566         !                                                        !--  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 ) 
     567         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, pt, pfu_ho ) 
     568         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 
    837569         !                                                        !--  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 
     570         DO jl = 1, jpl 
     571            DO jj = 2, jpjm1 
     572               DO ji = fs_2, fs_jpim1 
     573                  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)  & 
     574                     &                         - pt   (ji,jj,jl) * pdt * ( pu  (ji,jj) - pu  (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk 
     575                  zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 
     576               END DO 
     577            END DO 
     578         END DO 
     579         CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
    879580         ! 
    880581         !                                                        !--  ultimate interpolation of pt at v-point  --! 
    881582         IF( ll_hoxy ) THEN 
    882             CALL ultimate_y( kn_umx, pdt, zzt, pv, pvc, pt_v, pfv_ho ) 
     583            CALL ultimate_y( kn_umx, pdt, zzt, pv, pt_v, pfv_ho ) 
    883584         ELSE 
    884             CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) 
     585            CALL ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho ) 
    885586         ENDIF 
    886587         !                                                        !--  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 ) 
     588         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pt, pfv_ho ) 
     589         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 
    889590         !          
    890591         ! 
     
    892593         ! 
    893594         !                                                        !--  ultimate interpolation of pt at v-point  --! 
    894          CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) 
     595         CALL ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho ) 
    895596         !                                                        !--  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 ) 
     597         IF( kn_limiter == 2 )   CALL limiter_y( pdt, pv, pt, pfv_ho ) 
     598         IF( kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 
    898599         !                                                        !--  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 
     600         DO jl = 1, jpl 
     601            DO jj = 2, jpjm1 
     602               DO ji = fs_2, fs_jpim1 
     603                  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)  & 
     604                     &                         - pt   (ji,jj,jl) * pdt * ( pv  (ji,jj) - pv  (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk 
     605                  zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 
     606               END DO 
     607            END DO 
     608         END DO 
     609         CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 
    940610         ! 
    941611         !                                                        !--  ultimate interpolation of pt at u-point  --! 
    942612         IF( ll_hoxy ) THEN 
    943             CALL ultimate_x( kn_umx, pdt, zzt, pu, puc, pt_u, pfu_ho ) 
     613            CALL ultimate_x( kn_umx, pdt, zzt, pu, pt_u, pfu_ho ) 
    944614         ELSE 
    945             CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) 
     615            CALL ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho ) 
    946616         ENDIF 
    947617         !                                                        !--  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 ) 
     618         IF( kn_limiter == 2 )   CALL limiter_x( pdt, pu, pt, pfu_ho ) 
     619         IF( kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 
    950620         ! 
    951621         ! 
    952622      ENDIF 
    953623 
    954       
    955624      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 
     625         CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
    988626      ENDIF 
    989627      ! 
     
    991629 
    992630 
    993    SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) 
     631   SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho ) 
    994632      !!--------------------------------------------------------------------- 
    995633      !!                    ***  ROUTINE ultimate_x  *** 
     
    1002640      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    1003641      !!---------------------------------------------------------------------- 
    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 
     642      INTEGER                         , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2) 
     643      REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step 
     644      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pu        ! ice i-velocity component 
     645      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt        ! tracer fields 
    1009646      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_u      ! tracer at u-point  
    1010647      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfu_ho    ! high order flux  
    1011648      ! 
    1012649      INTEGER  ::   ji, jj, jl             ! dummy loop indices 
    1013       REAL(wp) ::   zcu, zdx2, zdx4    !   -      - 
     650      REAL(wp) ::   zcu, zdx2, zdx4        !   -      - 
    1014651      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztu1, ztu2, ztu3, ztu4 
    1015652      !!---------------------------------------------------------------------- 
     
    1076713                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    1077714                  zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    1078 !!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
     715!!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    1079716                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        & 
    1080717                     &                                               -              zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   & 
     
    1092729                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    1093730                  zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    1094 !!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
     731!!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    1095732                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        & 
    1096733                     &                                               -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   & 
     
    1108745                  zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    1109746                  zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    1110                   !!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj) 
     747!!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    1111748                  zdx4 = zdx2 * zdx2 
    1112749                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (            (                   pt  (ji+1,jj,jl) + pt  (ji,jj,jl)       & 
     
    1121758         ! 
    1122759      END SELECT 
    1123       !                                                     !-- High order flux in i-direction  --! 
     760      ! 
     761      ! if pt at u-point is negative then use the upstream value 
     762      !    this should not be necessary if a proper sea-ice mask is set in Ultimate 
     763      !    to degrade the order of the scheme when necessary (for ex. at the ice edge) 
    1124764      IF( ll_neg ) THEN 
    1125765         DO jl = 1, jpl 
     
    1134774         END DO 
    1135775      ENDIF 
    1136  
     776      !                                                     !-- High order flux in i-direction  --! 
    1137777      DO jl = 1, jpl 
    1138778         DO jj = 1, jpjm1 
    1139779            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 
     780               pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 
    1145781            END DO 
    1146782         END DO 
     
    1150786    
    1151787  
    1152    SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) 
     788   SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho ) 
    1153789      !!--------------------------------------------------------------------- 
    1154790      !!                    ***  ROUTINE ultimate_y  *** 
     
    1161797      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    1162798      !!---------------------------------------------------------------------- 
    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 
     799      INTEGER                         , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2) 
     800      REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step 
     801      REAL(wp), DIMENSION(:,:  )      , INTENT(in   ) ::   pv        ! ice j-velocity component 
     802      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   ) ::   pt        ! tracer fields 
    1168803      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pt_v      ! tracer at v-point  
    1169804      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   pfv_ho    ! high order flux  
    1170805      ! 
    1171       INTEGER  ::   ji, jj, jl       ! dummy loop indices 
     806      INTEGER  ::   ji, jj, jl         ! dummy loop indices 
    1172807      REAL(wp) ::   zcv, zdy2, zdy4    !   -      - 
    1173808      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztv1, ztv2, ztv3, ztv4 
     
    1235870                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    1236871                  zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    1237 !!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
     872!!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    1238873                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)       & 
    1239874                     &                                     -                        zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) )   & 
     
    1250885                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    1251886                  zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    1252 !!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
     887!!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    1253888                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                        ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)       & 
    1254889                     &                                               -          zcv * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) )   & 
     
    1265900                  zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    1266901                  zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    1267 !!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj) 
     902!!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    1268903                  zdy4 = zdy2 * zdy2 
    1269904                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)      & 
     
    1278913         ! 
    1279914      END SELECT 
    1280       !                                                     !-- High order flux in j-direction  --! 
     915      ! 
     916      ! if pt at v-point is negative then use the upstream value 
     917      !    this should not be necessary if a proper sea-ice mask is set in Ultimate 
     918      !    to degrade the order of the scheme when necessary (for ex. at the ice edge) 
    1281919      IF( ll_neg ) THEN 
    1282920         DO jl = 1, jpl 
     
    1291929         END DO 
    1292930      ENDIF 
    1293  
     931      !                                                     !-- High order flux in j-direction  --! 
    1294932      DO jl = 1, jpl 
    1295933         DO jj = 1, jpjm1 
    1296934            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 
     935               pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 
    1302936            END DO 
    1303937         END DO 
     
    1307941      
    1308942 
    1309    SUBROUTINE nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho ) 
     943   SUBROUTINE nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
    1310944      !!--------------------------------------------------------------------- 
    1311945      !!                    ***  ROUTINE nonosc  *** 
     
    1315949      !! 
    1316950      !! **  Method  :   ... ??? 
    1317       !!       warning : pt and pt_low must be masked, but the boundaries 
     951      !!       warning : pt and pt_ups must be masked, but the boundaries 
    1318952      !!       conditions on the fluxes are not necessary zalezak (1979) 
    1319953      !!       drange (1995) multi-dimensional forward-in-time and upstream- 
    1320954      !!       in-space based differencing for fluid 
    1321955      !!---------------------------------------------------------------------- 
    1322       REAL(wp)                     , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0) 
    1323       REAL(wp)                     , INTENT(in   ) ::   pdt              ! tracer time-step 
     956      REAL(wp)                   , INTENT(in   ) ::   pamsk            ! advection of concentration (1) or other tracers (0) 
     957      REAL(wp)                   , INTENT(in   ) ::   pdt              ! tracer time-step 
    1324958      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pu               ! ice i-velocity => u*e2 
    1325959      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   puc              ! ice i-velocity *A => u*e2*a 
    1326960      REAL(wp), DIMENSION (:,:  ), INTENT(in   ) ::   pv               ! ice j-velocity => v*e1 
    1327961      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 
     962      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   ptc, pt, pt_ups  ! before field & upstream guess of after field 
     963      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   pfv_ups, pfu_ups ! upstream flux 
    1330964      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux 
    1331965      ! 
    1332966      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    1333       REAL(wp) ::   zpos, zneg, zbig, zsml, z1_dt, zpos2, zneg2   ! local scalars 
    1334       REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign, zcoef        !   -      - 
     967      REAL(wp) ::   zpos, zneg, zbig, zsml, zup, zdo, z1_dt              ! local scalars 
     968      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zsign, zcoef, zzt      !   -      - 
    1335969      REAL(wp), DIMENSION(jpi,jpj    ) :: zbup, zbdo 
    1336       REAL(wp), DIMENSION(jpi,jpj,jpl) :: zbetup, zbetdo, zti_low, ztj_low, zzt 
     970      REAL(wp), DIMENSION(jpi,jpj,jpl) :: zbetup, zbetdo, zti_ups, ztj_ups 
    1337971      !!---------------------------------------------------------------------- 
    1338972      zbig = 1.e+40_wp 
    1339973      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  
    1424974       
    1425975      ! antidiffusive flux : high order minus low order 
     
    1428978         DO jj = 1, jpjm1 
    1429979            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) 
     980               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 
     981               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 
    1432982            END DO 
    1433983         END DO 
     
    1441991      !                        |      |      |    *   |     
    1442992      !                        |      |      |        |    * 
    1443       !            t_low :       i-1     i       i+1       i+2    
     993      !            t_ups :       i-1     i       i+1       i+2    
    1444994      IF( ll_prelimiter_zalesak ) THEN 
    1445995          
     
    1447997            DO jj = 2, jpjm1 
    1448998               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. ) 
     999                  zti_ups(ji,jj,jl)= pt_ups(ji+1,jj  ,jl) 
     1000                  ztj_ups(ji,jj,jl)= pt_ups(ji  ,jj+1,jl) 
     1001               END DO 
     1002            END DO 
     1003         END DO 
     1004         CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_ups, 'T', 1., ztj_ups, 'T', 1. ) 
    14551005 
    14561006         !! this does not work ?? 
    14571007         !!            DO jj = 2, jpjm1 
    14581008         !!               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) )           & 
     1009         !!                  IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_ups (ji+1,jj  ) - pt_ups (ji  ,jj) ) .AND.     & 
     1010         !!               &      SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_ups (ji  ,jj+1) - pt_ups (ji  ,jj) )           & 
    14611011         !!               &    ) 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) )         & 
     1012         !!                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_ups(ji+1,jj ) - zti_ups(ji  ,jj) ) .AND.   & 
     1013         !!               &         SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_ups(ji,jj+1 ) - ztj_ups(ji  ,jj) )         & 
    14641014         !!               &       ) THEN 
    14651015         !!                        pfu_ho(ji,jj) = 0.  ;   pfv_ho(ji,jj) = 0. 
    14661016         !!                     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) )        & 
     1017         !!                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_ups (ji  ,jj) - pt_ups (ji-1,jj  ) ) .AND.  & 
     1018         !!               &         SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_ups (ji  ,jj) - pt_ups (ji  ,jj-1) )        & 
    14691019         !!               &       )  THEN 
    14701020         !!                        pfu_ho(ji,jj) = 0.  ;   pfv_ho(ji,jj) = 0. 
     
    14771027            DO jj = 2, jpjm1 
    14781028               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 
     1029                  IF ( pfu_ho(ji,jj,jl) * ( pt_ups(ji+1,jj,jl) - pt_ups(ji,jj,jl) ) <= 0. .AND.  & 
     1030                     & pfv_ho(ji,jj,jl) * ( pt_ups(ji,jj+1,jl) - pt_ups(ji,jj,jl) ) <= 0. ) THEN 
    14811031                     ! 
    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. 
     1032                     IF(  pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj,jl) - zti_ups(ji,jj,jl) ) <= 0 .AND.  & 
     1033                        & pfv_ho(ji,jj,jl) * ( ztj_ups(ji,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0) THEN 
     1034                        pfu_ho(ji,jj,jl)=0. 
     1035                        pfv_ho(ji,jj,jl)=0. 
     1036                     ENDIF 
    14841037                     ! 
    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. 
     1038                     IF(  pfu_ho(ji,jj,jl) * ( pt_ups(ji  ,jj,jl) - pt_ups(ji-1,jj,jl) ) <= 0 .AND.  & 
     1039                        & pfv_ho(ji,jj,jl) * ( pt_ups(ji  ,jj,jl) - pt_ups(ji,jj-1,jl) ) <= 0) THEN 
     1040                        pfu_ho(ji,jj,jl)=0. 
     1041                        pfv_ho(ji,jj,jl)=0. 
     1042                     ENDIF 
    14871043                     ! 
    14881044                  ENDIF 
     
    14921048         CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond. 
    14931049 
    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  
    15231050      ENDIF 
    15241051 
     
    15261053      ! Search local extrema 
    15271054      ! -------------------- 
    1528       ! max/min of pt & pt_low with large negative/positive value (-/+zbig) outside ice cover 
     1055      ! max/min of pt & pt_ups with large negative/positive value (-/+zbig) outside ice cover 
    15291056      z1_dt = 1._wp / pdt 
    15301057      DO jl = 1, jpl 
     
    15321059         DO jj = 1, jpj 
    15331060            DO ji = 1, jpi 
    1534                IF    ( pt(ji,jj,jl) <= 0._wp .AND. pt_low(ji,jj,jl) <= 0._wp ) THEN 
     1061               IF    ( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN 
    15351062                  zbup(ji,jj) = -zbig 
    15361063                  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 
     1064               ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) > 0._wp ) THEN 
     1065                  zbup(ji,jj) = pt_ups(ji,jj,jl) 
     1066                  zbdo(ji,jj) = pt_ups(ji,jj,jl) 
     1067               ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN 
    15411068                  zbup(ji,jj) = pt(ji,jj,jl) 
    15421069                  zbdo(ji,jj) = pt(ji,jj,jl) 
    15431070               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) ) 
     1071                  zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_ups(ji,jj,jl) ) 
     1072                  zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_ups(ji,jj,jl) ) 
    15461073               ENDIF 
    15471074            END DO 
     
    15511078            DO ji = fs_2, fs_jpim1   ! vector opt. 
    15521079               ! 
    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 
     1080               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 
     1081               zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1) ) 
    15631082               ! 
    15641083               zpos = MAX( 0., pfu_ho(ji-1,jj,jl) ) - MIN( 0., pfu_ho(ji  ,jj,jl) ) &  ! positive/negative part of the flux 
     
    15671086                  & + MAX( 0., pfv_ho(ji,jj  ,jl) ) - MIN( 0., pfv_ho(ji,jj-1,jl) ) 
    15681087               ! 
    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. 
     1088               zpos = zpos - ( pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) & 
     1089                  &          ) * ( 1. - pamsk ) 
     1090               zneg = zneg + ( pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) & 
     1091                  &          ) * ( 1. - pamsk ) 
     1092               ! 
     1093               !                                  ! up & down beta terms 
     1094               IF( zpos > 0. ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt 
     1095               ELSE                 ; zbetup(ji,jj,jl) = 0. ! zbig 
    15761096               ENDIF 
    15771097               ! 
    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 
     1098               IF( zneg > 0. ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 
     1099               ELSE                 ; zbetdo(ji,jj,jl) = 0. ! zbig 
    15851100               ENDIF 
    15861101               ! 
     
    15891104               IF( zdo ==  zbig )   zbetdo(ji,jj,jl) = 0. ! zbig             
    15901105               ! 
    1591  
    1592 !!            IF( ji==26 .AND. jj==86) THEN 
    1593 !               WRITE(numout,*) '-----------------' 
    1594 !               WRITE(numout,*) 'zpos',zpos,zpos2 
    1595 !               WRITE(numout,*) 'zneg',zneg,zneg2 
    1596 !               WRITE(numout,*) 'puc/pu',ABS(puc(ji,jj))/MAX(epsi20, ABS(pu(ji,jj))) 
    1597 !               WRITE(numout,*) 'pvc/pv',ABS(pvc(ji,jj))/MAX(epsi20, ABS(pv(ji,jj))) 
    1598 !               WRITE(numout,*) 'pucm1/pu',ABS(puc(ji-1,jj))/MAX(epsi20, ABS(pu(ji-1,jj))) 
    1599 !               WRITE(numout,*) 'pvcm1/pv',ABS(pvc(ji,jj-1))/MAX(epsi20, ABS(pv(ji,jj-1))) 
    1600 !               WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)+pfu_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt 
    1601 !               WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)+pfv_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt 
    1602 !               WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)+pfu_low(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt 
    1603 !               WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)+pfv_low(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt 
    1604 !               WRITE(numout,*) 'pfu_low',pfu_low(ji,jj) * r1_e1e2t(ji,jj) * pdt 
    1605 !               WRITE(numout,*) 'pfv_low',pfv_low(ji,jj) * r1_e1e2t(ji,jj) * pdt 
    1606 !               WRITE(numout,*) 'pfu_lowm1',pfu_low(ji-1,jj) * r1_e1e2t(ji,jj) * pdt 
    1607 !               WRITE(numout,*) 'pfv_lowm1',pfv_low(ji,jj-1) * r1_e1e2t(ji,jj) * pdt 
    1608 !                
    1609 !               WRITE(numout,*) 'pt',pt(ji,jj) 
    1610 !               WRITE(numout,*) 'ptim1',pt(ji-1,jj) 
    1611 !               WRITE(numout,*) 'ptjm1',pt(ji,jj-1) 
    1612 !               WRITE(numout,*) 'pt_low',pt_low(ji,jj) 
    1613 !               WRITE(numout,*) 'zbetup',zbetup(ji,jj) 
    1614 !               WRITE(numout,*) 'zbetdo',zbetdo(ji,jj) 
    1615 !               WRITE(numout,*) 'zup',zup 
    1616 !               WRITE(numout,*) 'zdo',zdo 
    1617 !            ENDIF 
    1618             ! 
    16191106            END DO 
    16201107         END DO 
     
    16331120               ! 
    16341121               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  
     1122               ! 
     1123               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_ups(ji,jj,jl) 
     1124               ! 
    16441125            END DO 
    16451126         END DO 
     
    16521133               ! 
    16531134               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 
     1135               ! 
     1136               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_ups(ji,jj,jl) 
     1137               ! 
    16621138            END DO 
    16631139         END DO 
    16641140 
    16651141         ! 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 
     1142!!         DO jj = 2, jpjm1 
     1143!!            DO ji = 2, fs_jpim1   ! vector opt. 
     1144!!               zzt = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1145!!                  &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
     1146!!                  &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1147!!                  &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
     1148!!                  &         ) * tmask(ji,jj,1) 
     1149!!               IF( zzt < -epsi20 ) THEN 
     1150!!                  WRITE(numout,*) 'T<0 nonosc',zzt 
     1151!!               ENDIF 
     1152!!            END DO 
     1153!!         END DO 
    16831154 
    16841155      END DO 
     
    16881159   END SUBROUTINE nonosc_2d 
    16891160 
    1690    SUBROUTINE limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 
     1161   SUBROUTINE limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 
    16911162      !!--------------------------------------------------------------------- 
    16921163      !!                    ***  ROUTINE limiter_x  *** 
     
    16961167      REAL(wp)                   , INTENT(in   )           ::   pdt          ! tracer time-step 
    16971168      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 
    16991169      REAL(wp), DIMENSION (:,:,:), INTENT(in   )           ::   pt           ! ice tracer 
    17001170      REAL(wp), DIMENSION (:,:,:), INTENT(inout)           ::   pfu_ho       ! high order flux 
     
    17471217                  ELSE 
    17481218                     Cr = 0. 
    1749                      !IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm * 1.e20 
    1750                      !ELSE                        ;   Cr = Rjp * 1.e20 
    1751                      !ENDIF 
    17521219                  ENDIF 
    17531220 
     
    17871254   END SUBROUTINE limiter_x 
    17881255 
    1789    SUBROUTINE limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 
     1256   SUBROUTINE limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 
    17901257      !!--------------------------------------------------------------------- 
    17911258      !!                    ***  ROUTINE limiter_y  *** 
     
    17951262      REAL(wp)                   , INTENT(in   )           ::   pdt          ! tracer time-step 
    17961263      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 
    17981264      REAL(wp), DIMENSION (:,:,:), INTENT(in   )           ::   pt           ! ice tracer 
    17991265      REAL(wp), DIMENSION (:,:,:), INTENT(inout)           ::   pfv_ho       ! high order flux 
     
    18471313                  ELSE 
    18481314                     Cr = 0. 
    1849                      !IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm * 1.e20 
    1850                      !ELSE                        ;   Cr = Rjp * 1.e20 
    1851                      !ENDIF 
    18521315                  ENDIF 
    18531316 
Note: See TracChangeset for help on using the changeset viewer.