Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex	(revision 11406)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex	(revision 11407)
@@ -1329,8 +1329,8 @@
   \label{eq:Eqn_zad_Aimp_partition}
 Cu_{min} &= 0.15 \nonumber \\
-Cu_{max} &= 0.27 \nonumber \\
+Cu_{max} &= 0.3  \nonumber \\
 Cu_{cut} &= 2Cu_{max} - Cu_{min} \nonumber \\
 Fcu    &= 4Cu_{max}*(Cu_{max}-Cu_{min}) \nonumber \\
-Cf &=
+C\kern-0.14em f &=
      \begin{cases}
         0.0                                                        &\text{if $Cu \leq Cu_{min}$} \\
@@ -1340,6 +1340,4 @@
 \end{align}
 
-\noindent With these settings the coefficient ($Cf$) is shown in \autoref{fig:zad_Aimp_coeff}
-
 \begin{figure}[!t]
   \begin{center}
@@ -1347,5 +1345,5 @@
     \caption{
       \protect\label{fig:zad_Aimp_coeff}
-      The value of the partitioning coefficient ($Cf$) used to partition vertical velocities into parts to
+      The value of the partitioning coefficient ($C\kern-0.14em f$) used to partition vertical velocities into parts to
       be treated implicitly and explicitly for a range of typical Courant numbers (\forcode{ln_zad_Aimp=.true.})
     }
@@ -1359,22 +1357,38 @@
 \begin{align}
   \label{eq:Eqn_zad_Aimp_partition2}
-    w_{i_{ijk}} &= Cf_{ijk} w_{n_{ijk}}     \nonumber \\
-    w_{n_{ijk}} &= (1-Cf_{ijk}) w_{n_{ijk}}           
+    w_{i_{ijk}} &= C\kern-0.14em f_{ijk} w_{n_{ijk}}     \nonumber \\
+    w_{n_{ijk}} &= (1-C\kern-0.14em f_{ijk}) w_{n_{ijk}}           
 \end{align}
 
 \noindent Note that the coefficient is such that the treatment is never fully implicit;
 the three cases from \autoref{eq:Eqn_zad_Aimp_partition} can be considered as:
-fully-explicit; mixed explicit/implicit and mostly-implicit.  The $w_i$ component is added
-to the implicit solvers for the vertical mixing in \mdl{dynzdf} and \mdl{trazdf} in a
-similar way to \citep{shchepetkin_OM15}.  This is sufficient for the flux-limited
-advection scheme (\forcode{ln_traadv_mus}) but further intervention is required when using
-the flux-corrected scheme (\forcode{ln_traadv_fct}). For these schemes the implicit
-upstream fluxes must be added to both the monotonic guess and to the higher order solution
-when calculating the antidiffusive fluxes. The implicit vertical fluxes are then removed
-since they are added by the implicit solver later on.
+fully-explicit; mixed explicit/implicit and mostly-implicit.  With the settings shown the
+coefficient ($C\kern-0.14em f$) varies as shown in \autoref{fig:zad_Aimp_coeff}. Note with these values
+the $Cu_{cut}$ boundary between the mixed implicit-explicit treatment and 'mostly
+implicit' is 0.45 which is just below the stability limited given in
+\autoref{tab:zad_Aimp_CFLcrit}  for a 3rd order scheme.
+
+The $w_i$ component is added to the implicit solvers for the vertical mixing in
+\mdl{dynzdf} and \mdl{trazdf} in a similar way to \citep{shchepetkin_OM15}.  This is
+sufficient for the flux-limited advection scheme (\forcode{ln_traadv_mus}) but further
+intervention is required when using the flux-corrected scheme (\forcode{ln_traadv_fct}).
+For these schemes the implicit upstream fluxes must be added to both the monotonic guess
+and to the higher order solution when calculating the antidiffusive fluxes. The implicit
+vertical fluxes are then removed since they are added by the implicit solver later on.
 
 The adaptive-implicit vertical advection option is new to NEMO at v4.0 and has yet to be 
 used in a wide range of simulations. The following test simulation, however, does illustrate
 the potential benefits and will hopefully encourage further testing and feedback from users:
