Index: /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/cfgs/SHARED/field_def_nemo-oce.xml	(revision 9975)
+++ /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/cfgs/SHARED/field_def_nemo-oce.xml	(revision 9976)
@@ -78,4 +78,8 @@
         <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" />
         <!-- next variables available with key_diahth -->
         <field id="mlddzt"       long_name="Thermocline Depth (depth of max dT/dz)"         standard_name="depth_at_maximum_upward_derivative_of_sea_water_potential_temperature"             unit="m"                         />
Index: /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/cfgs/SHARED/namelist_ref	(revision 9975)
+++ /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/cfgs/SHARED/namelist_ref	(revision 9976)
@@ -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/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DOM/dom_oce.F90
===================================================================
--- /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DOM/dom_oce.F90	(revision 9975)
+++ /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DOM/dom_oce.F90	(revision 9976)
@@ -135,6 +135,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3w_0 ,   e3w_b ,   e3w_n            !: w- vert. scale factor [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3uw_0 ,  e3uw_b ,  e3uw_n , e3uw_a   !: uw-vert. scale factor [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3vw_0 ,  e3vw_b ,  e3vw_n , e3vw_a   !: vw-vert. scale factor [m]
 
    !                                                        !  ref.   ! before  !   now   !
@@ -267,5 +267,6 @@
          &      e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) ,         &
          &      e3uw_b(jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) ,         &               
-         &      e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) ,     STAT=ierr(5) )                       
+         &      e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) ,         &               
+         &      e3uw_a(jpi,jpj,jpk) , e3vw_a(jpi,jpj,jpk) ,     STAT=ierr(5) )                       
          !
       ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) ,                                           &
Index: /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DOM/domvvl.F90
===================================================================
--- /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DOM/domvvl.F90	(revision 9975)
+++ /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DOM/domvvl.F90	(revision 9976)
@@ -149,7 +149,9 @@
 
       ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
-      e3t_a(:,:,:) = e3t_n(:,:,:)
-      e3u_a(:,:,:) = e3u_n(:,:,:)
-      e3v_a(:,:,:) = e3v_n(:,:,:)
+      e3t_a (:,:,:) = e3t_n (:,:,:)
+      e3u_a (:,:,:) = e3u_n (:,:,:)
+      e3v_a (:,:,:) = e3v_n (:,:,:)
+      e3uw_a(:,:,:) = e3uw_n(:,:,:)
+      e3vw_a(:,:,:) = e3vw_n(:,:,:)
       !
       !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
@@ -534,10 +536,12 @@
       END IF
 
-      ! *********************************** !
-      ! After scale factors at u- v- points !
-      ! *********************************** !
+      ! ******************************************* !
+      ! After scale factors at u- v- uw- vw- points !
+      ! ******************************************* !
 
       CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' )
       CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' )
+      CALL dom_vvl_interpol( e3u_a(:,:,:), e3uw_a(:,:,:), 'UW' )
+      CALL dom_vvl_interpol( e3v_a(:,:,:), e3vw_a(:,:,:), 'VW' )
 
       ! *********************************** !
Index: /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DYN/dynzdf.F90
===================================================================
--- /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DYN/dynzdf.F90	(revision 9975)
+++ /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DYN/dynzdf.F90	(revision 9976)
@@ -71,4 +71,8 @@
       REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars
       REAL(wp) ::   zzws, ze3va        !   -      -
+      REAL(wp) ::   z1_e3un, z1_e3vn   !   -      -
+      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,34 +159,110 @@
       !                    !* 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
