Index: /NEMO/branches/2018/dev_r9947_SI3_advection/tests/ICEADV/MY_SRC/icedyn_adv_umx.F90
===================================================================
--- /NEMO/branches/2018/dev_r9947_SI3_advection/tests/ICEADV/MY_SRC/icedyn_adv_umx.F90	(revision 10280)
+++ /NEMO/branches/2018/dev_r9947_SI3_advection/tests/ICEADV/MY_SRC/icedyn_adv_umx.F90	(revision 10281)
@@ -35,4 +35,19 @@
    REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120
 
+   REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z1_a
+   
+   ! alternate directions for upstream
+   ! clem: it gives results for Lipscomb test that are the same as "ll_upsxy=false"
+   LOGICAL :: ll_upsxy = .FALSE.
+   
+   ! prelimiter
+   LOGICAL :: ll_prelimiter = .TRUE.
+   LOGICAL :: ll_prelimiter_zalesak = .TRUE.   ! from: Zalesak(1979) eq. 14 => better than Devore
+   LOGICAL :: ll_prelimiter_devore  = .FALSE.  ! from: Devore eq. 11
+
+   ! iterate on the limiter (only nonosc)
+   LOGICAL :: ll_limiter_it2 = .TRUE.
+   
+
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -77,5 +92,5 @@
       !                                     !    must be >= 0.01 and the best seems to be 0.1
       REAL(wp), DIMENSION(jpi,jpj) ::   zudy, zvdx, zcu_box, zcv_box, zfu, zfv
-      REAL(wp), DIMENSION(jpi,jpj) ::   z1_a, zh_i, zh_s, zs_i, zo_i, ze_i, ze_s, z1_ap, zh_ip 
+      REAL(wp), DIMENSION(jpi,jpj) ::   zh_i, zh_s, zs_i, zo_i, ze_i, ze_s, zh_ip 
       !!----------------------------------------------------------------------
       !
@@ -112,4 +127,5 @@
       END DO
 
+      IF(.NOT. ALLOCATED(z1_a )) ALLOCATE(z1_a (jpi,jpj))
       !---------------!
       !== advection ==!
@@ -149,10 +165,10 @@
             END DO
             !
-            IF ( ln_pnd_H12 ) THEN   ! melt ponds
+            IF ( ln_pnd_H12 ) THEN   ! melt ponds (must be the last ones to be advected because of z1_a...)
                ! to avoid a problem with the limiter nonosc when A gets close to 0
                pa_ip(:,:,jl) = pa_ip(:,:,jl) + zeps * tmask(:,:,1)
                !
-               WHERE( pa_ip(:,:,jl) > epsi20 )   ;   z1_ap(:,:) = 1._wp / pa_ip(:,:,jl)
-               ELSEWHERE                         ;   z1_ap(:,:) = 0.
+               WHERE( pa_ip(:,:,jl) > epsi20 )   ;   z1_a(:,:) = 1._wp / pa_ip(:,:,jl)
+               ELSEWHERE                         ;   z1_a(:,:) = 0.
                END WHERE
                !
@@ -213,26 +229,73 @@
       INTEGER  ::   ji, jj           ! dummy loop indices  
       REAL(wp) ::   ztra             ! local scalar
-!!clem      REAL(wp) ::   zeps = 1.e-02        ! local scalar
       INTEGER  ::   kn_limiter = 1   ! 1=nonosc ; 2=superbee ; 3=h3(rachid)
-      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ho , zfv_ho , zt_u, zt_v, zpt
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ho , zfv_ho , zt_u, zt_v, zpt, zz1_a
       REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfv_ups, zt_ups   ! only for nonosc 
       !!----------------------------------------------------------------------
       !
-      ! add a constant value to avoid problems with zeros
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zpt(ji,jj) = pt(ji,jj) !!clem + zeps * tmask(ji,jj,1)
-         END DO
-      END DO
-
       !  upstream (_ups) advection with initial mass fluxes
       ! ---------------------------------------------------
