Index: branches/2015/dev_r5776_UKMO2_OBS_efficiency_improvs/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90
===================================================================
--- branches/2015/dev_r5776_UKMO2_OBS_efficiency_improvs/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90	(revision 5795)
+++ branches/2015/dev_r5776_UKMO2_OBS_efficiency_improvs/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90	(revision 6041)
@@ -8,5 +8,5 @@
    !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
    !!----------------------------------------------------------------------
-#if defined key_zdftmx   ||   defined key_esopa
+#if defined key_zdftmx
    !!----------------------------------------------------------------------
    !!   'key_zdftmx'                                  Tidal vertical mixing
@@ -54,5 +54,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -105,17 +105,17 @@
       !!              Koch-Larrouy et al. 2007, GRL.
       !!----------------------------------------------------------------------
-      USE oce, zav_tide  =>   ua    ! use ua as workspace
-      !!
       INTEGER, INTENT(in) ::   kt   ! ocean time-step 
-      !!
+      !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   ztpc         ! scalar workspace
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zkz
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zkz
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zav_tide
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('zdf_tmx')
       !
-      CALL wrk_alloc( jpi,jpj, zkz )
-
+      CALL wrk_alloc( jpi,jpj,       zkz )
+      CALL wrk_alloc( jpi,jpj,jpk,   zav_tide )
+      !
       !                          ! ----------------------- !
       !                          !  Standard tidal mixing  !  (compute zav_tide)
@@ -136,22 +136,19 @@
 
       DO jk = 2, jpkm1     !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s
-         DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
-            DO ji = 1, jpi
-               zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) * wmask(ji,jj,jk)  !kz max = 300 cm2/s
-            END DO
-         END DO
+         zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. ) * wmask(:,:,jk)  !kz max = 300 cm2/s
       END DO
 
       IF( kt == nit000 ) THEN       !* check at first time-step: diagnose the energy consumed by zav_tide
-         ztpc = 0.e0
+         ztpc = 0._wp
          DO jk= 1, jpk
             DO jj= 1, jpj
                DO ji= 1, jpi
-                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj)   &
-                     &         * MAX( 0.e0, rn2(ji,jj,jk) ) * zav_tide(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
+                  ztpc = ztpc + fse3w(ji,jj,jk) * e1e2t(ji,jj)                  &
+                     &        * MAX( 0.e0, rn2(ji,jj,jk) ) * zav_tide(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
                END DO
             END DO
          END DO
          ztpc= rau0 / ( rn_tfe * rn_me ) * ztpc
+         IF( lk_mpp )   CALL mpp_sum( ztpc )
          IF(lwp) WRITE(numout,*) 
          IF(lwp) WRITE(numout,*) '          N Total power consumption by av_tide    : ztpc = ', ztpc * 1.e-12 ,'TW'
@@ -167,13 +164,6 @@
       !                          ! ----------------------- !
       DO jk = 2, jpkm1              !* update momentum & tracer diffusivity with tidal mixing
-         DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
-            DO ji = 1, jpi
-               avt(ji,jj,jk) = avt(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk)
-               avm(ji,jj,jk) = avm(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      
-      DO jk = 2, jpkm1              !* update momentum & tracer diffusivity with tidal mixing
+         avt(:,:,jk) = avt(:,:,jk) + zav_tide(:,:,jk) * wmask(:,:,jk)
+         avm(:,:,jk) = avm(:,:,jk) + zav_tide(:,:,jk) * wmask(:,:,jk)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1  ! vector opt.
@@ -190,5 +180,6 @@
       IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_tide , clinfo1=' tmx - av_tide: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk)
       !
-      CALL wrk_dealloc( jpi,jpj, zkz )
+      CALL wrk_dealloc( jpi,jpj,       zkz )
+      CALL wrk_dealloc( jpi,jpj,jpk,   zav_tide )
       !
       IF( nn_timing == 1 )  CALL timing_stop('zdf_tmx')
@@ -239,5 +230,4 @@
       DO jk = 1, jpkm1             
          zdn2dz     (:,:,jk) = rn2(:,:,jk) - rn2(:,:,jk+1)           ! Vertical profile of dN2/dz
-!CDIR NOVERRCHK
          zempba_3d_1(:,:,jk) = SQRT(  MAX( 0.e0, rn2(:,:,jk) )  )    !    -        -    of N
          zempba_3d_2(:,:,jk) =        MAX( 0.e0, rn2(:,:,jk) )       !    -        -    of N^2
@@ -248,6 +238,6 @@
       zsum2(:,:) = 0.e0
       DO jk= 2, jpk
-         zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1)
-         zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1)               
+         zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk) * wmask(:,:,jk)
+         zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk) * wmask(:,:,jk)               
       END DO
       DO jj = 1, jpj