+      IF( ln_zad_Aimp ) THEN   !!
+         IF( ln_dynadv_vec ) THEN      !==  Vector invariant advection  ==!
+            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.
+                        z1_e3un = 1._wp / e3u_n(ji,jj,jk)
+                        zzwi = ( ( avm   (ji+1,jj,jk  ) + avm   (ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
+                           &     / e3uw_a(ji  ,jj,jk  ) ) * z1_e3un * wumask(ji,jj,jk  )
+                        zzws = ( ( avm   (ji+1,jj,jk+1) + avm   (ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
+                           &     / e3uw_a(ji  ,jj,jk+1) ) * z1_e3un * wumask(ji,jj,jk+1)
+                        zWu  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji+1,jj,jk  )   &
+                           &               + wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1)   )
+                        zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWu, 0._wp ) * z1_e3un )
+                        zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWu, 0._wp ) * z1_e3un )
+                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWu, 0._wp ) + MIN( zWu, 0._wp ) ) * z1_e3un )
+                     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.
+                        z1_e3un = 1._wp / e3u_n(ji,jj,jk)
+                        zzwi = ( ( avm   (ji+1,jj,jk  ) + avm(ji,jj,jk  ) )   &
+                           &     / e3uw_a(ji  ,jj,jk  ) ) * z1_e3un * wumask(ji,jj,jk  )
+                        zzws = ( ( avm   (ji+1,jj,jk+1) + avm(ji,jj,jk+1) )   &
+                           &     / e3uw_a(ji  ,jj,jk+1) ) * z1_e3un * wumask(ji,jj,jk+1)
+                        zWu  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji+1,jj,jk  )   &
+                           &               + wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1)   )
+                        zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWu, 0._wp ) * z1_e3un )
+                        zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWu, 0._wp ) * z1_e3un )
+                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWu, 0._wp ) + MIN( zWu, 0._wp ) ) * z1_e3un )
+                     END DO
+                  END DO
+               END DO
+            END SELECT
+         ELSE                          !==  Flux form advection  ==!
+            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 = ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
+                           &         / ( ze3ua * e3uw_a(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
+                        zzws = ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
+                           &         / ( ze3ua * e3uw_a(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) = - zdt * ( zzwi +  MIN( zWui, 0._wp ) ) 
+                        zws(ji,jj,jk) = - zdt * ( zzws -  MAX( zWus, 0._wp ) )
+                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - 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 = ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_a(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
+                        zzws = ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_a(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) = - zdt * ( zzwi +  MIN( zWui, 0._wp ) ) 
+                        zws(ji,jj,jk) = - zdt * ( zzws -  MAX( zWus, 0._wp ) )
+                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWui, 0._wp ) + MIN( zWus, 0._wp ) )
+                     END DO
+                  END DO
+               END DO
+            END SELECT
+         ENDIF
+      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_a(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_a(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_a(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
+                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_a(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
+      ENDIF
       !
       DO jj = 2, jpjm1     !* Surface boundary conditions
@@ -274,34 +354,110 @@
       !                       !* 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
+      IF( ln_zad_Aimp ) THEN   !!
+         IF( ln_dynadv_vec ) THEN      !==  Vector invariant advection  ==!
+            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.
+                        z1_e3vn = 1._wp / e3v_n(ji,jj,jk)
+                        zzwi = ( ( avm   (ji,jj+1,jk  ) + avm   (ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
+                           &     / e3vw_a(ji,jj  ,jk  ) ) * z1_e3vn * wvmask(ji,jj,jk  )
+                        zzws = ( ( avm   (ji,jj+1,jk+1) + avm   (ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
+                           &     / e3vw_a(ji,jj  ,jk+1) ) * z1_e3vn * wvmask(ji,jj,jk+1)
+                        zWv  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji,jj+1,jk  )   &
+                           &               + wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1)   )
+                        zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWv, 0._wp ) * z1_e3vn )
+                        zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWv, 0._wp ) * z1_e3vn )
+                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWv, 0._wp ) + MIN( zWv, 0._wp ) ) * z1_e3vn )
+                     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.
+                        z1_e3vn = 1._wp / e3v_n(ji,jj,jk)
+                        zzwi = ( ( avm   (ji,jj+1,jk  ) + avm(ji,jj,jk  ) )   &
+                           &     / e3vw_a(ji,jj  ,jk  ) ) * z1_e3vn * wvmask(ji,jj,jk  )
+                        zzws = ( ( avm   (ji,jj+1,jk+1) + avm(ji,jj,jk+1) )   &
+                           &     / e3vw_a(ji  ,jj,jk+1) ) * z1_e3vn * wvmask(ji,jj,jk+1)
+                        zWv  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji,jj+1,jk  )   &
+                           &               + wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1)   )
+                        zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWv, 0._wp ) * z1_e3vn )
+                        zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWv, 0._wp ) * z1_e3vn )
+                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWv, 0._wp ) + MIN( zWv, 0._wp ) ) * z1_e3vn )
+                     END DO
+                  END DO
+               END DO
+            END SELECT
+         ELSE                          !==  Flux form advection  ==!
+            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 U-point
+                        zzwi = ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
+                           &         / ( ze3va * e3vw_a(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
+                        zzws = ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
+                           &         / ( ze3va * e3vw_a(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                        zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) )
+                        zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) )
+                        zwi(ji,jj,jk) = - zdt * ( zzwi +  MIN( zWvi, 0._wp ) ) 
+                        zws(ji,jj,jk) = - zdt * ( zzws -  MAX( zWvs, 0._wp ) )
+                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - 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 U-point
+                        zzwi = ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_a(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
+                        zzws = ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_a(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                        zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) )
+                        zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) )
+                        zwi(ji,jj,jk) = - zdt * ( zzwi +  MIN( zWvi, 0._wp ) ) 
+                        zws(ji,jj,jk) = - zdt * ( zzws -  MAX( zWvs, 0._wp ) )
+                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
+                     END DO
+                  END DO
+               END DO
+            END SELECT
+         ENDIF
+      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 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
+      ENDIF
       !
       DO jj = 2, jpjm1        !* Surface boundary conditions
