Index: /NEMO/branches/2019/fix_vvl_ticket1791/tests/VORTEX/MY_SRC/domvvl.F90
===================================================================
--- /NEMO/branches/2019/fix_vvl_ticket1791/tests/VORTEX/MY_SRC/domvvl.F90	(revision 11011)
+++ /NEMO/branches/2019/fix_vvl_ticket1791/tests/VORTEX/MY_SRC/domvvl.F90	(revision 11012)
@@ -139,5 +139,5 @@
       CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )    ! from T to V 
       CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )    ! from U to F
+      CALL dom_vvl_interpol( e3t_n(:,:,:), e3f_n(:,:,:), 'F' )    ! from U to F
       !                                ! Vertical interpolation of e3t,u,v 
       CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )  ! from T to W
@@ -629,5 +629,5 @@
       ! - JC - hu_b, hv_b, hur_b, hvr_b also
       
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  )
+      CALL dom_vvl_interpol( e3t_n(:,:,:), e3f_n(:,:,:), 'F'  )
       
       ! Vertical scale factor interpolations
@@ -691,7 +691,19 @@
       !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
       !
-      INTEGER ::   ji, jj, jk                                       ! dummy loop indices
+      INTEGER ::   ji, jj, jk, jkbot                                ! dummy loop indices
+      INTEGER ::   nmet                                             ! horizontal interpolation method
       REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F
-      !!----------------------------------------------------------------------
+      REAL(wp) ::  zmin, zsmall, zfractap                           ! Minimum thicknesses UVF-points
+      REAL(wp) ::  zdo, zup                                         ! Lower and upper interfaces depths anomalies
+      REAL(wp), DIMENSION(jpi,jpj) :: zs                            ! Surface interface depth anomaly
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw                        ! Interface depth anomaly
+      !!----------------------------------------------------------------------
+      !
+!      nmet = 0  ! Original method (Surely wrong)
+      nmet = 1  ! Interface interpolation
+!      nmet = 2  ! Internal interfaces interpolation only, spread barotropic increment
+      !
+      zsmall = 1.e-10    ! Minimum thickness at U or V points (m)
+      zfractap = 0.8_wp  ! Fraction of T-point thickness in the shallowest direction
       !
       IF(ln_wd_il) THEN
@@ -701,43 +713,238 @@
       END IF
       !
+      IF ( ( nmet==1 ).OR.( nmet==2 ) ) THEN
+         SELECT CASE ( pout )
+            !
+         CASE( 'U', 'V', 'F' )
+            !
+            ! Retrieve ssh:
+            zs(:,:) = 0._wp
+            DO jk = 1, jpkm1
+                zs(:,:) = zs(:,:) + pe3_in(:,:,jk)*tmask(:,:,jk)
+            END DO 
+            zs(:,:) = zs(:,:) - ht_0(:,:)
+            !
+            ! Interface depth:
+            zw(:,:,:) =  0._wp
+            DO jk=2,jpk
+               zw(:,:,jk) = zw(:,:,jk-1) + pe3_in(:,:,jk-1)*tmask(:,:,jk-1)
+            END DO 
+            !
+            DO jk=1,jpkm1
+               zw(:,:,jk) = zw(:,:,jk) - zs(:,:) *tmask(:,:,jk)
+            END DO
+            zw(:,:,jpk) = gdepw_0(:,:,jpk)
+            !
+            IF ( nmet==2 ) THEN        ! Consider "internal" interfaces only
+               !
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     DO jk=1,jpk
+                        zw(ji,jj,jk) = ((zw(ji,jj,jk) + zs(ji,jj)-gdepw_0(ji,jj,mikt(ji,jj)))              &
+                                     & * ht_0(ji,jj) / (ht_0(ji,jj) + zs(ji,jj) + 1._wp - ssmask(ji,jj))  &
+                                     &  + gdepw_0(ji,jj,mikt(ji,jj)) )* tmask(ji,jj,jk)
+                     END DO
+                  END DO
+               END DO 
+            ENDIF 
+            !
+         END SELECT
+      END IF
+      !
+      pe3_out(:,:,:) = 0.0_wp
+      !
       SELECT CASE ( pout )    !==  type of interpolation  ==!
          !
       CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
