Index: trunk/NEMO/OPA_SRC/ZDF/zdftmx.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdftmx.F90	(revision 1495)
+++ trunk/NEMO/OPA_SRC/ZDF/zdftmx.F90	(revision 1496)
@@ -17,9 +17,10 @@
    USE dom_oce         ! ocean space and time domain variables
    USE zdf_oce         ! ocean vertical physics variables
-   USE in_out_manager  ! I/O manager
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE eosbn2          ! ocean equation of state
    USE phycst          ! physical constants
    USE prtctl          ! Print control
+   USE in_out_manager  ! I/O manager
+   USE iom             ! I/O Manager
 
    IMPLICIT NONE
@@ -30,5 +31,5 @@
    LOGICAL, PUBLIC, PARAMETER ::   lk_zdftmx = .TRUE.    !: tidal mixing flag
 
-   !                                                    !!* Namelist  namtmx *
+   !                                                    !!* 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)
@@ -58,18 +59,31 @@
       !!                   
       !! ** Purpose :   add to the vertical mixing coefficients the effect of
-      !!      tidal mixing (Simmons et al 2004).
+      !!              tidal mixing (Simmons et al 2004).
       !!
       !! ** Method  : - tidal-induced vertical mixing is given by:
-      !!      Kz_tides = az_tmx / rn2, where az_tmx is a 2D coefficient set in
-      !!      zdf_tmx_init. This mixing is added to avt, avmu and avmv as
-      !!      follows:  avt  = avt  +    az_tmx
-      !!                avmu = avmu + mi(az_tmx)
-      !!                avmv = avmv + mj(az_tmx)
-      !!
-      !! ** Action  :   Update avt, avmu and avmv
+      !!                  Kz_tides = az_tmx / max( rn_n2min, N^2 )
+      !!              where az_tmx is a coefficient that specified the 3D space 
+      !!              distribution of the faction of tidal energy taht is used
+      !!              for mixing. Its expression is set in zdf_tmx_init routine,
+      !!              following Simmons et al. 2004.
+      !!                NB: a specific bounding procedure is performed on av_tide
+      !!              so that the input tidal energy is actually almost used. The
+      !!              basic maximum value is 60 cm2/s, but values of 300 cm2/s 
+      !!              can be reached in area where bottom stratification is too 
+      !!              weak.
+      !!
+      !!              - update av_tide in the Indonesian Through Flow area
+      !!              following Koch-Larrouy et al. (2007) parameterisation
+      !!              (see tmx_itf routine).
+      !!
+      !!              - update the model vertical eddy viscosity and diffusivity: 
+      !!                     avt  = avt  +    av_tides
+      !!                     avmu = avmu + mi(av_tides)
+      !!                     avmv = avmv + mj(av_tides)
+      !!
+      !! ** Action  :   avt, avmu, avmv   increased by tidal mixing
       !!
       !! References : Simmons et al. 2004, Ocean Modelling, 6, 3-4, 245-263.
-      !!              St Laurent et al. 2002, GRL, 29, 2106, 10.1029/2002GL015 633.
-      !!              Osborn 1982, JPO, 10, 83-89.
+      !!              Koch-Larrouy et al. 2007, GRL.
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step 
@@ -80,16 +94,19 @@
       !!----------------------------------------------------------------------
 
-      !                             ! Initialization (first time-step only)
-      IF( kt == nit000  )   CALL zdf_tmx_init
-
-      !                             ! First estimation (with n2 bound by rn_n2min) bounded by 60 cm2/s
+      !                        
+      IF( kt == nit000  )   CALL zdf_tmx_init      ! Initialization (first time-step only)
+
+      !                     ! ----------------------- !
+      !                     !  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(:,:,:) )  )
 
-      zkz(:,:) = 0.e0               ! Associated potential energy consummed over the whole water column
+      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 * av_tide(:,:,jk)* tmask(:,:,jk)
       END DO
 
-      DO jj = 1, jpj                ! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
+      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
             IF( zkz(ji,jj) /= 0.e0 )   zkz(ji,jj) = en_tmx(ji,jj) / zkz(ji,jj)
