Index: /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Chapters/Chap_DIA.tex
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Chapters/Chap_DIA.tex	(revision 5959)
+++ /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Chapters/Chap_DIA.tex	(revision 5960)
@@ -1522,4 +1522,13 @@
 \\
 
+\section{Courant numbers}
+Courant numbers provide a theoretical indication of the model's numerical stability. The advective Courant numbers can be calculated according to
+\begin{equation}
+\label{eq:CFL}
+C_u = |u|\frac{\rdt}{e_{1u}}, \quad C_v = |v|\frac{\rdt}{e_{2v}}, \quad C_w = |w|\frac{\rdt}{e_{3w}}
+\end{equation}
+in the zonal, meridional and vertical directions respectively. The vertical component is included although it is not strictly valid as the vertical velocity is calculated from the continuity equation rather than as a prognostic variable. Physically this represents the rate at which information is propogated across a grid cell. Values greater than 1 indicate that information is propagated across more than one grid cell in a single time step.
+
+The variables can be activated by setting the \np{nn\_diacfl} namelist parameter to 1 in the \ngn{namctl} namelist. The diagnostics will be written out to an ascii file named cfl\_diagnostics.ascii. In this file the maximum value of $C_u$, $C_v$, and $C_w$ are printed at each timestep along with the coordinates of where the maximum value occurs. At the end of the model run the maximum value of $C_u$, $C_v$, and $C_w$ for the whole model run is printed along with the coordinates of each. The maximum values from the run are also copied to the ocean.output file. 
 
 
Index: /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Namelist/namctl
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Namelist/namctl	(revision 5959)
+++ /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Namelist/namctl	(revision 5960)
@@ -13,3 +13,4 @@
                            !     (no physical validity of the results)
    nn_timing   =    0      !  timing by routine activated (=1) creates timing.output file, or not (=0)
+   nn_diacfl   =    0      !  Write out cfl diagnostics (=1) in cfl_diagnostics.ascii, or not (=0)
 /
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5959)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5960)
@@ -1028,4 +1028,5 @@
                            !     (no physical validity of the results)
    nn_timing   =    0      !  timing by routine activated (=1) creates timing.output file, or not (=0)
+   nn_diacfl   =    0      !  Write out cfl diagnostics (=1) in cfl_diagnostics.ascii, or not (=0)
 /
 !-----------------------------------------------------------------------
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90	(revision 5959)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90	(revision 5960)
@@ -148,5 +148,5 @@
       NAMELIST/namctl/ ln_ctl  , nn_print, nn_ictls, nn_ictle,   &
          &             nn_isplt, nn_jsplt, nn_jctls, nn_jctle,   &
-         &             nn_bench, nn_timing
+         &             nn_bench, nn_timing, nn_diacfl
       NAMELIST/namcfg/ cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, &
          &             jpizoom, jpjzoom, jperio, ln_use_jattr
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OOO_SRC/nemogcm.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OOO_SRC/nemogcm.F90	(revision 5959)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OOO_SRC/nemogcm.F90	(revision 5960)
@@ -127,5 +127,5 @@
       NAMELIST/namctl/ ln_ctl, nn_print, nn_ictls, nn_ictle,   &
          &             nn_isplt, nn_jsplt, nn_jctls, nn_jctle,   &