-         DO jk = 1, jpk
+         IF (nmet==0) THEN
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
+                        &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
+                        &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
+                  END DO
+               END DO
+            END DO
+         ELSE
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
-                     &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
-                     &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
-               END DO
-            END DO
-         END DO
+                  jkbot = mbku(ji,jj)
+                  zdo = hu_0(ji,jj)
+                  DO jk=jkbot,1,-1
+                     zup =   0.5_wp * ( e1e2t(ji  ,jj)*(zw(ji  ,jj,jk)-gdepw_0(ji  ,jj,jk))                     &
+                         &            + e1e2t(ji+1,jj)*(zw(ji+1,jj,jk)-gdepw_0(ji+1,jj,jk)) ) * r1_e1e2u(ji,jj) &
+                         & + 0.5_wp * (gdepw_0(ji  ,jj,jk)+gdepw_0(ji+1,jj,jk))
+                     !
+                     ! If there is a step, taper bottom interface:
+                     zmin = 0._wp
+                     IF (zup > zdo) THEN                     
+                        IF ( ht_0(ji+1,jj) < ht_0(ji,jj) ) THEN
+                           zmin = zfractap * (zw(ji+1,jj,jk+1)-zw(ji+1,jj,jk))
+                        ELSE
+                           zmin = zfractap * (zw(ji  ,jj,jk+1)-zw(ji  ,jj,jk))
+                        ENDIF
+                     ENDIF       
+                     zmin = MAX(zmin, zsmall)        
+                     zup  = MIN(zup, zdo-zmin)
+                     pe3_out(ji,jj,jk) = (zdo - zup) * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
+                                       &               - umask(ji,jj,jk) * e3u_0(ji,jj,jk) 
+                     zdo = zup
+                  END DO
+               END DO
+            END DO
+         END IF
+         ! 
+         IF (nmet==2) THEN           ! Spread sea level anomaly
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  DO jk=1,jpk
+                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                   &
+                           &               + ( pe3_out(ji,jj,jk) + e3u_0(ji,jj,jk) )               & 
+                           &               / ( hu_0(ji,jj) + 1._wp - ssumask(ji,jj) )              &
+                           &               * 0.5_wp * r1_e1e2u(ji,jj)                              &
+                           &               * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )        &
+                           &               * ( e1e2t(ji,jj)*zs(ji,jj) + e1e2t(ji+1,jj)*zs(ji+1,jj) )       
+                  END DO
+               END DO
+            END DO
+            !
+         ENDIF
+         !
          CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp )
          pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
          !
       CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
-         DO jk = 1, jpk
+         IF (nmet==0) THEN
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
+                        &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
+                        &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
+                  END DO
+               END DO
+            END DO
+         ELSE
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
-                     &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
-                     &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
-               END DO
-            END DO
-         END DO
+                  jkbot = mbkv(ji,jj)
+                  zdo = hv_0(ji,jj)
+                  DO jk=jkbot,1,-1
+                     zup =   0.5_wp * ( e1e2t(ji  ,jj)*(zw(ji,jj  ,jk)-gdepw_0(ji,jj  ,jk))                     &
+                         &            + e1e2t(ji+1,jj)*(zw(ji,jj+1,jk)-gdepw_0(ji,jj+1,jk)) ) * r1_e1e2v(ji,jj) &
+                         & + 0.5_wp * (gdepw_0(ji,jj  ,jk)+gdepw_0(ji,jj+1,jk))
+                     !
+                     ! If there is a step, taper bottom interface:
+                     zmin = 0._wp
+                     IF (zup > zdo) THEN                     
+                        IF ( ht_0(ji,jj+1) < ht_0(ji,jj) ) THEN
+                           zmin = zfractap * (zw(ji,jj+1,jk+1)-zw(ji,jj+1,jk))
+                        ELSE
+                           zmin = zfractap * (zw(ji,jj  ,jk+1)-zw(ji,jj  ,jk))
+                        ENDIF
+                     ENDIF 
+                     zmin = MAX(zmin, zsmall)        
+                     zup  = MIN(zup, zdo-zmin)
+                     pe3_out(ji,jj,jk) = (zdo - zup) * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
+                                       &               - vmask(ji,jj,jk) * e3v_0(ji,jj,jk) 
+                     zdo = zup
+                  END DO
+               END DO
+            END DO
+         END IF
+         !
+         IF (nmet==2) THEN           ! Spread sea level anomaly
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  DO jk=1,jpk
+                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                                       &
+                           &               + ( pe3_out(ji,jj,jk) + e3v_0(ji,jj,jk) )                                   & 
+                           &               / ( hv_0(ji,jj) + 1._wp - ssvmask(ji,jj) )                                  &
+                           &               * 0.5_wp * r1_e1e2v(ji,jj)                                                  &
+                                           * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )                            &
+                           &               * ( e1e2t(ji,jj)*zs(ji,jj) + e1e2t(ji,jj+1)*zs(ji,jj+1) )       
+                  END DO
+               END DO
+            END DO
+            !
+         ENDIF
+         !
          CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp )
          pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
          !