@@ -285,5 +275,5 @@
       zkz(:,:) = 0.e0               ! Associated potential energy consummed over the whole water column
       DO jk = 2, jpkm1
-         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zavt_itf(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1)
+         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zavt_itf(:,:,jk) * wmask(:,:,jk)
       END DO
 
@@ -295,5 +285,5 @@
 
       DO jk = 2, jpkm1              ! Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zavt_itf bound by 300 cm2/s
-         zavt_itf(:,:,jk) = zavt_itf(:,:,jk) * MIN( zkz(:,:), 120./10. ) * tmask(:,:,jk) * tmask(:,:,jk-1)   ! kz max = 120 cm2/s
+         zavt_itf(:,:,jk) = zavt_itf(:,:,jk) * MIN( zkz(:,:), 120./10. ) * wmask(:,:,jk)   ! kz max = 120 cm2/s
       END DO
 
@@ -303,9 +293,10 @@
             DO jj= 1, jpj
                DO ji= 1, jpi
-                  ztpc = ztpc + e1t(ji,jj) * e2t(ji,jj) * fse3w(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) )   &
-                     &                     * zavt_itf(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
+                  ztpc = ztpc + e1e2t(ji,jj) * fse3w(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) )   &
+                     &                       * zavt_itf(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
                END DO
             END DO
          END DO
+         IF( lk_mpp )   CALL mpp_sum( ztpc )
          ztpc= rau0 * ztpc / ( rn_me * rn_tfe_itf )
          IF(lwp) WRITE(numout,*) '          N Total power consumption by zavt_itf: ztpc = ', ztpc * 1.e-12 ,'TW'
@@ -361,15 +352,13 @@
       !!              Koch-Larrouy et al. 2007, GRL.
       !!----------------------------------------------------------------------
-      USE oce     ,         zav_tide =>  ua         ! ua used as workspace
-      !!
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   inum         ! local integer
       INTEGER  ::   ios
       REAL(wp) ::   ztpc, ze_z   ! local scalars
-      REAL(wp), DIMENSION(:,:)  , POINTER ::  zem2, zek1   ! read M2 and K1 tidal energy
-      REAL(wp), DIMENSION(:,:)  , POINTER ::  zkz          ! total M2, K1 and S2 tidal energy
-      REAL(wp), DIMENSION(:,:)  , POINTER ::  zfact        ! used for vertical structure function
-      REAL(wp), DIMENSION(:,:)  , POINTER ::  zhdep        ! Ocean depth 
-      REAL(wp), DIMENSION(:,:,:), POINTER ::  zpc      ! power consumption
+      REAL(wp), DIMENSION(:,:)  , POINTER ::  zem2, zek1     ! read M2 and K1 tidal energy
+      REAL(wp), DIMENSION(:,:)  , POINTER ::  zkz            ! total M2, K1 and S2 tidal energy
+      REAL(wp), DIMENSION(:,:)  , POINTER ::  zfact          ! used for vertical structure function
+      REAL(wp), DIMENSION(:,:)  , POINTER ::  zhdep          ! Ocean depth 
+      REAL(wp), DIMENSION(:,:,:), POINTER ::  zpc, zav_tide  ! power consumption
       !!
       NAMELIST/namzdf_tmx/ rn_htmx, rn_n2min, rn_tfe, rn_me, ln_tmx_itf, rn_tfe_itf
@@ -378,17 +367,17 @@
       IF( nn_timing == 1 )  CALL timing_start('zdf_tmx_init')
       !
-      CALL wrk_alloc( jpi,jpj, zem2, zek1, zkz, zfact, zhdep )
-      CALL wrk_alloc( jpi,jpj,jpk, zpc )
-      
-      REWIND( numnam_ref )              ! Namelist namzdf_tmx in reference namelist : Tidal Mixing
+      CALL wrk_alloc( jpi,jpj,       zem2, zek1, zkz, zfact, zhdep )
+      CALL wrk_alloc( jpi,jpj,jpk,   zpc, zav_tide )
+      !
+      REWIND( numnam_ref )             ! Namelist namzdf_tmx in reference namelist : Tidal Mixing
       READ  ( numnam_ref, namzdf_tmx, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in reference namelist', lwp )
-
-      REWIND( numnam_cfg )              ! Namelist namzdf_tmx in configuration namelist : Tidal Mixing
+      !
+      REWIND( numnam_cfg )             ! Namelist namzdf_tmx in configuration namelist : Tidal Mixing
       READ  ( numnam_cfg, namzdf_tmx, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in configuration namelist', lwp )
       IF(lwm) WRITE ( numond, namzdf_tmx )
-
-      IF(lwp) THEN                   ! Control print
+      !
+      IF(lwp) THEN                     ! Control print
          WRITE(numout,*)
          WRITE(numout,*) 'zdf_tmx_init : tidal mixing'
@@ -402,33 +391,31 @@
          WRITE(numout,*) '      ITF tidal dissipation efficiency      = ', rn_tfe_itf
       ENDIF
-
-      !                              ! allocate tmx arrays
+      !                                ! allocate tmx arrays
       IF( zdf_tmx_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tmx_init : unable to allocate tmx arrays' )
 
-      IF( ln_tmx_itf ) THEN          ! read the Indonesian Through Flow mask
+      IF( ln_tmx_itf ) THEN            ! read the Indonesian Through Flow mask
          CALL iom_open('mask_itf',inum)
          CALL iom_get (inum, jpdom_data, 'tmaskitf',mask_itf,1) ! 
          CALL iom_close(inum)
       ENDIF
-
-      ! read M2 tidal energy flux : W/m2  ( zem2 < 0 )
+      !                                ! read M2 tidal energy flux : W/m2  ( zem2 < 0 )
       CALL iom_open('M2rowdrg',inum)
       CALL iom_get (inum, jpdom_data, 'field',zem2,1) ! 
       CALL iom_close(inum)
-
-      ! read K1 tidal energy flux : W/m2  ( zek1 < 0 )
+      !                                ! read K1 tidal energy flux : W/m2  ( zek1 < 0 )
       CALL iom_open('K1rowdrg',inum)
       CALL iom_get (inum, jpdom_data, 'field',zek1,1) ! 
       CALL iom_close(inum)
- 
-      ! Total tidal energy ( M2, S2 and K1  with S2=(1/2)^2 * M2 )
-      ! only the energy available for mixing is taken into account,
-      ! (mixing efficiency tidal dissipation efficiency)
+      !                                ! Total tidal energy ( M2, S2 and K1  with S2=(1/2)^2 * M2 )
+      !                                ! only the energy available for mixing is taken into account,
+      !                                ! (mixing efficiency tidal dissipation efficiency)
       en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * ssmask(:,:)
 
 !============
 !TG: Bug for VVL? Should this section be moved out of _init and be updated at every timestep?
-      ! Vertical structure (az_tmx)
-      DO jj = 1, jpj                ! part independent of the level
+!!gm : you are right, but tidal mixing acts in deep ocean (H>500m) where e3 is O(100m)
+!!     the error is thus ~1% which I feel comfortable with, compared to uncertainties in tidal energy dissipation.
+      !                                ! Vertical structure (az_tmx)
+      DO jj = 1, jpj                         ! part independent of the level
          DO ji = 1, jpi
             zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean
@@ -437,5 +424,5 @@
          END DO
       END DO
-      DO jk= 1, jpk                 ! complete with the level-dependent part
+      DO jk= 1, jpk                          ! complete with the level-dependent part
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -445,5 +432,5 @@
       END DO
 !===========
-
+      !
       IF( nprint == 1 .AND. lwp ) THEN
          ! Control print
@@ -454,29 +441,25 @@
             zav_tide(:,:,jk) = az_tmx(:,:,jk) / MAX( rn_n2min, rn2(:,:,jk) )
          END DO
-
-         ztpc = 0.e0
+         !
+         ztpc = 0._wp
          zpc(:,:,:) = MAX(rn_n2min,rn2(:,:,:)) * zav_tide(:,:,:)
          DO jk= 2, jpkm1
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
+                  ztpc = ztpc + fse3w(ji,jj,jk) * e1e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
                END DO
             END DO
          END DO
+         IF( lk_mpp )   CALL mpp_sum( ztpc )
          ztpc= rau0 * 1/(rn_tfe * rn_me) * ztpc
-
+         !
          WRITE(numout,*) 
          WRITE(numout,*) '          Total power consumption of the tidally driven part of Kz : ztpc = ', ztpc * 1.e-12 ,'TW'
-
-
+         !
          ! control print 2
          zav_tide(:,:,:) = MIN( zav_tide(:,:,:), 60.e-4 )   
-         zkz(:,:) = 0.e0
+         zkz(:,:) = 0._wp
          DO jk = 2, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zkz(ji,jj) = zkz(ji,jj) + fse3w(ji,jj,jk) * MAX(0.e0, rn2(ji,jj,jk)) * rau0 * zav_tide(ji,jj,jk) * wmask(ji,jj,jk)
-               END DO
-            END DO
+               zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX(0.e0, rn2(:,:,jk)) * rau0 * zav_tide(:,:,jk) * wmask(:,:,jk)
          END DO
          ! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz
@@ -497,31 +480,28 @@
          END DO
          WRITE(numout,*) '          Min de zkz ', ztpc, ' Max = ', maxval(zkz(:,:) )
-
+         !
          DO jk = 2, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) * wmask(ji,jj,jk)  !kz max = 300 cm2/s
-               END DO
-            END DO
-         END DO
-         ztpc = 0.e0
+            zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. ) * wmask(:,:,jk)  !kz max = 300 cm2/s
+         END DO
+         ztpc = 0._wp
          zpc(:,:,:) = Max(0.e0,rn2(:,:,:)) * zav_tide(:,:,:)
          DO jk= 1, jpk
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
+                  ztpc = ztpc + fse3w(ji,jj,jk) * e1e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
                END DO
             END DO
          END DO
+         IF( lk_mpp )   CALL mpp_sum( ztpc )
          ztpc= rau0 * 1/(rn_tfe * rn_me) * ztpc
          WRITE(numout,*) '          2 Total power consumption of the tidally driven part of Kz : ztpc = ', ztpc * 1.e-12 ,'TW'
-
+!!gm bug mpp  in these diagnostics
          DO jk = 1, jpk
-            ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zav_tide(:,:,jk)     * tmask_i(:,:) )   &
-               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) )
-            ztpc = 1.E50
+            ze_z =                  SUM( e1e2t(:,:) * zav_tide(:,:,jk) * tmask_i(:,:) )   &
+               &     / MAX( 1.e-20, SUM( e1e2t(:,:) * wmask   (:,:,jk) * tmask_i(:,:) ) )
+            ztpc = 1.e50
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  IF( zav_tide(ji,jj,jk) /= 0.e0 )   ztpc =Min( ztpc, zav_tide(ji,jj,jk) )
+                  IF( zav_tide(ji,jj,jk) /= 0.e0 )   ztpc = MIN( ztpc, zav_tide(ji,jj,jk) )
                END DO
             END DO
