Index: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90
===================================================================
--- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90	(revision 8185)
+++ branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90	(revision 8568)
@@ -1,10 +1,10 @@
 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
-   !! 
+   !!======================================================================
+   !! History :  3.4  !  2010-03  (E. Blockley)  Original code
+   !!            3.6  !  2014-06  (T. Graham)  Removed CPP key & Updated to vn3.6
+   !!            4.0  !  2017-09  (G. Madec)  style + comments
    !!----------------------------------------------------------------------
    !!   dia_cfl        : Compute and output Courant numbers at each timestep
@@ -12,8 +12,9 @@
    USE oce             ! ocean dynamics and active tracers
    USE dom_oce         ! ocean space and time domain
+   USE domvvl          ! 
+   !
    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
 
@@ -21,12 +22,14 @@
    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
+   CHARACTER(LEN=50) :: clname="cfl_diagnostics.ascii"    ! ascii filename
+   INTEGER           :: numcfl                            ! outfile unit
+   !
+   INTEGER, DIMENSION(3) ::   nCu_loc, nCv_loc, nCw_loc   ! U, V, and W run max locations in the global domain
+   REAL(wp)              ::   rCu_max, rCv_max, rCw_max   ! associated run max Courant number 
 
-   INTEGER  :: numcfl                                       ! outfile unit
-   CHARACTER(LEN=50) :: clname="cfl_diagnostics.ascii"      ! ascii filename
+!!gm CAUTION: need to declare these arrays here, otherwise the calculation fails in multi-proc !
+!!gm          8 don't understand why.
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace
+!!gm end
 
    PUBLIC   dia_cfl       ! routine called by step.F90
@@ -40,8 +43,5 @@
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
-
 CONTAINS
-
 
    SUBROUTINE dia_cfl ( kt )
