Index: trunk/NEMO/OPA_SRC/ZDF/zdftmx.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdftmx.F90	(revision 1517)
+++ trunk/NEMO/OPA_SRC/ZDF/zdftmx.F90	(revision 1518)
@@ -27,19 +27,19 @@
    PRIVATE
 
-   PUBLIC zdf_tmx          ! called in step module 
+   PUBLIC   zdf_tmx    ! called in step module 
 
    LOGICAL, PUBLIC, PARAMETER ::   lk_zdftmx = .TRUE.    !: tidal mixing flag
 
-   !                                                    !!* Namelist  namtmx : tidal mixing *
-   REAL(wp) ::  rn_htmx  = 500.                          ! vertical decay scale for turbulence (meters)
-   REAL(wp) ::  rn_n2min = 1.e-8                         ! threshold of the Brunt-Vaisala frequency (s-1)
-   REAL(wp) ::  rn_tfe   = 1./3.                         ! tidal dissipation efficiency (St Laurent et al. 2002)
-   REAL(wp) ::  rn_tfe_itf  = 1.                         ! ITF tidal dissipation efficiency (St Laurent et al. 2002)
-   REAL(wp) ::  rn_me    = 0.2                           ! mixing efficiency (Osborn 1980)
-
-   REAL(wp), DIMENSION(jpi,jpj) ::   en_tmx              ! energy available for tidal mixing (W/m2)
-   REAL(wp), DIMENSION(jpi,jpj) ::   mask_itf            ! mask to use over Indonesian area
-
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   az_tmx          ! coefficient used to evaluate the tidal induced Kz
+   !                                  !!* Namelist  namtmx : tidal mixing *
+   REAL(wp) ::  rn_htmx    = 500.      ! vertical decay scale for turbulence (meters)
+   REAL(wp) ::  rn_n2min   = 1.e-8     ! threshold of the Brunt-Vaisala frequency (s-1)
+   REAL(wp) ::  rn_tfe     = 1./3.     ! tidal dissipation efficiency (St Laurent et al. 2002)
+   REAL(wp) ::  rn_me      = 0.2       ! mixing efficiency (Osborn 1980)
+   LOGICAL  ::  ln_tmx_itf = .TRUE.    ! Indonesian Through Flow (ITF): Koch-Larrouy et al. (2007) parameterization
+   REAL(wp) ::  rn_tfe_itf = 1.        ! ITF tidal dissipation efficiency (St Laurent et al. 2002)
+
+   REAL(wp), DIMENSION(jpi,jpj)     ::   en_tmx     ! energy available for tidal mixing (W/m2)
+   REAL(wp), DIMENSION(jpi,jpj)     ::   mask_itf   ! mask to use over Indonesian area
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   az_tmx     ! coefficient used to evaluate the tidal induced Kz
 
    !! * Substitutions
@@ -97,7 +97,7 @@
       IF( kt == nit000  )   CALL zdf_tmx_init      ! Initialization (first time-step only)
 
-      !                     ! ----------------------- !
-      !                     !  Standard tidal mixing  !  (compute av_tide)
-      !                     ! ----------------------- !
+      !                                        ! ----------------------- !
+      !                                        !  Standard tidal mixing  !  (compute av_tide)
+      !                                        ! ----------------------- !
       !                             !* First estimation (with n2 bound by rn_n2min) bounded by 60 cm2/s
       av_tide(:,:,:) = MIN(  60.e-4, az_tmx(:,:,:) / MAX( rn_n2min, rn2(:,:,:) )  )
@@ -133,11 +133,11 @@
       ENDIF
        
-      !                     ! ----------------------- !
-      CALL tmx_itf( kt )    !    ITF  tidal mixing    !  (update av_tide)
-      !                     ! ----------------------- !
-
-      !                     ! ----------------------- !
-      !                     !   Update  mixing coefs  !                          
-      !                     ! ----------------------- !
+      !                                        ! ----------------------- !
+      IF( ln_tmx_itf )   CALL tmx_itf( kt )    !    ITF  tidal mixing    !  (update av_tide)
+      !                                        ! ----------------------- !
+
+      !                                        ! ----------------------- !
+      !                                        !   Update  mixing coefs  !                          
+      !                                        ! ----------------------- !
       DO jk = 2, jpkm1              !* update momentum & tracer diffusivity with tidal mixing
          avt(:,:,jk) = avt(:,:,jk) + av_tide(:,:,jk)