@@ -530,5 +510,5 @@
          END DO
 
-         WRITE(numout,*) '          e_tide : ', SUM( e1t*e2t*en_tmx ) / ( rn_tfe * rn_me ) * 1.e-12, 'TW'
+         WRITE(numout,*) '          e_tide : ', SUM( e1e2t*en_tmx ) / ( rn_tfe * rn_me ) * 1.e-12, 'TW'
          WRITE(numout,*) 
          WRITE(numout,*) '          Initial profile of tidal vertical mixing'
@@ -539,16 +519,17 @@
                END DO
             END DO
-            ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zkz(:,:)     * tmask_i(:,:) )   &
-               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) )
+            ze_z =                  SUM( e1e2t(:,:) * zkz  (:,:)    * tmask_i(:,:) )   &
+               &     / MAX( 1.e-20, SUM( e1e2t(:,:) * wmask(:,:,jk) * tmask_i(:,:) ) )
             WRITE(numout,*) '                jk= ', jk,'   ', ze_z * 1.e4,' cm2/s'
          END DO
          DO jk = 1, jpk
             zkz(:,:) = az_tmx(:,:,jk) /rn_n2min
-            ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zkz(:,:)     * tmask_i(:,:) )   &
-               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) )
+            ze_z =                  SUM( e1e2t(:,:) * zkz  (:,:)    * tmask_i(:,:) )   &
+               &     / MAX( 1.e-20, SUM( e1e2t(:,:) * wmask(:,:,jk) * tmask_i(:,:) ) )
             WRITE(numout,*) 
             WRITE(numout,*) '          N2 min - jk= ', jk,'   ', ze_z * 1.e4,' cm2/s min= ',MINVAL(zkz)*1.e4,   &
                &       'max= ', MAXVAL(zkz)*1.e4, ' cm2/s'
          END DO
+!!gm  end bug mpp
          !
       ENDIF