@@ -52,116 +52,94 @@
       !!               and output to ascii file 'cfl_diagnostics.ascii'
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !
+      INTEGER :: ji, jj, jk   ! dummy loop indices
+      REAL(wp)::   z2dt, zCu_max, zCv_max, zCw_max       ! local scalars
+      INTEGER , DIMENSION(3)           ::   iloc_u , iloc_v , iloc_w , iloc   ! workspace
+!!gm this does not work      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace
+      !!----------------------------------------------------------------------
+      !
+      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   ;    z2dt = rdt
+      ELSE                                        ;    z2dt = rdt * 2._wp
+      ENDIF
+      !
+      !                
+      DO jk = 1, jpk       ! calculate Courant numbers
+         DO jj = 1, jpj
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * z2dt / e1u  (ji,jj)      ! for i-direction
+               zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * z2dt / e2v  (ji,jj)      ! for j-direction
+               zCw_cfl(ji,jj,jk) = ABS( wn(ji,jj,jk) ) * z2dt / e3w_n(ji,jj,jk)   ! for k-direction
+            END DO
+         END DO         
+      END DO
+      !
+      !                    ! calculate maximum values and locations
+      IF( lk_mpp ) THEN
+         CALL mpp_maxloc( zCu_cfl, umask, zCu_max, iloc_u(1), iloc_u(2), iloc_u(3) )
+         CALL mpp_maxloc( zCv_cfl, vmask, zCv_max, iloc_v(1), iloc_v(2), iloc_v(3) )
+         CALL mpp_maxloc( zCw_cfl, wmask, zCw_max, iloc_w(1), iloc_w(2), iloc_w(3) )
+      ELSE
+         iloc = MAXLOC( ABS( zcu_cfl(:,:,:) ) )
+         iloc_u(1) = iloc(1) + nimpp - 1
+         iloc_u(2) = iloc(2) + njmpp - 1
+         iloc_u(3) = iloc(3)
+         zCu_max = zCu_cfl(iloc(1),iloc(2),iloc(3))
+         !
+         iloc = MAXLOC( ABS( zcv_cfl(:,:,:) ) )
+         iloc_v(1) = iloc(1) + nimpp - 1
+         iloc_v(2) = iloc(2) + njmpp - 1
+         iloc_v(3) = iloc(3)
+         zCv_max = zCv_cfl(iloc(1),iloc(2),iloc(3))
+         !
+         iloc = MAXLOC( ABS( zcw_cfl(:,:,:) ) )
+         iloc_w(1) = iloc(1) + nimpp - 1
+         iloc_w(2) = iloc(2) + njmpp - 1
+         iloc_w(3) = iloc(3)
+         zCw_max = zCw_cfl(iloc(1),iloc(2),iloc(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, iloc_u(1), iloc_u(2), iloc_u(3)
+         WRITE(numcfl,FMT='(11x,     a6,5x,f6.4,1x,i4,1x,i4,1x,i4)')     'Max Cv', zCv_max, iloc_v(1), iloc_v(2), iloc_v(3)
+         WRITE(numcfl,FMT='(11x,     a6,5x,f6.4,1x,i4,1x,i4,1x,i4)')     'Max Cw', zCw_max, iloc_w(1), iloc_w(2), iloc_w(3)
+      ENDIF
+      !
+      !                    ! update maximum Courant numbers from whole run if applicable
+      IF( zCu_max > rCu_max ) THEN   ;   rCu_max = zCu_max   ;   nCu_loc(:) = iloc_u(:)   ;   ENDIF
+      IF( zCv_max > rCv_max ) THEN   ;   rCv_max = zCv_max   ;   nCv_loc(:) = iloc_v(:)   ;   ENDIF
+      IF( zCw_max > rCw_max ) THEN   ;   rCw_max = zCw_max   ;   nCw_loc(:) = iloc_w(:)   ;   ENDIF
 
-      INTEGER, INTENT(in) ::  kt                            ! ocean time-step index
+      !                    ! 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', rCu_max, nCu_loc(1), nCu_loc(2), nCu_loc(3)
+         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCu_max
+         WRITE(numcfl,*) '******************************************'
+         WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cv', rCv_max, nCv_loc(1), nCv_loc(2), nCv_loc(3)
+         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCv_max
+         WRITE(numcfl,*) '******************************************'
+         WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cw', rCw_max, nCw_loc(1), nCw_loc(2), nCw_loc(3)
+         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCw_max
+         CLOSE( numcfl ) 
+         !
+         ! to ocean output
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_cfl : Maximum Courant number information for the run '
+         WRITE(numout,*) '~~~~~~~'
+         WRITE(numout,*) '   Max Cu = ', rCu_max, ' at (i,j,k) = (',nCu_loc(1),nCu_loc(2),nCu_loc(3),') => dt/C = ', z2dt/rCu_max
+         WRITE(numout,*) '   Max Cv = ', rCv_max, ' at (i,j,k) = (',nCv_loc(1),nCv_loc(2),nCv_loc(3),') => dt/C = ', z2dt/rCv_max
+         WRITE(numout,*) '   Max Cw = ', rCw_max, ' at (i,j,k) = (',nCw_loc(1),nCw_loc(2),nCw_loc(3),') => dt/C = ', z2dt/rCw_max
+      ENDIF
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('dia_cfl')
+      !
+   END SUBROUTINE dia_cfl
 
-      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/e3w_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
@@ -171,36 +149,27 @@
       !! ** 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')
-
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_cfl : Outputting CFL diagnostics to ',TRIM(clname), ' file'
+         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
-
+      !
+      rCu_max = 0._wp
+      rCv_max = 0._wp
+      rCw_max = 0._wp
+      !
+!!gm required to work
+      ALLOCATE ( zCu_cfl(jpi,jpj,jpk), zCv_cfl(jpi,jpj,jpk), zCw_cfl(jpi,jpj,jpk) )
+!!gm end
+      !      
    END SUBROUTINE dia_cfl_init
 
+   !!======================================================================
 END MODULE diacfl