@@ -180,115 +180,50 @@
       INTEGER  ::   ji, jj, jk    ! dummy loop indices
       REAL(wp) ::   zcoef, ztpc   ! temporary scalar
-      REAL(wp), DIMENSION(jpi,jpj) ::   zkz   ! temporary 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj) ::   zsum1 , zsum2 , zsum    
-      REAL(wp), DIMENSION(jpi,jpj) ::   z1sum1, z1sum2, z1sum  
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zempba_3d_1, zempba_3d_2   ! 3D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zempba_3d  , zdn2dz        !  -      -
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zavt_itf                   !  -      -
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zkz                        ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zsum1 , zsum2 , zsum       !  -      -
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zempba_3d_1, zempba_3d_2   ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zempba_3d  , zdn2dz        !  -      -
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zavt_itf                   !  -      -
       !!----------------------------------------------------------------------
 
       !                             ! compute the form function using N2 at each time step
-!!gm better/faster coding  (not tmask as rn2 is already masked)
-!      zempba_3d_1(:,:,jpk) = 0.e0
-!      zempba_3d_2(:,:,jpk) = 0.e0
-!      DO jk = 1, jpkm1             
-!         zdn2dz     (:,:,jk) = rn2(:,:,jk) - rn2(:,:,jk+1)           ! Vertical profile of dN2/dz
+      zempba_3d_1(:,:,jpk) = 0.e0
+      zempba_3d_2(:,:,jpk) = 0.e0
+      DO jk = 1, jpkm1             
+         zdn2dz     (:,:,jk) = rn2(:,:,jk) - rn2(:,:,jk+1)           ! Vertical profile of dN2/dz
 !CDIR NOVERRCHK
-!         zempba_3d_1(:,:,jk) = SQRT(  MAX( O.e0, rn2(:,:,jk) )  )    !    -        -    of N
-!         zempba_3d_2(:,:,jk) =        MAX( O.e0, rn2(:,:,jk) )       !    -        -    of N^2
-!      END DO
-!!gm replacing:
-      zempba_3d_1(:,:,:) = 0.e0
-      zempba_3d_2(:,:,:) = 0.e0
-   
-      DO jk = 1, jpk-1              ! Vertical profile dN2/dz
-         DO jj = 1, jpj               
-            DO ji = 1, jpi
-               zdn2dz(ji,jj,jk) =(rn2(ji,jj,jk) - rn2(ji,jj,jk+1) )
-               IF( rn2(ji,jj,jk) > 0 ) THEN 
-                  zempba_3d_1(ji,jj,jk)=SQRT(rn2(ji,jj,jk))*tmask(ji,jj,jk)
-                  zempba_3d_2(ji,jj,jk)=rn2(ji,jj,jk)*tmask(ji,jj,jk)
-               ENDIF
-            END DO
-         END DO
-      END DO
-!!gm end
-
-!!gm better/faster coding  (suppression of z1sum, z1sum1 and z1sum2)
-!      zsum (:,:) = 0.e0
-!      zsum1(:,:) = 0.e0
-!      zsum2(:,:) = 0.e0
-!      DO jk= 2, jpk
-!         zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk)
-!         zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk)                
-!      END DO
-!      DO jj = 1, jpj
-!         DO ji = 1, jpi
-!            IF( zsum1(ji,jj) /= 0.e0 )   zsum1(ji,jj) = 1.e0 / zsum1(ji,jj)
-!            IF( zsum2(ji,jj) /= 0.e0 )   zsum2(ji,jj) = 1.e0 / zsum2(ji,jj)                
-!         END DO
-!      END DO
-!
-!      DO jk= 1, jpk
-!         DO jj = 1, jpj
-!            DO ji = 1, jpi
-!               zcoef = 0.5 - SIGN( 0.5, zdn2dz(ji,jj,jk) )       ! =0 if dN2/dz > 0, =1 otherwise 
-!               ztpc  = zempba_3d_1(ji,jj,jk) * z1sum1(ji,jj) *        zcoef     &
-!                  &  + zempba_3d_2(ji,jj,jk) * z1sum2(ji,jj) * ( 1. - zcoef )
-!               !
-!               zempba_3d(ji,jj,jk) =               ztpc 
-!               zsum     (ji,jj)    = zsum(ji,jj) + ztpc * fse3w(ji,jj,jk)
-!            END DO
-!         END DO
-!       END DO
-!       DO jj = 1, jpj
-!          DO ji = 1, jpi
-!             IF( zsum(ji,jj) > 0.e0 )   zsum(ji,jj) = en_tmx(:,:) / zsum(ji,jj)                
-!          END DO
-!       END DO
-!
-!      !                             ! first estimation bounded by 10 cm2/s (with n2 bounded by rn_n2min) 
-!      zcoef = rn_tfe_itf / ( rn_tfe * rau0 )
-!      DO jk = 1, jpk
-!         zavt_itf(:,:,jk) = MIN(  10.e-4, zcoef * zsum(:,:) * zempba_3d(:,:,jk)   &
-!            &                                      / MAX( rn_n2min, rn2(:,:,jk) ) * tmask(:,:,jk)  )
-!      END DO           
-!!gm replacing:
-
-      zsum1(:,:)=0.e0
-      zsum2(:,:)=0.e0
+         zempba_3d_1(:,:,jk) = SQRT(  MAX( 0.e0, rn2(:,:,jk) )  )    !    -        -    of N
+         zempba_3d_2(:,:,jk) =        MAX( 0.e0, rn2(:,:,jk) )       !    -        -    of N^2
+      END DO
+      !
+      zsum (:,:) = 0.e0
+      zsum1(:,:) = 0.e0
+      zsum2(:,:) = 0.e0
       DO jk= 2, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zsum1(ji,jj)=zsum1(ji,jj)+zempba_3d_1(ji,jj,jk)*fse3w(ji,jj,jk)