+
+\begin{figure}[!t]
+  \begin{center}
+    \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_overflow_frames}
+    \caption{
+      \protect\label{fig:zad_Aimp_overflow_frames}
+      A time-series of temperature vertical cross-sections for the OVERFLOW test case. These results are for the default
+      settings with \forcode{nn_rdt=10.0} and without adaptive implicit vertical advection (\forcode{ln_zad_Aimp=.false.}).
+    }
+  \end{center}
+\end{figure}
 
 \subsection{Adaptive-implicit vertical advection in the OVERFLOW test-case}
@@ -1389,6 +1403,6 @@
      ln_zdfnpc     = .true.
      ln_traadv_fct = .true.
-        nn_fct_h   =  4
-        nn_fct_v   =  4
+        nn_fct_h   =  2
+        nn_fct_v   =  2
 \end{verbatim}
 
@@ -1404,26 +1418,45 @@
 candidate, therefore, for use of the adaptive-implicit vertical advection scheme.
 
-The results with \forcode{ln_zad_Aimp=.true.} and a variety of model timesteps are shown
-in \autoref{fig:zad_Aimp_overflow_all_rdt} (together with the equivalent frames from the
-base run).  In this simple example the use of the adaptive-implicit vertcal advection
-scheme has enabled a 12x increase in the model timestep without significantly altering the
-solution (although at this extreme the plume is more diffuse and has not travelled so far).
-Notably, the solution with and without the scheme is slightly different even with
-\forcode{nn_rdt = 10.}; suggesting that the base run was close enough to instability to
-trigger the scheme despite completing successfully.  To assist in diagnosing how active
-the scheme is, in both location and time, the 3D implicit and explicit components of the
-vertical velocity are available via XIOS as \texttt{wimp} and \texttt{wexp} respectively.
-Likewise, the partitioning coefficient ($Cf$) is also available as \texttt{wi\_cff}.
-
-\begin{figure}[!t]
-  \begin{center}
-    \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_overflow_frames}
-    \caption{
-      \protect\label{fig:zad_Aimp_overflow_frames}
-      A time-series of temperature vertical cross-sections for the OVERFLOW test case. These results are for the default
-      settings with \forcode{nn_rdt=10.0} and without adaptive implicit vertical advection (\forcode{ln_zad_Aimp=.false.}).
-    }
-  \end{center}
-\end{figure}
+The results with \forcode{ln_zad_Aimp=.true.} and a variety of model timesteps
+are shown in \autoref{fig:zad_Aimp_overflow_all_rdt} (together with the equivalent
+frames from the base run).  In this simple example the use of the adaptive-implicit
+vertcal advection scheme has enabled a 12x increase in the model timestep without
+significantly altering the solution (although at this extreme the plume is more diffuse
+and has not travelled so far).  Notably, the solution with and without the scheme is
+slightly different even with \forcode{nn_rdt = 10.}; suggesting that the base run was
+close enough to instability to trigger the scheme despite completing successfully.
+To assist in diagnosing how active the scheme is, in both location and time, the 3D
+implicit and explicit components of the vertical velocity are available via XIOS as
+\texttt{wimp} and \texttt{wexp} respectively.  Likewise, the partitioning coefficient
+($C\kern-0.14em f$) is also available as \texttt{wi\_cff}. For a quick oversight of
+the schemes activity the global maximum values of the absolute implicit component
+of the vertical velocity and the partitioning coefficient are written to the netCDF
+version of the run statistics file (\texttt{run.stat.nc}) if this is active (see
+\autoref{sec:MISC_opt} for activation details).
+
+\autoref{fig:zad_Aimp_maxCf} shows examples of the maximum partitioning coefficient for
+the various overflow tests.  Note that the adaptive-implicit vertical advection scheme is
+active even in the base run with \forcode{nn_rdt=10.0s} adding to the evidence that the
+test case is close to stability limits even with this value. At the larger timesteps, the
+vertical velocity is treated mostly implicitly at some location throughout the run. The
+oscillatory nature of this measure appears to be linked to the progress of the plume front
+as each cusp is associated with the location of the maximum shifting to the adjacent cell.
+This is illustrated in \autoref{fig:zad_Aimp_maxCf_loc} where the i- and k- locations of the
+maximum have been overlaid for the base run case.
+
+\medskip
+\noindent Only limited tests have been performed in more realistic configurations. In the
+ORCA2\_ICE\_PISCES reference configuration the scheme does activate and passes
+restartability and reproducibility tests but it is unable to improve the model's stability
+enough to allow an increase in the model time-step. A view of the time-series of maximum
+partitioning coefficient (not shown here)  suggests that the default time-step of 5400s is
+already pushing at stability limits, especially in the initial start-up phase. The
+time-series does not, however, exhibit any of the 'cuspiness' found with the overflow
+tests.
+
+\medskip
+\noindent A short test with an eORCA1 configuration promises more since a test using a
+time-step of 3600s remains stable with \forcode{ln_zad_Aimp=.true.} whereas the
+time-step is limited to 2700s without.
 
 \begin{figure}[!t]