-      ! fluxes
-      DO jj = 1, jpjm1
-         DO ji = 1, fs_jpim1
-            zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * zpt(ji+1,jj)
-            zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * zpt(ji,jj+1)
-         END DO
-      END DO
+      IF( .NOT. ll_upsxy ) THEN
+
+         ! fluxes
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1
+               zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj)
+               zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1)
+            END DO
+         END DO
+
+      ELSE
+         !
+         IF    ( pamsk == 1. ) THEN ; zz1_a(:,:) = 1._wp     ;        ! 1   if advection of A
+         ELSEIF( pamsk == 0. ) THEN ; zz1_a(:,:) = z1_a(:,:) ; ENDIF  ! 1/A if advection of V
+         !
+         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
+            ! flux in x-direction
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1
+                  zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj)
+               END DO
+            END DO
+            ! first guess of tracer content from u-flux
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) * zz1_a(ji,jj)
+               END DO
+            END DO
+            CALL lbc_lnk( zpt, 'T', 1. )
+            !
+            ! flux in y-direction
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1
+                  zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * zpt(ji,jj+1)
+               END DO
+            END DO
+            !
+         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==!
+            ! flux in y-direction
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1
+                  zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1)
+               END DO
+            END DO
+            ! first guess of tracer content from v-flux
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zpt(ji,jj) = ( ptc(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) * zz1_a(ji,jj) 
+               END DO
+            END DO
+            CALL lbc_lnk( zpt, 'T', 1. )
+            !
+            ! flux in x-direction
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1
+                  zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * zpt(ji+1,jj)
+               END DO
+            END DO
+            !
+         ENDIF
+         
+      ENDIF
       ! guess after content field with upstream scheme
       DO jj = 2, jpjm1
@@ -251,10 +314,10 @@
       CASE ( 20 )                          !== centered second order ==!
          !
-         CALL cen2( kn_limiter, jt, kt, pdt, zpt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho,  &
+         CALL cen2( kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho,  &
             &       zt_ups, zfu_ups, zfv_ups )
          !
       CASE ( 1:5 )                         !== 1st to 5th order ULTIMATE-MACHO scheme ==!
          !
-         CALL macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, zpt, pu, pv, puc, pvc, pubox, pvbox, ptc, zt_u, zt_v, zfu_ho, zfv_ho,  &
+         CALL macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, zt_u, zt_v, zfu_ho, zfv_ho,  &
             &        zt_ups, zfu_ups, zfv_ups )
          !
@@ -266,6 +329,6 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1
-               pfu(ji,jj) =  zfu_ho(ji,jj) !!clem - zeps * puc(ji,jj)
-               pfv(ji,jj) =  zfv_ho(ji,jj) !!clem - zeps * pvc(ji,jj)
+               pfu(ji,jj) =  zfu_ho(ji,jj)
+               pfv(ji,jj) =  zfv_ho(ji,jj)
             END DO
          END DO
@@ -278,6 +341,5 @@
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1 
-            ztra       = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) )         &  !        Div(uaH) or        Div(ua)
-!!clem               &           + ( puc   (ji,jj) - puc   (ji-1,jj) + pvc   (ji,jj) - pvc   (ji,jj-1) ) * zeps  &  ! epsi * Div(ua)  or epsi * Div(u) 
+            ztra       = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) )  &  !        Div(uaH) or        Div(ua)
                &         ) * r1_e1e2t(ji,jj)  
             ptc(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
@@ -430,5 +492,6 @@
       !
       INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp), DIMENSION(jpi,jpj) ::   zzt
+      REAL(wp) ::   ztra
+      REAL(wp), DIMENSION(jpi,jpj) ::   zzt, zzfu_ho, zzfv_ho
       !!----------------------------------------------------------------------
       !
@@ -478,5 +541,24 @@
          !
       ENDIF
-      IF( kn_limiter == 1 )   CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+      IF( kn_limiter == 1 ) THEN
+         IF( .NOT. ll_limiter_it2 ) THEN
+            CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+         ELSE
+            zzfu_ho(:,:) = pfu_ho(:,:)
+            zzfv_ho(:,:) = pfv_ho(:,:)
+            ! 1st iteration of nonosc (limit the flux with the upstream solution)
+            CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho )
+            ! guess after content field with high order
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  ztra       = - ( zzfu_ho(ji,jj) - zzfu_ho(ji-1,jj) + zzfv_ho(ji,jj) - zzfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj)
+                  zzt(ji,jj) =   ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
+               END DO
+            END DO
+            CALL lbc_lnk( zzt, 'T', 1. )
+            ! 2nd iteration of nonosc (limit the flux with the limited high order solution)
+            CALL nonosc_2d ( pdt, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho )
+         ENDIF
+      ENDIF
       !
    END SUBROUTINE macho