-               zsum2(ji,jj)=zsum2(ji,jj)+zempba_3d_2(ji,jj,jk)*fse3w(ji,jj,jk)                
-            END DO
-         END DO
-      END DO
-
-      z1sum1(:,:)=0.e0
-      z1sum2(:,:)=0.e0
+         zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk)
+         zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk)                
+      END DO
       DO jj = 1, jpj
-        DO ji = 1, jpi
-         IF (zsum1(ji,jj) /= 0) z1sum1(ji,jj)=1.0/zsum1(ji,jj)
-         IF (zsum2(ji,jj) /= 0) z1sum2(ji,jj)=1.0/zsum2(ji,jj)                
-        END DO
-      END DO
-
-      zsum (:,:)=0.e0
-      z1sum(:,:)=0.e0
+         DO ji = 1, jpi
+            IF( zsum1(ji,jj) /= 0.e0 )   zsum1(ji,jj) = 1.e0 / zsum1(ji,jj)
+            IF( zsum2(ji,jj) /= 0.e0 )   zsum2(ji,jj) = 1.e0 / zsum2(ji,jj)                
+         END DO
+      END DO
+
       DO jk= 1, jpk
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zempba_3d(ji,jj,jk)=zempba_3d_1(ji,jj,jk)*z1sum1(ji,jj)*0.5*(1-SIGN(1.0,zdn2dz(ji,jj,jk)))       &
-                  &               +zempba_3d_2(ji,jj,jk)*z1sum2(ji,jj)*0.5*(1+SIGN(1.0,zdn2dz(ji,jj,jk)))
-               zsum(ji,jj)=zsum(ji,jj)+zempba_3d(ji,jj,jk)*fse3w(ji,jj,jk)
+               zcoef = 0.5 - SIGN( 0.5, zdn2dz(ji,jj,jk) )       ! =0 if dN2/dz > 0, =1 otherwise 
+               ztpc  = zempba_3d_1(ji,jj,jk) * zsum1(ji,jj) *        zcoef     &
+                  &  + zempba_3d_2(ji,jj,jk) * zsum2(ji,jj) * ( 1. - zcoef )
+               !
+               zempba_3d(ji,jj,jk) =               ztpc 
+               zsum     (ji,jj)    = zsum(ji,jj) + ztpc * fse3w(ji,jj,jk)
             END DO
          END DO
        END DO