@@ -1440,4 +1473,26 @@
 \end{figure}
 
+\begin{figure}[!t]
+  \begin{center}
+    \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_maxCf}
+    \caption{
+      \protect\label{fig:zad_Aimp_maxCf}
+      The maximum partitioning coefficient during a series of test runs with increasing model timestep length.
+      At the larger timesteps, the vertical velocity is treated mostly implicitly at some location throughout 
+      the run.
+    }
+  \end{center}
+\end{figure}
+
+\begin{figure}[!t]
+  \begin{center}
+    \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_maxCf_loc}
+    \caption{
+      \protect\label{fig:zad_Aimp_maxCf_loc}
+      The maximum partitioning coefficient for the  \forcode{nn_rdt=10.0s} case overlaid with  information on the gridcell i- and k-
+      locations of the maximum value. 
+    }
+  \end{center}
+\end{figure}
 
 % ================================================================
Index: /NEMO/trunk/src/OCE/DYN/sshwzv.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/sshwzv.F90	(revision 11406)
+++ /NEMO/trunk/src/OCE/DYN/sshwzv.F90	(revision 11407)
@@ -286,5 +286,5 @@
       REAL(wp)             ::   zCu, zcff, z1_e3t                     ! 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_max = 0.30_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
@@ -346,5 +346,5 @@
                   wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk)
                   !
-                  Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient
+                  Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient below and in stp_ctl
                END DO
             END DO
@@ -353,5 +353,5 @@
       ELSE
          ! Fully explicit everywhere
-         Cu_adv(:,:,:) = 0._wp                          ! Reuse array to output coefficient
+         Cu_adv(:,:,:) = 0._wp                          ! Reuse array to output coefficient below and in stp_ctl
          wi    (:,:,:) = 0._wp
       ENDIF
Index: /NEMO/trunk/src/OCE/TRA/traadv_fct.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/traadv_fct.F90	(revision 11406)
+++ /NEMO/trunk/src/OCE/TRA/traadv_fct.F90	(revision 11407)
@@ -21,4 +21,5 @@
    USE diaar5         ! AR5 diagnostics
    USE phycst  , ONLY : rau0_rcp
+   USE zdf_oce , ONLY : ln_zad_Aimp
    !
    USE in_out_manager ! I/O manager
@@ -86,4 +87,6 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdx, ztrdy, ztrdz, zptry
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zwinf, zwdia, zwsup
+      LOGICAL  ::   ll_zAimp                                 ! flag to apply adaptive implicit vertical advection
       !!----------------------------------------------------------------------
       !
@@ -97,4 +100,5 @@
       l_hst = .FALSE.
       l_ptr = .FALSE.
+      ll_zAimp = .FALSE.
       IF( ( cdtype =='TRA' .AND. l_trdtra  ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) )       l_trd = .TRUE.
       IF(   cdtype =='TRA' .AND. ln_diaptr )                                                l_ptr = .TRUE. 
@@ -116,4 +120,22 @@
       !
       zwi(:,:,:) = 0._wp        
