Index: /NEMO/trunk/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- /NEMO/trunk/cfgs/SHARED/field_def_nemo-oce.xml	(revision 10363)
+++ /NEMO/trunk/cfgs/SHARED/field_def_nemo-oce.xml	(revision 10364)
@@ -77,4 +77,10 @@
         <!-- variables available with MLE -->
         <field id="Lf_NHpf"      long_name="MLE: Lf = N H / f"   unit="m" />
+
+        <!-- next variables available with ln_zad_Aimp=.true. -->
+        <field id="Courant"      long_name="Courant number"                                                                 unit="#"   grid_ref="grid_T_3D" />
+        <field id="wimp"         long_name="Implicit vertical velocity"                                                     unit="m/s" grid_ref="grid_T_3D" />
+        <field id="wexp"         long_name="Explicit vertical velocity"                                                     unit="m/s" grid_ref="grid_T_3D" />
+        <field id="wi_cff"       long_name="Fraction of implicit vertical velocity"                                         unit="#"   grid_ref="grid_T_3D" />
 
         <!-- next variables available with key_diahth -->
Index: /NEMO/trunk/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/trunk/cfgs/SHARED/namelist_ref	(revision 10363)
+++ /NEMO/trunk/cfgs/SHARED/namelist_ref	(revision 10364)
@@ -972,4 +972,7 @@
 &namzdf        !   vertical physics manager                             (default: NO selection)
 !-----------------------------------------------------------------------
+   !                       ! adaptive-implicit vertical advection
+   ln_zad_Aimp = .false.      !  Courant number dependent scheme (Shchepetkin 2015)
+   !
    !                       ! type of vertical closure (required)
    ln_zdfcst   = .false.      !  constant mixing
Index: /NEMO/trunk/src/OCE/DYN/dynzdf.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/dynzdf.F90	(revision 10363)
+++ /NEMO/trunk/src/OCE/DYN/dynzdf.F90	(revision 10364)
@@ -71,4 +71,8 @@
       REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars
       REAL(wp) ::   zzws, ze3va        !   -      -
+      REAL(wp) ::   z1_e3ua, z1_e3va   !   -      -
+      REAL(wp) ::   zWu , zWv          !   -      -
+      REAL(wp) ::   zWui, zWvi         !   -      -
+      REAL(wp) ::   zWus, zWvs         !   -      -
       REAL(wp), DIMENSION(jpi,jpj,jpk)        ::  zwi, zwd, zws   ! 3D workspace 
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv   !  -      -
@@ -155,41 +159,88 @@
       !                    !* Matrix construction
       zdt = r2dt * 0.5
-      SELECT CASE( nldf_dyn )
-      CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1 
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point
-                  zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
-                     &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
-                  zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
-                     &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
-                  zwi(ji,jj,jk) = zzwi
-                  zws(ji,jj,jk) = zzws
-                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
-               END DO
-            END DO
-         END DO
-      CASE DEFAULT               ! iso-level lateral mixing
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1 
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point
-                  zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
-                  zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
-                  zwi(ji,jj,jk) = zzwi
-                  zws(ji,jj,jk) = zzws
-                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
-               END DO
-            END DO
-         END DO
-      END SELECT
-      !
-      DO jj = 2, jpjm1     !* Surface boundary conditions
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            zwi(ji,jj,1) = 0._wp
-            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
-         END DO
-      END DO
+      IF( ln_zad_Aimp ) THEN   !!
+         SELECT CASE( nldf_dyn )
+         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1 
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
+                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
+                        &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
+                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
+                        &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+                     zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) )
+                     zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
+                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 
+                     zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
+                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
+                  END DO
+               END DO
+            END DO
+         CASE DEFAULT               ! iso-level lateral mixing
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1 
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
+                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
+                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+                     zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) )
+                     zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
+                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )
+                     zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
+                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
+                  END DO
+               END DO
+            END DO
+         END SELECT
+         DO jj = 2, jpjm1     !* Surface boundary conditions
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zwi(ji,jj,1) = 0._wp
+               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)
+               zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw_n(ji,jj,2) ) * wumask(ji,jj,2)
+               zWus = 0.5_wp * ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) )
+               zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp )
+               zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) )
+            END DO
+         END DO
+      ELSE
+         SELECT CASE( nldf_dyn )
+         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1 
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
+                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
+                        &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
+                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
+                        &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+                     zwi(ji,jj,jk) = zzwi
+                     zws(ji,jj,jk) = zzws
+                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
+                  END DO
+               END DO
+            END DO
+         CASE DEFAULT               ! iso-level lateral mixing
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1 
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
+                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
+                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+                     zwi(ji,jj,jk) = zzwi
+                     zws(ji,jj,jk) = zzws
+                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
+                  END DO
+               END DO
+            END DO
+         END SELECT
+         DO jj = 2, jpjm1     !* Surface boundary conditions
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zwi(ji,jj,1) = 0._wp
+               zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
+            END DO
+         END DO
+      ENDIF
+      !
       !
       !              !==  Apply semi-implicit bottom friction  ==!