@@ -733,13 +815,13 @@
      
 
-   SUBROUTINE nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )
+   SUBROUTINE nonosc_2d( pdt, ptc, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE nonosc  ***
       !!     
-      !! **  Purpose :   compute monotonic tracer fluxes from the upstream 
+       !! **  Purpose :   compute monotonic tracer fluxes from the upstream 
       !!       scheme and the before field by a nonoscillatory algorithm 
       !!
       !! **  Method  :   ... ???
-      !!       warning : ptc and pt_ups must be masked, but the boundaries
+      !!       warning : ptc and pt_low must be masked, but the boundaries
       !!       conditions on the fluxes are not necessary zalezak (1979)
       !!       drange (1995) multi-dimensional forward-in-time and upstream-
@@ -747,12 +829,12 @@
       !!----------------------------------------------------------------------
       REAL(wp)                     , INTENT(in   ) ::   pdt              ! tracer time-step
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   ptc, pt_ups      ! before field & upstream guess of after field
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pfv_ups, pfu_ups ! upstream flux
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   ptc, pt_low      ! before field & upstream guess of after field
+      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pfv_low, pfu_low ! upstream flux
       REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   pfv_ho, pfu_ho   ! monotonic flux
       !
       INTEGER  ::   ji, jj    ! dummy loop indices
       REAL(wp) ::   zpos, zneg, zbig, zsml, z1_dt   ! local scalars
-      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo         !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zdiv
+      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign        !   -      -
+      REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zdiv, zti_low, ztj_low
       !!----------------------------------------------------------------------
       zbig = 1.e+40_wp
@@ -771,20 +853,74 @@
       DO jj = 1, jpjm1
          DO ji = 1, fs_jpim1   ! vector opt.
-            pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_ups(ji,jj)
-            pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_ups(ji,jj)
-         END DO
-      END DO
-
+            pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_low(ji,jj)
+            pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_low(ji,jj)
+         END DO
+      END DO
+
+      ! extreme case where pfu_ho has to be zero
+      ! ----------------------------------------
+      !                                    pfu_ho
+      !                           *         --->
+      !                        |      |  *   |        | 
+      !                        |      |      |    *   |    
+      !                        |      |      |        |    *
+      !            t_low :       i-1     i       i+1       i+2   
+      IF( ll_prelimiter ) THEN
+
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1 
+               zti_low(ji,jj)= pt_low(ji+1,jj  )
+               ztj_low(ji,jj)= pt_low(ji  ,jj+1)
+            END DO
+         END DO
+         CALL lbc_lnk_multi( zti_low, 'T', 1., ztj_low, 'T', 1. )
+
+         IF( ll_prelimiter_zalesak ) THEN
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  IF(    SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji+1,jj) - pt_low (ji  ,jj) ) ) THEN
+                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_low(ji+1,jj) - zti_low(ji  ,jj) ) )   pfu_ho(ji,jj) = 0.
+                     IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji  ,jj) - pt_low (ji-1,jj) ) )   pfu_ho(ji,jj) = 0.
+                  ENDIF
+                  IF(    SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji,jj+1) - pt_low (ji,jj  ) ) ) THEN
+                     IF( SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_low(ji,jj+1) - ztj_low(ji,jj  ) ) )   pfv_ho(ji,jj) = 0.
+                     IF( SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji,jj  ) - pt_low (ji,jj-1) ) )   pfv_ho(ji,jj) = 0.
+                  ENDIF
+               END DO
+            END DO
+            CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond.
+
+         ELSEIF( ll_prelimiter_devore ) THEN
+            z1_dt = 1._wp / pdt
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  zsign = SIGN( 1., pt_low(ji+1,jj) - pt_low(ji,jj) )
+                  pfu_ho(ji,jj) =  zsign * MAX( 0. , MIN( ABS(pfu_ho(ji,jj)) , &
+                     &                          zsign * ( pt_low (ji  ,jj) - pt_low (ji-1,jj) ) * e1e2t(ji  ,jj) * z1_dt , &
+                     &                          zsign * ( zti_low(ji+1,jj) - zti_low(ji  ,jj) ) * e1e2t(ji+1,jj) * z1_dt ) )
+                  
+                  zsign = SIGN( 1., pt_low(ji,jj+1) - pt_low(ji,jj) )
+                  pfv_ho(ji,jj) =  zsign * MAX( 0. , MIN( ABS(pfv_ho(ji,jj)) , &
+                     &                          zsign * ( pt_low (ji,jj  ) - pt_low (ji,jj-1) ) * e1e2t(ji,jj  ) * z1_dt , &
+                     &                          zsign * ( ztj_low(ji,jj+1) - ztj_low(ji,jj  ) ) * e1e2t(ji,jj+1) * z1_dt ) )
+               END DO
+            END DO
+            CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond.
+            
+         ENDIF
+         
+      ENDIF
+      
       ! Search local extrema
       ! --------------------
-      ! max/min of ptc & pt_ups with large negative/positive value (-/+zbig) outside ice cover
+      ! max/min of ptc & pt_low with large negative/positive value (-/+zbig) outside ice cover
       DO jj = 1, jpj