-
        DO jj = 1, jpj
           DO ji = 1, jpi
-             IF( zsum(ji,jj) > 0 ) z1sum(ji,jj)=1.0/zsum(ji,jj)                
+             IF( zsum(ji,jj) > 0.e0 )   zsum(ji,jj) = 1.e0 / zsum(ji,jj)                
           END DO
        END DO
@@ -296,9 +231,8 @@
       !                             ! first estimation bounded by 10 cm2/s (with n2 bounded by rn_n2min) 
       zcoef = rn_tfe_itf / ( rn_tfe * rau0 )
-      DO jk= 1, jpk
-         zavt_itf(:,:,jk) = MIN(  10.e-4, zcoef * en_tmx(:,:) * z1sum(:,:) * zempba_3d(:,:,jk)   &
-            &                                   / MAX( rn_n2min, rn2(:,:,jk) ) * tmask(:,:,jk)  )
+      DO jk = 1, jpk
+         zavt_itf(:,:,jk) = MIN(  10.e-4, zcoef * en_tmx(:,:) * zsum(:,:) * zempba_3d(:,:,jk)   &
+            &                                      / MAX( rn_n2min, rn2(:,:,jk) ) * tmask(:,:,jk)  )
       END DO           
-!!gm end
 
       zkz(:,:) = 0.e0               ! Associated potential energy consummed over the whole water column
@@ -375,14 +309,14 @@
       !!              Koch-Larrouy et al. 2007, GRL.
       !!----------------------------------------------------------------------
-      INTEGER ::   ji, jj, jk                           ! dummy loop indices
-      INTEGER ::   inum                                 ! temporary logical unit
-      REAL(wp) ::   ztpc, ze_z                          ! total power consumption
-      REAL(wp), DIMENSION(jpi,jpj) ::  zem2, zek1       ! read M2 and K1 tidal energy
-      REAL(wp), DIMENSION(jpi,jpj) ::  zkz              ! total M2, K1 and S2 tidal energy
-      REAL(wp), DIMENSION(jpi,jpj) ::  zfact            ! used for vertical structure function
-      REAL(wp), DIMENSION(jpi,jpj) ::  zhdep                ! Ocean depth 
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zpc          ! power consumption
-      !!
-      NAMELIST/namtmx/ rn_htmx, rn_n2min, rn_tfe, rn_tfe_itf, rn_me
+      INTEGER ::   ji, jj, jk    ! dummy loop indices
+      INTEGER ::   inum          ! temporary logical unit
+      REAL(wp) ::   ztpc, ze_z   ! total power consumption
+      REAL(wp), DIMENSION(jpi,jpj) ::  zem2, zek1   ! read M2 and K1 tidal energy
+      REAL(wp), DIMENSION(jpi,jpj) ::  zkz          ! total M2, K1 and S2 tidal energy
+      REAL(wp), DIMENSION(jpi,jpj) ::  zfact        ! used for vertical structure function
+      REAL(wp), DIMENSION(jpi,jpj) ::  zhdep        ! Ocean depth 
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zpc      ! power consumption
+      !!
+      NAMELIST/namtmx/ rn_htmx, rn_n2min, rn_tfe, rn_me, ln_tmx_itf, rn_tfe_itf
       !!----------------------------------------------------------------------
 
@@ -404,13 +338,15 @@
          WRITE(numout,*) '             Brunt-Vaisala frequency threshold     = ', rn_n2min
          WRITE(numout,*) '             Tidal dissipation efficiency          = ', rn_tfe
+         WRITE(numout,*) '             Mixing efficiency                     = ', rn_me
+         WRITE(numout,*) '             ITF specific parameterisation         = ', ln_tmx_itf
          WRITE(numout,*) '             ITF tidal dissipation efficiency      = ', rn_tfe_itf
-         WRITE(numout,*) '             Mixing efficiency                     = ', rn_me
          WRITE(numout,*)
       ENDIF
 
-      ! read mask kz banda
-      CALL iom_open('mask_itf',inum)
-      CALL iom_get (inum, jpdom_data, 'tmaskitf',mask_itf,1) ! 
-      CALL iom_close(inum)
+      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 )