@@ -274,41 +325,88 @@
       !                       !* Matrix construction
       zdt = r2dt * 0.5
-      SELECT CASE( nldf_dyn )
-      CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1   
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point
-                  zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
-                     &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
-                  zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
-                     &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
-                  zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
-                  zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
-                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
-               END DO
-            END DO
-         END DO
-      CASE DEFAULT               ! iso-level lateral mixing
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1   
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point
-                  zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
-                  zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
-                  zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
-                  zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
-                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
-               END DO
-            END DO
-         END DO
-      END SELECT
-      !
-      DO jj = 2, jpjm1        !* Surface boundary conditions
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            zwi(ji,jj,1) = 0._wp
-            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
-         END DO
-      END DO
+      IF( ln_zad_Aimp ) THEN   !!
+         SELECT CASE( nldf_dyn )
+         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv)
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1 
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
+                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
+                        &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
+                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
+                        &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                     zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  )
+                     zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1)
+                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp )
+                     zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp )
+                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
+                  END DO
+               END DO
+            END DO
+         CASE DEFAULT               ! iso-level lateral mixing
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1 
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
+                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
+                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                     zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  )
+                     zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1)
+                     zwi(ji,jj,jk) = zzwi  + zdt * MIN( zWvi, 0._wp )
+                     zws(ji,jj,jk) = zzws  - zdt * MAX( zWvs, 0._wp )
+                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
+                  END DO
+               END DO
+            END DO
+         END SELECT
+         DO jj = 2, jpjm1     !* Surface boundary conditions
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zwi(ji,jj,1) = 0._wp
+               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)
+               zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw_n(ji,jj,2) ) * wvmask(ji,jj,2)
+               zWvs = 0.5_wp * ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) )
+               zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp )
+               zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) )
+            END DO
+         END DO
+      ELSE
+         SELECT CASE( nldf_dyn )
+         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1   
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
+                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
+                        &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
+                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
+                        &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                     zwi(ji,jj,jk) = zzwi
+                     zws(ji,jj,jk) = zzws
+                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
+                  END DO
+               END DO
+            END DO
+         CASE DEFAULT               ! iso-level lateral mixing
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1   
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
+                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
+                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                     zwi(ji,jj,jk) = zzwi
+                     zws(ji,jj,jk) = zzws
+                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
+                  END DO
+               END DO
+            END DO
+         END SELECT
+         DO jj = 2, jpjm1        !* Surface boundary conditions
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zwi(ji,jj,1) = 0._wp
+               zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
+            END DO
+         END DO
+      ENDIF
+      !
       !              !==  Apply semi-implicit top/bottom friction  ==!
       !
Index: /NEMO/trunk/src/OCE/DYN/sshwzv.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/sshwzv.F90	(revision 10363)
+++ /NEMO/trunk/src/OCE/DYN/sshwzv.F90	(revision 10364)
@@ -28,4 +28,5 @@
 #endif
    !
+   USE iom 
    USE in_out_manager ! I/O manager
    USE restart        ! only for lrst_oce