+      !
+      ! If adaptive vertical advection, check if it is needed on this PE at this time
+      IF( ln_zad_Aimp ) THEN
+         IF( MAXVAL( ABS( wi(:,:,:) ) ) > 0._wp ) ll_zAimp = .TRUE.
+      END IF
+      ! If active adaptive vertical advection, build tridiagonal matrix
+      IF( ll_zAimp ) THEN
+         ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk))
+         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.)
+                  zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) / e3t_a(ji,jj,jk)
+                  zwinf(ji,jj,jk) =  p2dt * MIN( wi(ji,jj,jk  ) , 0._wp ) / e3t_a(ji,jj,jk)
+                  zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t_a(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      END IF
       !
       DO jn = 1, kjpt            !==  loop over the tracers  ==!
@@ -169,4 +191,29 @@
             END DO
          END DO
+         
+         IF ( ll_zAimp ) THEN
+            CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 )
+            !
+            ztu(:,:,1) = 0._wp; ztu(:,:,jpk) = 0._wp
+            DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
+                     zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
+                     ztu(ji,jj,jk) =  0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.  
+                     pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztu(ji,jj,jk) - ztu(ji  ,jj  ,jk+1) ) &
+                        &                                  * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            zwz(:,:,:) = zwz(:,:,:) + ztu(:,:,:)
+            !
+         END IF
          !                
          IF( l_trd .OR. l_hst )  THEN             ! trend diagnostics (contribution of upstream fluxes)
@@ -280,4 +327,31 @@
          CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1., zwx, 'U', -1. , zwy, 'V', -1.,  zwz, 'W',  1. )
          !
+         !         
+         IF ( ll_zAimp ) THEN
+            DO jk = 1, jpkm1     !* trend and after field with monotonic scheme
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     !                             ! total intermediate advective trends
+                     ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                        &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
+                        &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
+                     ztu(ji,jj,jk)  = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            !
+            CALL tridia_solver( zwdia, zwsup, zwinf, ztu, ztu , 0 )
+            !
+            ztu(:,:,1) = 0._wp
+            DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask)
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
+                     zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
+                     zwz(ji,jj,jk) =  zwz(ji,jj,jk) + 0.5 * e1e2t(ji,jj) * ( zfp_wk * ztu(ji,jj,jk) + zfm_wk * ztu(ji,jj,jk-1) ) * wmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         END IF
          !        !==  monotonicity algorithm  ==!
          !
@@ -289,11 +363,33 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.  
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
-                     &                                   + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
-                     &                                   + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) &
-                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
+                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
+                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra / e3t_n(ji,jj,jk)
+                  zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !
+         IF ( ll_zAimp ) THEN
+            !
+            DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask)
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
+                     zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
+                     zwz(ji,jj,jk) = - 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.  
+                     pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) &
+                        &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         END IF         
          !
          IF( l_trd .OR. l_hst ) THEN   ! trend diagnostics // heat/salt transport
@@ -318,4 +414,7 @@
       END DO                     ! end of tracer loop
       !
+      IF ( ll_zAimp ) THEN
+         DEALLOCATE( zwdia, zwinf, zwsup )
+      ENDIF
       IF( l_trd .OR. l_hst ) THEN 
          DEALLOCATE( ztrdx, ztrdy, ztrdz )
Index: /NEMO/trunk/src/OCE/stpctl.F90
===================================================================
--- /NEMO/trunk/src/OCE/stpctl.F90	(revision 11406)
+++ /NEMO/trunk/src/OCE/stpctl.F90	(revision 11407)
@@ -96,5 +96,5 @@
             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  )
+               istatus = NF90_DEF_VAR( idrun,       'Cf_max', NF90_DOUBLE, (/ idtime /), idc1  )
             ENDIF
             istatus = NF90_ENDDEF(idrun)
@@ -123,5 +123,5 @@
       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
+         zmax(9) = MAXVAL(   Cu_adv(:,:,:)   , mask = tmask(:,:,:) == 1._wp ) ! partitioning coeff. max
       ENDIF
       !