-         &             nn_bench, nn_timing
+         &             nn_bench, nn_timing, nn_diacfl
       NAMELIST/namcfg/ cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, &
          &             jpizoom, jpjzoom, jperio, ln_use_jattr
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90	(revision 5960)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90	(revision 5960)
@@ -0,0 +1,207 @@
+MODULE diacfl
+   !!==============================================================================
+   !!                       ***  MODULE  diacfl  ***
+   !! Output CFL diagnostics to ascii file
+   !!==============================================================================
+   !! History :  1.0  !  2010-03  (E. Blockley)  Original code
+   !!                 !  2014-06  (T Graham) Removed CPP key & Updated to vn3.6
+   !! 
+   !!----------------------------------------------------------------------
+   !!   dia_cfl        : Compute and output Courant numbers at each timestep
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and active tracers
+   USE dom_oce         ! ocean space and time domain
+   USE lib_mpp         ! distribued memory computing
+   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
+   USE in_out_manager  ! I/O manager
+   USE domvvl     
+   USE timing          ! Performance output
+
+   IMPLICIT NONE
+   PRIVATE
+
+   REAL(wp) :: cu_max, cv_max, cw_max                      ! Run max U Courant number 
+   INTEGER, DIMENSION(3) :: cu_loc, cv_loc, cw_loc         ! Run max locations
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zcu_cfl           ! Courant number arrays
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zcv_cfl           ! Courant number arrays
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zcw_cfl           ! Courant number arrays
+
+   INTEGER  :: numcfl                                       ! outfile unit
+   CHARACTER(LEN=50) :: clname="cfl_diagnostics.ascii"      ! ascii filename
+
+   PUBLIC   dia_cfl       ! routine called by step.F90
+   PUBLIC   dia_cfl_init  ! routine called by nemogcm
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+
+CONTAINS
+
+
+   SUBROUTINE dia_cfl ( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_cfl  ***
+      !!
+      !! ** Purpose :  Compute the Courant numbers Cu=u*dt/dx and Cv=v*dt/dy
+      !!               and output to ascii file 'cfl_diagnostics.ascii'
+      !!----------------------------------------------------------------------
+
+      INTEGER, INTENT(in) ::  kt                            ! ocean time-step index
+
+      REAL(wp) :: zcu_max, zcv_max, zcw_max                 ! max Courant numbers per timestep
+      INTEGER, DIMENSION(3) :: zcu_loc, zcv_loc, zcw_loc    ! max Courant number locations
+
+      REAL(wp) :: dt                                        ! temporary scalars
+      INTEGER, DIMENSION(3) :: zlocu, zlocv, zlocw          ! temporary arrays 
+      INTEGER  :: ji, jj, jk                                ! dummy loop indices
+
+      
+      IF( nn_diacfl == 1) THEN
+         IF( nn_timing == 1 )   CALL timing_start('dia_cfl')
+         ! setup timestep multiplier to account for initial Eulerian timestep
+         IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    dt = rdt
+         ELSE                                        ;    dt = rdt * 2.0
+         ENDIF
+
+             ! calculate Courant numbers
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, fs_jpim1   ! vector opt.
+
+                  ! Courant number for x-direction (zonal current)
+                  zcu_cfl(ji,jj,jk) = ABS(un(ji,jj,jk))*dt/e1u(ji,jj)
+
+                  ! Courant number for y-direction (meridional current)
+                  zcv_cfl(ji,jj,jk) = ABS(vn(ji,jj,jk))*dt/e2v(ji,jj)
+
+                  ! Courant number for z-direction (vertical current)
+                  zcw_cfl(ji,jj,jk) = ABS(wn(ji,jj,jk))*dt/fse3w_n(ji,jj,jk)
+               END DO
+            END DO         
+         END DO
+
+         ! calculate maximum values and locations
+         IF( lk_mpp ) THEN
+            CALL mpp_maxloc(zcu_cfl,umask,zcu_max, zcu_loc(1), zcu_loc(2), zcu_loc(3))
+            CALL mpp_maxloc(zcv_cfl,vmask,zcv_max, zcv_loc(1), zcv_loc(2), zcv_loc(3))
+            CALL mpp_maxloc(zcw_cfl,tmask,zcw_max, zcw_loc(1), zcw_loc(2), zcw_loc(3))
+         ELSE
+            zlocu = MAXLOC( ABS( zcu_cfl(:,:,:) ) )
+            zcu_loc(1) = zlocu(1) + nimpp - 1
+            zcu_loc(2) = zlocu(2) + njmpp - 1
+            zcu_loc(3) = zlocu(3)
+            zcu_max = zcu_cfl(zcu_loc(1),zcu_loc(2),zcu_loc(3))
+
+            zlocv = MAXLOC( ABS( zcv_cfl(:,:,:) ) )
+            zcv_loc(1) = zlocv(1) + nimpp - 1
+            zcv_loc(2) = zlocv(2) + njmpp - 1
+            zcv_loc(3) = zlocv(3)
+            zcv_max = zcv_cfl(zcv_loc(1),zcv_loc(2),zcv_loc(3))
+
+            zlocw = MAXLOC( ABS( zcw_cfl(:,:,:) ) )
+            zcw_loc(1) = zlocw(1) + nimpp - 1
+            zcw_loc(2) = zlocw(2) + njmpp - 1
+            zcw_loc(3) = zlocw(3)
+            zcw_max = zcw_cfl(zcw_loc(1),zcw_loc(2),zcw_loc(3))
+         ENDIF
+      
+         ! write out to file
+         IF( lwp ) THEN
+            WRITE(numcfl,FMT='(2x,i4,5x,a6,5x,f6.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zcu_max, zcu_loc(1), zcu_loc(2), zcu_loc(3)
+            WRITE(numcfl,FMT='(11x,a6,5x,f6.4,1x,i4,1x,i4,1x,i4)') 'Max Cv', zcv_max, zcv_loc(1), zcv_loc(2), zcv_loc(3)
+            WRITE(numcfl,FMT='(11x,a6,5x,f6.4,1x,i4,1x,i4,1x,i4)') 'Max Cw', zcw_max, zcw_loc(1), zcw_loc(2), zcw_loc(3)
+         ENDIF
+
+         ! update maximum Courant numbers from whole run if applicable
+         IF( zcu_max > cu_max ) THEN
+            cu_max = zcu_max
+            cu_loc = zcu_loc
+         ENDIF
+         IF( zcv_max > cv_max ) THEN
+            cv_max = zcv_max
+            cv_loc = zcv_loc
+         ENDIF
+         IF( zcw_max > cw_max ) THEN
+            cw_max = zcw_max
+            cw_loc = zcw_loc
+         ENDIF
+
+         ! at end of run output max Cu and Cv and close ascii file
+         IF( kt == nitend .AND. lwp ) THEN
+            ! to ascii file
+            WRITE(numcfl,*) '******************************************'
+            WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cu', cu_max, cu_loc(1), cu_loc(2), cu_loc(3)
+            WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cu_max)
+            WRITE(numcfl,*) '******************************************'
+            WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cv', cv_max, cv_loc(1), cv_loc(2), cv_loc(3)
+            WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cv_max)
+            WRITE(numcfl,*) '******************************************'
+            WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cw', cw_max, cw_loc(1), cw_loc(2), cw_loc(3)
+            WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cw_max)
+            CLOSE( numcfl ) 
+
+            ! to ocean output
+            WRITE(numout,*)
+            WRITE(numout,*) 'dia_cfl     : Maximum Courant number information for the run:'
+            WRITE(numout,*) '~~~~~~~~~~~~'
+            WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cu', cu_max, 'at (i, j, k) = (', cu_loc(1), cu_loc(2), cu_loc(3), ')'
+            WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cu_max)
+            WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cv', cv_max, 'at (i, j, k) = (', cv_loc(1), cv_loc(2), cv_loc(3), ')'
+            WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cv_max)
+            WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cw', cw_max, 'at (i, j, k) = (', cw_loc(1), cw_loc(2), cw_loc(3), ')'
+            WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cw_max)
+
+         ENDIF
+
+         IF( nn_timing == 1 )   CALL timing_stop('dia_cfl')
+      ENDIF
+
+   END SUBROUTINE dia_cfl
+
+   SUBROUTINE dia_cfl_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_cfl_init  ***
+      !!                   
+      !! ** Purpose :   create output file, initialise arrays
+      !!----------------------------------------------------------------------
+
+
+      IF( nn_diacfl == 1 ) THEN
+         IF( nn_timing == 1 )   CALL timing_start('dia_cfl_init')
+
+         cu_max=0.0
+         cv_max=0.0
+         cw_max=0.0
+
+         ALLOCATE( zcu_cfl(jpi, jpj, jpk), zcv_cfl(jpi, jpj, jpk), zcw_cfl(jpi, jpj, jpk) )
+
+         zcu_cfl(:,:,:)=0.0
+         zcv_cfl(:,:,:)=0.0
+         zcw_cfl(:,:,:)=0.0
+
+         IF( lwp ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'dia_cfl     : Outputting CFL diagnostics to '//TRIM(clname)
+            WRITE(numout,*) '~~~~~~~~~~~~'
+            WRITE(numout,*)
+
+            ! create output ascii file
+            CALL ctl_opn( numcfl, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, 1 )
+            WRITE(numcfl,*) 'Timestep  Direction  Max C     i    j    k'
+            WRITE(numcfl,*) '******************************************'
+         ENDIF
+
+         IF( nn_timing == 1 )   CALL timing_stop('dia_cfl_init')
+
+      ENDIF
+
+   END SUBROUTINE dia_cfl_init
+
+END MODULE diacfl
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 5959)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 5960)
@@ -99,4 +99,5 @@
    LOGICAL ::   ln_ctl       !: run control for debugging
    INTEGER ::   nn_timing    !: run control for timing