@@ -41,4 +42,5 @@
    PUBLIC   ssh_nxt    ! called by step.F90
    PUBLIC   wzv        ! called by step.F90
+   PUBLIC   wAimp      ! called by step.F90
    PUBLIC   ssh_swp    ! called by step.F90
 
@@ -265,4 +267,92 @@
    END SUBROUTINE ssh_swp
 
+   SUBROUTINE wAimp( kt )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE wAimp  ***
+      !!                   
+      !! ** Purpose :   compute the Courant number and partition vertical velocity
+      !!                if a proportion needs to be treated implicitly
+      !!
+      !! ** Method  : - 
+      !!
+      !! ** action  :   wn      : now vertical velocity (to be handled explicitly)
+      !!            :   wi      : now vertical velocity (for implicit treatment)
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! time step
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp)             ::   zCu, zcff, z1_e3w                     ! local scalars
+      REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters
+      REAL(wp) , PARAMETER ::   Cu_max = 0.27                         ! local parameters
+      REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters
+      REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('wAimp')
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
+         IF(lwp) WRITE(numout,*) '~~~~~ '
+         !
+         Cu_adv(:,:,jpk) = 0._wp              ! bottom value : Cu_adv=0 (set once for all)
+      ENDIF
+      !
+      DO jk = 1, jpkm1            ! calculate Courant numbers
+         DO jj = 2, jpjm1
+            DO ji = 2, fs_jpim1   ! vector opt.
+               z1_e3w = 1._wp / e3w_n(ji,jj,jk)
+               Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )    &
+                  &                      + ( MAX( e2u(ji  ,jj)*e3uw_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   &
+                  &                          MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   &
+                  &                        * r1_e1e2t(ji,jj)                                                  &
+                  &                      + ( MAX( e1v(ji,jj  )*e3vw_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   &
+                  &                          MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   &
+                  &                        * r1_e1e2t(ji,jj)                                                  &
+                  &                      ) * z1_e3w
+            END DO
+         END DO
+      END DO
+      !
+      CALL iom_put("Courant",Cu_adv)
+      !
+      wi(:,:,:) = 0._wp                                 ! Includes top and bottom values set to zero
+      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere
+         DO jk = 1, jpkm1                               ! or scan Courant criterion and partition
+            DO jj = 2, jpjm1                            ! w where necessary
+               DO ji = 2, fs_jpim1   ! vector opt.
+                  !
+                  zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) )
+                  !
+                  IF( zCu < Cu_min ) THEN               !<-- Fully explicit
+                     zcff = 0._wp
+                  ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit
+                     zcff = ( zCu - Cu_min )**2
+                     zcff = zcff / ( Fcu + zcff )
+                  ELSE                                  !<-- Mostly implicit
+                     zcff = ( zCu - Cu_max )/ zCu
+                  ENDIF
+                  zcff = MIN(1._wp, zcff)
+                  !
+                  wi(ji,jj,jk) =           zcff   * wn(ji,jj,jk)
+                  wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk)
+                  !
+                  Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient
+               END DO
+            END DO
+         END DO
+      ELSE
+         ! Fully explicit everywhere
+         Cu_adv = 0.0_wp                                ! Reuse array to output coefficient
+      ENDIF
+      CALL iom_put("wimp",wi) 
+      CALL iom_put("wi_cff",Cu_adv)
+      CALL iom_put("wexp",wn)
+      !
+      IF( ln_timing )   CALL timing_stop('wAimp')
+      !
+   END SUBROUTINE wAimp
    !!======================================================================
 END MODULE sshwzv
Index: /NEMO/trunk/src/OCE/TRA/trazdf.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/trazdf.F90	(revision 10363)
+++ /NEMO/trunk/src/OCE/TRA/trazdf.F90	(revision 10364)
@@ -136,5 +136,5 @@
       !
       INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp) ::  zrhs             ! local scalars
+      REAL(wp) ::  zrhs, zzwi, zzws ! local scalars
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwi, zwt, zwd, zws
       !!---------------------------------------------------------------------