-         DO ji = fs_2, fs_jpim1
-            IF( ptc(ji,jj) == 0._wp .AND. pt_ups(ji,jj) == 0._wp .AND. zdiv(ji,jj) == 0._wp ) THEN
+         DO ji = 1, jpi
+            IF( ptc(ji,jj) == 0._wp .AND. pt_low(ji,jj) == 0._wp .AND. zdiv(ji,jj) == 0._wp ) THEN
                zbup(ji,jj) = -zbig
                zbdo(ji,jj) =  zbig
             ELSE               
-               zbup(ji,jj) = MAX( ptc(ji,jj) , pt_ups(ji,jj) )
-               zbdo(ji,jj) = MIN( ptc(ji,jj) , pt_ups(ji,jj) )
+               zbup(ji,jj) = MAX( ptc(ji,jj) , pt_low(ji,jj) )
+               zbdo(ji,jj) = MIN( ptc(ji,jj) , pt_low(ji,jj) )
             ENDIF
          END DO
@@ -795,25 +931,26 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             !
-            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
-            zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj), zbdo(ji+1,jj), zbdo(ji,jj-1), zbdo(ji,jj+1) )
-            !
-            zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji  ,jj) ) &
-               & + MAX( 0., pfv_ho(ji,jj-1) ) - MIN( 0., pfv_ho(ji,jj  ) )  ! positive/negative part of the flux
+            !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
+            !zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1) )
+            !
+            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
+               &                     zbup(ji-1,jj-1), zbup(ji+1,jj+1), zbup(ji+1,jj-1), zbup(ji-1,jj+1)  )
+            zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj  ), zbdo(ji+1,jj  ), zbdo(ji  ,jj-1), zbdo(ji  ,jj+1), &
+               &                     zbdo(ji-1,jj-1), zbdo(ji+1,jj+1), zbdo(ji+1,jj-1), zbdo(ji-1,jj+1)  )
+            !
+            zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji  ,jj) ) &  ! positive/negative part of the flux
+               & + MAX( 0., pfv_ho(ji,jj-1) ) - MIN( 0., pfv_ho(ji,jj  ) )
             zneg = MAX( 0., pfu_ho(ji  ,jj) ) - MIN( 0., pfu_ho(ji-1,jj) ) &
                & + MAX( 0., pfv_ho(ji,jj  ) ) - MIN( 0., pfv_ho(ji,jj-1) )
             !
             !                                  ! up & down beta terms
-!!clem            zbetup(ji,jj) = ( zup - pt_ups(ji,jj) ) / ( zpos + zsml ) * e1e2t(ji,jj) * z1_dt
-!!clem            zbetdo(ji,jj) = ( pt_ups(ji,jj) - zdo ) / ( zneg + zsml ) * e1e2t(ji,jj) * z1_dt
-            IF( zpos >= epsi20 ) THEN
-               zbetup(ji,jj) = ( zup - pt_ups(ji,jj) ) / zpos * e1e2t(ji,jj) * z1_dt
-            ELSE
-               zbetup(ji,jj) = zbig
+!!clem            zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / ( zpos + zsml ) * e1e2t(ji,jj) * z1_dt
+!!clem            zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / ( zneg + zsml ) * e1e2t(ji,jj) * z1_dt
+            IF( zpos >= epsi20 ) THEN ; zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / zpos * e1e2t(ji,jj) * z1_dt
+            ELSE                      ; zbetup(ji,jj) = zbig
             ENDIF
             !
-            IF( zneg >= epsi20 ) THEN
-               zbetdo(ji,jj) = ( pt_ups(ji,jj) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
-            ELSE
-               zbetdo(ji,jj) = zbig
+            IF( zneg >= epsi20 ) THEN ; zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt
+            ELSE                      ; zbetdo(ji,jj) = zbig
             ENDIF
             !
@@ -830,5 +967,5 @@
             zcu = 0.5  + SIGN( 0.5 , pfu_ho(ji,jj) )
             !
-            pfu_ho(ji,jj) = pfu_ho(ji,jj) * ( zcu * zau + ( 1._wp - zcu ) * zbu ) + pfu_ups(ji,jj)
+            pfu_ho(ji,jj) = pfu_ho(ji,jj) * ( zcu * zau + ( 1._wp - zcu ) * zbu ) + pfu_low(ji,jj)
          END DO
       END DO
@@ -840,5 +977,5 @@
             zcv = 0.5  + SIGN( 0.5 , pfv_ho(ji,jj) )
             !
-            pfv_ho(ji,jj) = pfv_ho(ji,jj) * ( zcv * zav + ( 1._wp - zcv ) * zbv ) + pfv_ups(ji,jj)
+            pfv_ho(ji,jj) = pfv_ho(ji,jj) * ( zcv * zav + ( 1._wp - zcv ) * zbv ) + pfv_low(ji,jj)
          END DO
       END DO