Index: /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DYN/sshwzv.F90
===================================================================
--- /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DYN/sshwzv.F90	(revision 9975)
+++ /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DYN/sshwzv.F90	(revision 9976)
@@ -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,96 @@
    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)             ::   zcff_min, zcff_max                    ! local scalars
+      REAL(wp) , PARAMETER ::   Cu_min = 0.6_wp                       ! local parameters
+      REAL(wp) , PARAMETER ::   Cu_max = 0.95_wp                      ! 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,*) '~~~~~ '
+         !
+         wi    (:,:,jpk) = 0._wp              ! bottom value : wi=0 (set once for all)
+         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 ) )   &
+                  &                        / e2t(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 ) )   &
+                  &                        / e1t(ji,jj)                                                       &
+                  &                      ) * z1_e3w
+            END DO
+         END DO
+      END DO
+      !
+      CALL iom_put("Courant",Cu_adv)
+      !
+      wi(:,:,:) = 0._wp
+      zcff_min=999. ; zcff_max = -999.
+      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                                  !<-- Fully implicit
+                     zcff = ( zCu - Cu_max )/ zCu
+                  ENDIF
+                  !
+                  !zcff = 0._wp
+                  wi(ji,jj,jk) =           zcff   * wn(ji,jj,jk)
+                  wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk)
+                  zcff_min = MIN( zcff, zcff_min )
+                  zcff_max = MAX( zcff, zcff_max )
+                  !
+               END DO
+            END DO
+         END DO
+      ELSE
+         ! Fully explicit everywhere
+         write(*,*) 'NoWi ',kt
+      ENDIF
+         write(*,*) 'ZCFF ',kt,zcff_min,zcff_max
+      CALL iom_put("wimp",wi)
+      CALL iom_put("wexp",wn)
+      !
+      IF( ln_timing )   CALL timing_stop('wAimp')
+      !
+   END SUBROUTINE wAimp
    !!======================================================================
 END MODULE sshwzv
Index: /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/TRA/trazdf.F90
===================================================================
--- /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/TRA/trazdf.F90	(revision 9975)
+++ /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/TRA/trazdf.F90	(revision 9976)
@@ -177,14 +177,26 @@
             !
             ! 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.
+                        zwi(ji,jj,jk) = - p2dt * ( zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  )  + MIN( wi(ji,jj,jk  ) , 0._wp ) )
+                        zws(ji,jj,jk) = - p2dt * ( zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1)  - MAX( wi(ji,jj,jk+1) , 0._wp ) )
+                        zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk)   &
+                           &                   - p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( 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/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/ZDF/zdf_oce.F90
===================================================================
--- /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/ZDF/zdf_oce.F90	(revision 9975)
+++ /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/ZDF/zdf_oce.F90	(revision 9976)
@@ -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/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/ZDF/zdfosm.F90
===================================================================
--- /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/ZDF/zdfosm.F90	(revision 9975)
+++ /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/ZDF/zdfosm.F90	(revision 9976)
@@ -32,5 +32,5 @@
 
    !!----------------------------------------------------------------------
-   !!   'key_zdfosm'                                             OSMOSIS scheme
+   !!   'ln_zdfosm'                                             OSMOSIS scheme
    !!----------------------------------------------------------------------
    !!   zdf_osm       : update momentum and tracer Kz from osm scheme
Index: /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/ZDF/zdfphy.F90
===================================================================
--- /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/ZDF/zdfphy.F90	(revision 9975)
+++ /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/ZDF/zdfphy.F90	(revision 9976)
@@ -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/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/oce.F90
===================================================================
--- /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/oce.F90	(revision 9975)
+++ /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/oce.F90	(revision 9976)
@@ -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/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/step.F90
===================================================================
--- /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/step.F90	(revision 9975)
+++ /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/step.F90	(revision 9976)
@@ -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
       