@@ -177,14 +177,28 @@
             !
             ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked)
-            DO jk = 1, jpkm1
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1   ! vector opt.
-!!gm BUG  I think, use e3w_a instead of e3w_n, not sure of that
-                     zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  )
-                     zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1)
-                     zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk)
-                 END DO
-               END DO
-            END DO
+            IF( ln_zad_Aimp ) THEN         ! Adaptive implicit vertical advection
+               DO jk = 1, jpkm1
+                  DO jj = 2, jpjm1
+                     DO ji = fs_2, fs_jpim1   ! vector opt. (ensure same order of calculation as below if wi=0.)
+                        zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  )
+                        zzws = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1)
+                        zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zzwi - zzws   &
+                           &                 + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )
+                        zwi(ji,jj,jk) = zzwi + p2dt *   MIN( wi(ji,jj,jk  ) , 0._wp )
+                        zws(ji,jj,jk) = zzws - p2dt *   MAX( wi(ji,jj,jk+1) , 0._wp )
+                    END DO
+                  END DO
+               END DO
+            ELSE
+               DO jk = 1, jpkm1
+                  DO jj = 2, jpjm1
+                     DO ji = fs_2, fs_jpim1   ! vector opt.
+                        zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk)
+                        zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1)
+                        zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk)
+                    END DO
+                  END DO
+               END DO
+            ENDIF
             !
             !! Matrix inversion from the first level
Index: /NEMO/trunk/src/OCE/ZDF/zdf_oce.F90
===================================================================
--- /NEMO/trunk/src/OCE/ZDF/zdf_oce.F90	(revision 10363)
+++ /NEMO/trunk/src/OCE/ZDF/zdf_oce.F90	(revision 10364)
@@ -18,4 +18,6 @@
 
    !                            !!* namelist namzdf: vertical physics *
+   !                             ! Adaptive-implicit vertical advection flag
+   LOGICAL , PUBLIC ::   ln_zad_Aimp !: adaptive (Courant number-based) implicit vertical advection
    !                             ! vertical closure scheme flags
    LOGICAL , PUBLIC ::   ln_zdfcst   !: constant coefficients
Index: /NEMO/trunk/src/OCE/ZDF/zdfosm.F90
===================================================================
--- /NEMO/trunk/src/OCE/ZDF/zdfosm.F90	(revision 10363)
+++ /NEMO/trunk/src/OCE/ZDF/zdfosm.F90	(revision 10364)
@@ -32,5 +32,5 @@
 
    !!----------------------------------------------------------------------
-   !!   'key_zdfosm'                                             OSMOSIS scheme
+   !!   'ln_zdfosm'                                             OSMOSIS scheme
    !!----------------------------------------------------------------------
    !!   zdf_osm       : update momentum and tracer Kz from osm scheme
Index: /NEMO/trunk/src/OCE/ZDF/zdfphy.F90
===================================================================
--- /NEMO/trunk/src/OCE/ZDF/zdfphy.F90	(revision 10363)
+++ /NEMO/trunk/src/OCE/ZDF/zdfphy.F90	(revision 10364)
@@ -80,4 +80,5 @@
          &             ln_zdfswm,                                    &     ! surface  wave-induced mixing
          &             ln_zdfiwm,                                    &     ! internal  -      -      -
+         &             ln_zad_Aimp,                                  &     ! apdative-implicit vertical advection
          &             rn_avm0, rn_avt0, nn_avb, nn_havtb                  ! coefficients
       !!----------------------------------------------------------------------
@@ -101,4 +102,6 @@
       IF(lwp) THEN                      ! Parameter print
          WRITE(numout,*) '   Namelist namzdf : set vertical mixing mixing parameters'
+         WRITE(numout,*) '      adaptive-implicit vertical advection'
+         WRITE(numout,*) '         Courant number targeted application   ln_zad_Aimp = ', ln_zad_Aimp
          WRITE(numout,*) '      vertical closure scheme'
          WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdfcst = ', ln_zdfcst