-      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
-         DO jk = 1, jpk
+      CASE( 'F' )                   !* from T-point to F-point : hor. surface weighted mean
+         IF (nmet==0) THEN
+            DO jk=1,jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     pe3_out(ji,jj,jk) = 0.25_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
+                           &                     *    r1_e1e2f(ji,jj)                                                  &
+                           &                     * (  e1e2t(ji  ,jj  ) * ( pe3_in(ji  ,jj  ,jk)-e3t_0(ji  ,jj  ,jk) )  & 
+                           &                        + e1e2t(ji  ,jj+1) * ( pe3_in(ji  ,jj+1,jk)-e3t_0(ji  ,jj+1,jk) )  &
+                           &                        + e1e2t(ji+1,jj  ) * ( pe3_in(ji+1,jj  ,jk)-e3t_0(ji+1,jj  ,jk) )  & 
+                           &                        + e1e2t(ji+1,jj+1) * ( pe3_in(ji+1,jj+1,jk)-e3t_0(ji+1,jj+1,jk) ) )                                                  
+                  END DO
+               END DO
+            END DO
+         ELSE
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
-                     &                       *    r1_e1e2f(ji,jj)                                                  &
-                     &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
-                     &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
-               END DO
-            END DO
-         END DO
+                  jkbot = MIN(mbku(ji,jj), mbku(ji,jj+1)) 
+                  zdo = hf_0(ji,jj)
+                  DO jk=jkbot,1,-1
+                     zup =     0.25_wp * (  e1e2t(ji  ,jj  ) * (zw(ji  ,jj  ,jk)-gdepw_0(ji  ,jj  ,jk))                         & 
+                           &              + e1e2t(ji+1,jj  ) * (zw(ji+1,jj  ,jk)-gdepw_0(ji+1,jj  ,jk))                         & 
+                           &              + e1e2t(ji  ,jj+1) * (zw(ji  ,jj+1,jk)-gdepw_0(ji  ,jj+1,jk))                         & 
+                           &              + e1e2t(ji+1,jj+1) * (zw(ji+1,jj+1,jk)-gdepw_0(ji+1,jj+1,jk))  ) *    r1_e1e2f(ji,jj) &
+                           & + 0.25_wp * ( gdepw_0(ji  ,jj  ,jk) + gdepw_0(ji+1,jj  ,jk) &
+                           &              +gdepw_0(ji  ,jj+1,jk) + gdepw_0(ji+1,jj+1,jk) )
+                     !
+                     ! If there is a step, taper bottom interface:
+                     zmin = 0._wp
+                     IF (zup > zdo) THEN 
+                        IF ( hu_0(ji,jj+1) < hu_0(ji,jj) ) THEN
+                           IF ( ht_0(ji+1,jj+1) < ht_0(ji  ,jj+1) ) THEN
+                              zmin = zfractap * (zw(ji+1,jj+1,jk+1)-zw(ji+1,jj+1,jk))
+                           ELSE
+                              zmin = zfractap * (zw(ji  ,jj+1,jk+1)-zw(ji  ,jj+1,jk))
+                           ENDIF
+                        ELSE
+                           IF ( ht_0(ji+1,jj  ) < ht_0(ji  ,jj  ) ) THEN
+                              zmin = zfractap * (zw(ji+1,jj  ,jk+1)-zw(ji+1,jj  ,jk))
+                           ELSE
+                              zmin = zfractap * (zw(ji  ,jj  ,jk+1)-zw(ji  ,jj  ,jk))
+                           ENDIF
+                        ENDIF
+                     ENDIF
+                     zmin = MAX(zmin, zsmall)        
+                     zup  = MIN(zup, zdo-zmin)
+                     !
+                     pe3_out(ji,jj,jk) = ( zdo - zup ) & 
+                                      & *( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
+                                      &  - umask(ji,jj,jk) * umask(ji,jj+1,jk) * e3f_0(ji,jj,jk) 
+                     zdo = zup
+                  END DO
+               END DO
+            END DO
+         END IF
+         !
+         IF (nmet==2) THEN           ! Spread sea level anomaly
+            !
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  DO jk=1,jpk
+                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                           &
+                           &               + ( pe3_out(ji,jj,jk) + e3f_0(ji,jj,jk) )                       & 
+                           &               / ( hf_0(ji,jj) + 1._wp - ssumask(ji,jj)*ssumask(ji,jj+1) )     &
+                           &               * 0.25_wp * r1_e1e2f(ji,jj)                                     & 
+                           &               * ( umask(ji,jj,jk)*umask(ji,jj+1,jk)*(1.0_wp - zlnwd) + zlnwd )&
+                           &               * ( e1e2t(ji  ,jj)*zs(ji  ,jj) + e1e2t(ji  ,jj+1)*zs(ji  ,jj+1) &
+                           &                  +e1e2t(ji+1,jj)*zs(ji+1,jj) + e1e2t(ji+1,jj+1)*zs(ji+1,jj+1) )               
+                  END DO
+               END DO
+            END DO
+         END IF
+         !
          CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp )
          pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
@@ -929,5 +1136,5 @@
                !
                ! usr_def_istate called here only to get sshb, that is needed to initialize e3t_b and e3t_n
-               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  )  
+               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  )
                ! usr_def_istate will be called again in istate_init to initialize ts(bn), ssh(bn), u(bn) and v(bn)
                !