+   INTEGER ::   nn_diacfl    !: flag whether to create CFL diagnostics
    INTEGER ::   nn_print     !: level of print (0 no print)
    INTEGER ::   nn_ictls     !: Start i indice for the SUM control
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 5959)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 5960)
@@ -66,4 +66,5 @@
    USE diadct          ! sections transports           (dia_dct_init routine)
    USE diaobs          ! Observation diagnostics       (dia_obs_init routine)
+   USE diacfl          ! CFL diagnostics               (dia_cfl_init routine)
    USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
    USE step            ! NEMO time-stepping                 (stp     routine)
@@ -231,5 +232,5 @@
       NAMELIST/namctl/ ln_ctl  , nn_print, nn_ictls, nn_ictle,   &
          &             nn_isplt, nn_jsplt, nn_jctls, nn_jctle,   &
-         &             nn_bench, nn_timing
+         &             nn_bench, nn_timing, nn_diacfl
       NAMELIST/namcfg/ cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, &
          &             jpizoom, jpjzoom, jperio, ln_use_jattr
@@ -458,4 +459,5 @@
       !                                      ! Diagnostics
       IF( lk_floats     )   CALL     flo_init   ! drifting Floats
+                            CALL dia_cfl_init   ! Initialise CFL diagnostics
       IF( lk_diaar5     )   CALL dia_ar5_init   ! ar5 diag
                             CALL dia_ptr_init   ! Poleward TRansports initialization
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5959)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5960)
@@ -234,11 +234,12 @@
       ! diagnostics and outputs             
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( lk_floats  )   CALL flo_stp       ( kstp )  ! drifting Floats
-      IF( lk_diahth  )   CALL dia_hth       ( kstp )  ! Thermocline depth (20 degres isotherm depth)
-      IF(.NOT.ln_cpl )   CALL dia_fwb       ( kstp )  ! Fresh water budget diagnostics
-      IF( lk_diadct  )   CALL dia_dct       ( kstp )  ! Transports
-      IF( lk_diaar5  )   CALL dia_ar5       ( kstp )  ! ar5 diag
-      IF( lk_diaharm )   CALL dia_harm      ( kstp )  ! Tidal harmonic analysis
-                         CALL dia_wri       ( kstp )  ! ocean model: outputs
+      IF( lk_floats  )   CALL flo_stp( kstp )         ! drifting Floats
+      IF( nn_diacfl == 1 )   CALL dia_cfl( kstp )         ! Courant number diagnostics
+      IF( lk_diahth  )   CALL dia_hth( kstp )         ! Thermocline depth (20 degres isotherm depth)
+      IF(.NOT.ln_cpl )   CALL dia_fwb( kstp )         ! Fresh water budget diagnostics
+      IF( lk_diadct  )   CALL dia_dct( kstp )         ! Transports
+      IF( lk_diaar5  )   CALL dia_ar5( kstp )         ! ar5 diag
+      IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis
+                         CALL dia_wri( kstp )         ! ocean model: outputs
       !
       IF( ln_crs     )   CALL crs_fld       ( kstp )  ! ocean model: online field coarsening & output
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 5959)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 5960)
@@ -82,4 +82,5 @@
    USE diahsb           ! heat, salt and volume budgets    (dia_hsb routine)
    USE diaharm
+   USE diacfl
    USE flo_oce          ! floats variables
    USE floats           ! floats computation               (flo_stp routine)
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAS_SRC/nemogcm.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAS_SRC/nemogcm.F90	(revision 5959)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAS_SRC/nemogcm.F90	(revision 5960)
@@ -184,5 +184,5 @@
       NAMELIST/namctl/ ln_ctl  , nn_print, nn_ictls, nn_ictle,   &
          &             nn_isplt, nn_jsplt, nn_jctls, nn_jctle,   &
-         &             nn_bench, nn_timing
+         &             nn_bench, nn_timing, nn_diacfl
       NAMELIST/namcfg/ cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, &
          &             jpizoom, jpjzoom, jperio, ln_use_jattr