@@ -127,4 +130,9 @@
       ENDIF
 
+      IF( ln_zad_Aimp ) THEN
+         IF( zdf_phy_alloc() /= 0 )   &
+        &       CALL ctl_stop( 'STOP', 'zdf_phy_init : unable to allocate adaptive-implicit z-advection arrays' )
+         wi(:,:,:) = 0._wp
+      ENDIF
       !                          !==  Background eddy viscosity and diffusivity  ==!
       IF( nn_avb == 0 ) THEN             ! Define avmb, avtb from namelist parameter
@@ -316,4 +324,13 @@
       !
    END SUBROUTINE zdf_phy
+   INTEGER FUNCTION zdf_phy_alloc()
+      !!----------------------------------------------------------------------
+      !!                 ***  FUNCTION zdf_phy_alloc  ***
+      !!----------------------------------------------------------------------
+     ! Allocate wi array (declared in oce.F90) for use with the adaptive-implicit vertical velocity option
+     ALLOCATE(     wi(jpi,jpj,jpk), Cu_adv(jpi,jpj,jpk),  STAT= zdf_phy_alloc )
+     IF( zdf_phy_alloc /= 0 )   CALL ctl_warn('zdf_phy_alloc: failed to allocate ln_zad_Aimp=T required arrays')
+     IF( lk_mpp             )   CALL mpp_sum ( zdf_phy_alloc )
+   END FUNCTION zdf_phy_alloc
 
    !!======================================================================