@@ -97,15 +114,15 @@
       END DO
 
-      DO jk = 2, jpkm1              ! Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> av_tide bound by 300 cm2/s
+      DO jk = 2, jpkm1              !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> av_tide bound by 300 cm2/s
          av_tide(:,:,jk) = av_tide(:,:,jk) * MIN( zkz(:,:), 30./6. )   !kz max = 300 cm2/s
       END DO
 
-      IF( kt == nit000 ) THEN       ! diagnose the nergy consumed by av_tide
+      IF( kt == nit000 ) THEN       !* check at first time-step: diagnose the energy consumed by av_tide
          ztpc = 0.e0
          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))*av_tide(ji,jj,jk)*tmask(ji,jj,jk)*tmask_i(ji,jj)
+                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj)   &
+                     &         * MAX( 0.e0, rn2(ji,jj,jk) ) * av_tide(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
                END DO
             END DO
@@ -113,5 +130,5 @@
          ztpc= rau0 / ( rn_tfe * rn_me ) * ztpc
          IF(lwp) WRITE(numout,*) 
-         IF(lwp) WRITE(numout,*) '          N Total power consumption by av_tide : ztpc = ', ztpc * 1.e-12 ,'TW'
+         IF(lwp) WRITE(numout,*) '          N Total power consumption by av_tide    : ztpc = ', ztpc * 1.e-12 ,'TW'
       ENDIF
        
@@ -132,14 +149,7 @@
          END DO
       END DO
-      CALL lbc_lnk( avmu, 'U', 1. )   
-      CALL lbc_lnk( avmv, 'V', 1. )
-       
-      CALL tmx_itf( kt )            ! Indonesian throughflow tidal vertical mixing 
-
-      IF(ln_ctl) THEN
-         CALL prt_ctl_info( clinfo1='tmx  - t: ', ivar1=INT(SUM( avt (:,:,:) ) ) )
-         CALL prt_ctl_info( clinfo1='tmx  - u: ', ivar1=INT(SUM( avmu(:,:,:) ) ),   &
-            &               clinfo2='     - v: ', ivar2=INT(SUM( avmv(:,:,:) ) ) )
-      ENDIF
+      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! lateral boundary condition
+
+      IF(ln_ctl)   CALL prt_ctl(tab3d_1=av_tide ,clinfo1=' tmx - av_tide: ',tab3d_2=avt,clinfo2=' avt: ',ovlap=1,kdim=jpk)
       !
    END SUBROUTINE zdf_tmx
@@ -150,17 +160,17 @@
       !!                  ***  ROUTINE tmx_itf  ***
       !!                   
-      !! ** Purpose :   add to the vertical mixing coefficients the effect of
-      !!      tidal mixing (st laurent et al. 2004) + specificities for the ITF 
-      !!      region (vertical profil of energy (P. Bouruet Aubertot + Theo Gerkema)
-      !!      , q=1 all the energy remains confined)
-      !!
-      !! ** Method  : - tidal-induced vertical mixing is given by:
-      !!      Kz_tides = az_tmx / rn2, where az_tmx is a 2D coefficient set in
-      !!      zdf_tmx_init. This mixing is added to avt, avmu and avmv as
-      !!      follows:  avt  = avt  +    az_tmx
-      !!                avmu = avmu + mi(az_tmx)
-      !!                avmv = avmv + mj(az_tmx)
-      !!
-      !! ** Action  :   Update avt, avmu and avmv
+      !! ** Purpose :   modify the vertical eddy diffusivity coefficients 
+      !!              (av_tide) in the Indonesian Through Flow area (ITF).
+      !!
+      !! ** Method  : - Following Koch-Larrouy et al. (2007), in the ITF defined
+      !!                by msk_itf (read in a file, see tmx_init), the tidal
+      !!                mixing coefficient is computed with :
+      !!                  * q=1 (i.e. all the tidal energy remains trapped in
+      !!                         the area and thus is used for mixing)
+      !!                  * the vertical distribution of the tifal energy is a
+      !!                    proportional to N above the thermocline (d(N^2)/dz > 0)
+      !!                    and to N^2 below the thermocline (d(N^2)/dz < 0)
+      !!
+      !! ** Action  :   av_tide   updated in the ITF area (msk_itf)
       !!
       !! References :  Koch-Larrouy et al. 2007, GRL 
@@ -335,22 +345,34 @@
       !!                     
       !! ** Purpose :   Initialization of the vertical tidal mixing, Reading
-      !!      of M2 and K1 tidal energy in nc files
-      !!
-      !! ** Method  :   Read the namtmx namelist and check the parameters
-      !!      called at the first timestep (nit000), 
-      !!              - Read the M2 and K1 tidal energy in NetCDF files.
+      !!              of M2 and K1 tidal energy in nc files
+      !!
+      !! ** Method  : - Read the namtmx namelist and check the parameters
+      !!
+      !!              - Read the input data in NetCDF files :
+      !!              M2 and K1 tidal energy. The total tidal energy, en_tmx, 
+      !!              is the sum of M2, K1 and S2 energy where S2 is assumed 
+      !!              to be: S2=(1/2)^2 * M2
+      !!              mask_itf, a mask array that determine where substituing 
+      !!              the standard Simmons et al. (2005) formulation with the
+      !!              one of Koch_Larrouy et al. (2007).
+      !!
       !!              - Compute az_tmx, a 3D coefficient that allows to compute
-      !!      the tidal-induced vertical mixing as follows:
-      !!         Kz_tides = avt0 + az_tmx / max( rn_n2min, N^2 )
-      !!      where az_tmx is given by.....
+      !!             the standard tidal-induced vertical mixing as follows:
+      !!                  Kz_tides = az_tmx / max( rn_n2min, N^2 )
+      !!             with az_tmx a bottom intensified coefficient is given by:
+      !!                 az_tmx(z) = en_tmx / ( rau0 * rn_htmx ) * EXP( -(H-z)/rn_htmx )
+      !!                                                  / ( 1. - EXP( - H   /rn_htmx ) ) 
+      !!             where rn_htmx the characteristic length scale of the bottom 
+      !!             intensification, en_tmx the tidal energy, and H the ocean depth
       !!
       !! ** input   :   - Namlist namtmx
-      !!                - NetCDF file  : M2_ORCA2.nc and K1_ORCA2.nc
+      !!                - NetCDF file : M2_ORCA2.nc, K1_ORCA2.nc, and mask_itf.nc
       !!
       !! ** Action  : - Increase by 1 the nstop flag is setting problem encounter
       !!              - defined az_tmx used to compute tidal-induced mixing
-      !!----------------------------------------------------------------------
-      USE iom                                           ! to read input files
-
+      !!
+      !! References : Simmons et al. 2004, Ocean Modelling, 6, 3-4, 245-263.
+      !!              Koch-Larrouy et al. 2007, GRL.
+      !!----------------------------------------------------------------------
       INTEGER ::   ji, jj, jk                           ! dummy loop indices
       INTEGER ::   inum                                 ! temporary logical unit
@@ -361,6 +383,5 @@
       REAL(wp), DIMENSION(jpi,jpj) ::  zhdep                ! Ocean depth 
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zpc          ! power consumption
-
-      !                                   ! * tidal mixing namelist * (namtmx)
+      !!
       NAMELIST/namtmx/ rn_htmx, rn_n2min, rn_tfe, rn_tfe_itf, rn_me
       !!----------------------------------------------------------------------
@@ -408,8 +429,6 @@
       en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * tmask(:,:,1)
 
-      ! Bound total energy required to be <= 1 W/m2
-      ! Vertical structure 
-      ! part independent of the level
-      DO jj = 1, jpj
+      ! Vertical structure (az_tmx)
+      DO jj = 1, jpj                ! part independent of the level
          DO ji = 1, jpi
             zhdep(ji,jj) = fsdepw(ji,jj,mbathy(ji,jj))       ! depth of the ocean
