Index: /NEMO/trunk/src/OCE/DYN/sshwzv.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/sshwzv.F90	(revision 11413)
+++ /NEMO/trunk/src/OCE/DYN/sshwzv.F90	(revision 11414)
@@ -9,4 +9,5 @@
    !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
    !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work
+   !!            4.0  !  2018-12  (A. Coward) add mixed implicit/explicit advection
    !!----------------------------------------------------------------------
 
@@ -279,5 +280,7 @@
       !!            :   wi      : now vertical velocity (for implicit treatment)
       !!
-      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
+      !! Reference  : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent
+      !!              implicit scheme for vertical advection in oceanic modeling. 
+      !!              Ocean Modelling, 91, 38-69.
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! time step
@@ -300,20 +303,40 @@
       ENDIF
       !
-      !
-      DO jk = 1, jpkm1            ! calculate Courant numbers
-         DO jj = 2, jpjm1
-            DO ji = 2, fs_jpim1   ! vector opt.
-               z1_e3t = 1._wp / e3t_n(ji,jj,jk)
-               Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )   &  ! 2*rdt and not r2dt (for restartability)
-                  &                             + ( MAX( e2u(ji  ,jj)*e3u_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   &
-                  &                                 MIN( e2u(ji-1,jj)*e3u_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   &
-                  &                               * r1_e1e2t(ji,jj)                                                 &
-                  &                             + ( MAX( e1v(ji,jj  )*e3v_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   &
-                  &                                 MIN( e1v(ji,jj-1)*e3v_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   &
-                  &                               * r1_e1e2t(ji,jj)                                                 &
-                  &                             ) * z1_e3t
+      ! Calculate Courant numbers
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = 2, fs_jpim1   ! vector opt.
+                  z1_e3t = 1._wp / e3t_n(ji,jj,jk)
+                  ! 2*rdt and not r2dt (for restartability)
+                  Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )                       &  
+                     &                             + ( MAX( e2u(ji  ,jj)*e3u_n(ji  ,jj,jk)*un(ji  ,jj,jk) + un_td(ji  ,jj,jk), 0._wp ) -   &
+                     &                                 MIN( e2u(ji-1,jj)*e3u_n(ji-1,jj,jk)*un(ji-1,jj,jk) + un_td(ji-1,jj,jk), 0._wp ) )   &
+                     &                               * r1_e1e2t(ji,jj)                                                                     &
+                     &                             + ( MAX( e1v(ji,jj  )*e3v_n(ji,jj  ,jk)*vn(ji,jj  ,jk) + vn_td(ji,jj  ,jk), 0._wp ) -   &
+                     &                                 MIN( e1v(ji,jj-1)*e3v_n(ji,jj-1,jk)*vn(ji,jj-1,jk) + vn_td(ji,jj-1,jk), 0._wp ) )   &
+                     &                               * r1_e1e2t(ji,jj)                                                                     &
+                     &                             ) * z1_e3t
+               END DO
             END DO
          END DO
-      END DO
+      ELSE
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = 2, fs_jpim1   ! vector opt.
+                  z1_e3t = 1._wp / e3t_n(ji,jj,jk)
+                  ! 2*rdt and not r2dt (for restartability)
+                  Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )   & 
+                     &                             + ( MAX( e2u(ji  ,jj)*e3u_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   &
+                     &                                 MIN( e2u(ji-1,jj)*e3u_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   &
+                     &                               * r1_e1e2t(ji,jj)                                                 &
+                     &                             + ( MAX( e1v(ji,jj  )*e3v_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   &
+                     &                                 MIN( e1v(ji,jj-1)*e3v_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   &
+                     &                               * r1_e1e2t(ji,jj)                                                 &
+                     &                             ) * z1_e3t
+               END DO
+            END DO
+         END DO
+      ENDIF
       CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1. )
       !