Index: /NEMO/trunk/src/OCE/oce.F90
===================================================================
--- /NEMO/trunk/src/OCE/oce.F90	(revision 10363)
+++ /NEMO/trunk/src/OCE/oce.F90	(revision 10364)
@@ -22,4 +22,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   vb   ,  vn    , va     !: j-horizontal velocity        [m/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           wn             !: vertical velocity            [m/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           wi             !: vertical vel. (adaptive-implicit) [m/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           hdivn          !: horizontal divergence        [s-1]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   tsb  ,  tsn   , tsa    !: 4D T-S fields                  [Celsius,psu] 
@@ -29,4 +30,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rhd    !: in situ density anomalie rhd=(rho-rau0)/rau0  [no units]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rhop   !: potential volumic mass                           [kg/m3]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   Cu_adv                   !: vertical Courant number (adaptive-implicit)
 
    !! free surface                                      !  before  ! now    ! after  !
Index: /NEMO/trunk/src/OCE/step.F90
===================================================================
--- /NEMO/trunk/src/OCE/step.F90	(revision 10363)
+++ /NEMO/trunk/src/OCE/step.F90	(revision 10364)
@@ -160,6 +160,7 @@
 
                             CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_hor)
-      IF(.NOT.ln_linssh )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors 
+      IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors 
                             CALL wzv           ( kstp )  ! now cross-level velocity 
+      IF( ln_zad_Aimp )     CALL wAimp         ( kstp )  ! Adaptive-implicit vertical advection partitioning
                             CALL eos    ( tsn, rhd, rhop, gdept_n(:,:,:) )  ! now in situ density for hpg computation
                             
@@ -198,4 +199,5 @@
          IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component)
                             CALL wzv        ( kstp )              ! now cross-level velocity 
+         IF( ln_zad_Aimp )  CALL wAimp      ( kstp )  ! Adaptive-implicit vertical advection partitioning
       ENDIF
       
Index: /NEMO/trunk/src/OCE/stpctl.F90
===================================================================
--- /NEMO/trunk/src/OCE/stpctl.F90	(revision 10363)
+++ /NEMO/trunk/src/OCE/stpctl.F90	(revision 10364)
@@ -24,4 +24,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp         ! distributed memory computing
+   USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables
    USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy
 
@@ -32,5 +33,5 @@
    PUBLIC stp_ctl           ! routine called by step.F90
 
-   INTEGER  ::   idrun, idtime, idssh, idu, ids1, ids2, istatus
+   INTEGER  ::   idrun, idtime, idssh, idu, ids1, ids2, idt1, idt2, idc1, idw1, istatus
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -68,5 +69,5 @@
       INTEGER , DIMENSION(3) ::   ilocu, ilocs1, ilocs2
       INTEGER , DIMENSION(2) ::   iloch
-      REAL(wp), DIMENSION(5) ::   zmax
+      REAL(wp), DIMENSION(9) ::   zmax
       CHARACTER(len=20) :: clname
       !!----------------------------------------------------------------------
@@ -90,6 +91,13 @@
             istatus = NF90_DEF_VAR( idrun,       's_min', NF90_DOUBLE, (/ idtime /), ids1  )
             istatus = NF90_DEF_VAR( idrun,       's_max', NF90_DOUBLE, (/ idtime /), ids2  )
+            istatus = NF90_DEF_VAR( idrun,       't_min', NF90_DOUBLE, (/ idtime /), idt1  )
+            istatus = NF90_DEF_VAR( idrun,       't_max', NF90_DOUBLE, (/ idtime /), idt2  )
+            IF( ln_zad_Aimp ) THEN
+               istatus = NF90_DEF_VAR( idrun,   'abs_wi_max', NF90_DOUBLE, (/ idtime /), idw1  )
+               istatus = NF90_DEF_VAR( idrun,       'Cu_max', NF90_DOUBLE, (/ idtime /), idc1  )
+            ENDIF
             istatus = NF90_ENDDEF(idrun)
          ENDIF
+         zmax(8:9) = 0._wp    ! initialise to zero in case ln_zad_Aimp option is not in use
          
       ENDIF
@@ -109,10 +117,16 @@
       zmax(3) = MAXVAL( -tsn(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp )   ! minus salinity max
       zmax(4) = MAXVAL(  tsn(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp )   !       salinity max
-      zmax(5) = REAL( nstop , wp )                                            ! stop indicator
+      zmax(5) = MAXVAL( -tsn(:,:,:,jp_tem) , mask = tmask(:,:,:) == 1._wp )   ! minus temperature max
+      zmax(6) = MAXVAL(  tsn(:,:,:,jp_tem) , mask = tmask(:,:,:) == 1._wp )   !       temperature max
+      zmax(7) = REAL( nstop , wp )                                            ! stop indicator
+      IF( ln_zad_Aimp ) THEN
+         zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = wmask(:,:,:) == 1._wp ) ! implicit vertical vel. max
+         zmax(9) = MAXVAL(   Cu_adv(:,:,:)   , mask = tmask(:,:,:) == 1._wp ) !       cell Courant no. max
+      ENDIF
       !
       IF( lk_mpp ) THEN
-         CALL mpp_max_multiple( zmax(:), 5 )    ! max over the global domain
+         CALL mpp_max_multiple( zmax(:), 9 )    ! max over the global domain
          !
-         nstop = NINT( zmax(5) )                 ! nstop indicator sheared among all local domains
+         nstop = NINT( zmax(7) )                 ! nstop indicator sheared among all local domains
       ENDIF
       !
@@ -172,4 +186,10 @@
          istatus = NF90_PUT_VAR( idrun,  ids1, (/-zmax(3)/), (/kt/), (/1/) )
          istatus = NF90_PUT_VAR( idrun,  ids2, (/ zmax(4)/), (/kt/), (/1/) )
+         istatus = NF90_PUT_VAR( idrun,  idt1, (/-zmax(5)/), (/kt/), (/1/) )
+         istatus = NF90_PUT_VAR( idrun,  idt2, (/ zmax(6)/), (/kt/), (/1/) )
+         IF( ln_zad_Aimp ) THEN
+            istatus = NF90_PUT_VAR( idrun,  idw1, (/ zmax(8)/), (/kt/), (/1/) )
+            istatus = NF90_PUT_VAR( idrun,  idc1, (/ zmax(9)/), (/kt/), (/1/) )
+         ENDIF
          IF( MOD( kt , 100 ) == 0 ) istatus = NF90_SYNC(idrun)
          IF( kt == nitend         ) istatus = NF90_CLOSE(idrun)
