Index: branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90
===================================================================
--- branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90	(revision 7171)
+++ branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90	(revision 7340)
@@ -6,10 +6,7 @@
    !! History :  3.6  !  2014-10  (E. Clementi)  Original code
    !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
-   !!   qiao_init      
    !!   zdf_qiao        : compute Qiao parameters
    !!----------------------------------------------------------------------
 
-   USE iom             ! I/O manager library
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! distribued memory computing library
@@ -18,19 +15,16 @@
    USE sbcwave         ! wave module
    USE dom_oce
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)  
    
-   !!----------------------------------------------------------------------
-   !!   qiao_init       : compute QBv: Qiao terms to be added to vertical eddy
-   !!                     diffusivity and viscosity coefficients 
-   !!----------------------------------------------------------------------
-
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   zdf_qiao    ! routine called in zdf_ric
+   PUBLIC zdf_qiao    ! routine called in step
 
-   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:)     :: QBv, QBvu, QBvv
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: qbv, qbvu, qbvv
 
    !! * Substitutions
 #  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
@@ -45,8 +39,8 @@
       !!                     ***  ROUTINE zdf_qiao ***
       !!
-      !! ** Purpose :Compute the Qiao term (QBv) to be added to
+      !! ** Purpose :Compute the Qiao term (qbv) to be added to
       !!             vertical viscosity and diffusivity coeffs.  
       !!
-      !! ** Method  :QBv = alpha * A * Us(0) * exp (3 * k * z)
+      !! ** Method  :qbv = alpha * A * Us(0) * exp (3 * k * z)
       !!             
       !! ** action  :Compute the Qiao wave dependent term 
@@ -56,68 +50,79 @@
       INTEGER, INTENT( in  ) ::  kt   ! ocean time step
       !
-      INTEGER                ::  jj, ji, jk
+      INTEGER :: jj, ji, jk   ! dummy loop indices
       !!---------------------------------------------------------------------
       !
-      !
-      !                                         ! -------------------- !
       IF( kt == nit000 ) THEN                   ! First call kt=nit000 !
-         ALLOCATE(QBv(jpi,jpj,jpk))             ! -------------------- !
-         ALLOCATE(QBvu(jpi,jpj,jpk))
-         ALLOCATE(QBvv(jpi,jpj,jpk))
+         IF( .NOT. ( ln_wave .AND. ln_sdw ) )   &
+            &   CALL ctl_stop ( 'Ask for wave Qiao enhanced turbulence but ln_wave   &
+            &                    and ln_sdw have to be activated')
+         IF( zdf_qiao_alloc() /= 0 )   &
+            &   CALL ctl_stop( 'STOP', 'zdf_qiao : unable to allocate arrays' )
       ENDIF
 
-      QBv (:,:,:) = 0.0
-      QBvu(:,:,:) = 0.0
-      QBvv(:,:,:) = 0.0
-
       !
-      ! Compute the Qiao term Bv (QBv) to be added to
+      ! Compute the Qiao term Bv (qbv) to be added to
       ! vertical viscosity and diffusivity
-      ! QBv = alpha * A * Us(0) * exp (3 * k * z)
+      ! qbv = alpha * A * Us(0) * exp (3 * k * z)
       ! alpha here is set to 1
       !---------------------------------------------------------------------------------
       !
-      IF ( ln_wave ) THEN
-         DO jk = 1, jpk
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                  QBv(ji,jj,jk) = 1.0 * 0.353553 * swh(ji,jj) * tsd2d(ji,jj) *       &
-               &              exp(3.0 * wnum(ji,jj) *                                &                     
-               &              (-MIN( fsdept(ji  ,jj  ,jk) , fsdept(ji+1,jj  ,jk),    &
-               &                     fsdept(ji  ,jj+1,jk) , fsdept(ji+1,jj+1,jk))))
-               END DO
+      DO jk = 1, jpk
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1
+               qbv(ji,jj,jk) = 1.0 * 0.353553 * swh(ji,jj) * tsd2d(ji,jj) *           &
+            &                  EXP(3.0 * wnum(ji,jj) *                                &                     
+            &                  (-MIN( fsdepw(ji  ,jj  ,jk), fsdepw(ji+1,jj  ,jk),     &
+            &                         fsdepw(ji  ,jj+1,jk), fsdepw(ji+1,jj+1,jk))))   &
+            &                          * wmask(ji,jj,jk)
             END DO
          END DO
-
-         QBv(jpi,:,:)=QBv(jpim1,:,:)
-         QBv(:,jpj,:)=QBv(:,jpjm1,:)
-
-         !
-         ! Interpolate Qiao parameter QBv into the grid_U and grid_V
-         !-------------------------------------------------
-         !
-         DO jk = 1, jpk
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                  QBvu(ji,jj,jk) = 0.5 *  umask(ji,jj,jk)  *               &
-               &           ( QBv(ji  ,jj,jk) * tmask(ji  ,jj,jk)           &
-               &           + QBv(ji+1,jj,jk) * tmask(ji+1,jj,jk) )
-                  QBvv(ji,jj,jk) = 0.5 *  vmask(ji,jj,jk)  *               &
-               &           ( QBv(ji,jj  ,jk) * tmask(ji,jj  ,jk)           &
-               &           + QBv(ji,jj+1,jk) * tmask(ji,jj+1,jk) )
-               END DO
+      END DO
+      !
+      CALL lbc_lnk( qbv, 'W', 1. )   ! Lateral boundary conditions
+         
+      !
+      ! Interpolate Qiao parameter qbv into the grid_U and grid_V
+      !----------------------------------------------------------
+      !
+      DO jk = 1, jpk
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1
+               qbvu(ji,jj,jk) = 0.5 * wumask(ji,jj,jk)  *              &  
+            &                  ( qbv(ji,jj,jk) + qbv(ji+1,jj  ,jk) )
+               qbvv(ji,jj,jk) = 0.5 * wvmask(ji,jj,jk)  *              &
+            &                  ( qbv(ji,jj,jk) + qbv(ji  ,jj+1,jk) )
             END DO
          END DO
-         ! 
-         QBvu(jpi,:,:)=QBvu(jpim1,:,:)
-         QBvu(:,jpj,:)=QBvu(:,jpjm1,:)
-         QBvv(jpi,:,:)=QBvv(jpim1,:,:)
-         QBvv(:,jpj,:)=QBvv(:,jpjm1,:)
+      END DO
+      ! 
+      CALL lbc_lnk( qbvu, 'U', 1. ) ; CALL lbc_lnk( qbvv, 'V', 1. )   ! Lateral boundary conditions
 
-        ELSE
-           CALL ctl_stop( 'STOP', 'To use Qiao formulation you have to set: ln_wave=.true.')
-        ENDIF
-        !
+      ! Enhance vertical mixing coeff.         
+      !-------------------------------
+      !
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               avmu(ji,jj,jk) = ( avmu(ji,jj,jk) + qbvu(ji,jj,jk) ) * umask(ji,jj,jk)
+               avmv(ji,jj,jk) = ( avmv(ji,jj,jk) + qbvv(ji,jj,jk) ) * vmask(ji,jj,jk)
+               avt (ji,jj,jk) = ( avt (ji,jj,jk) + qbv (ji,jj,jk) ) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
    END SUBROUTINE zdf_qiao
+
+   INTEGER FUNCTION zdf_qiao_alloc()
+      !!----------------------------------------------------------------------
+      !!                ***  FUNCTION zdf_qiao_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( qbv(jpi,jpj,jpk), qbvu(jpi,jpj,jpk), qbvv(jpi,jpj,jpk),   &
+         &      STAT = zdf_qiao_alloc )
+      !
+      IF( lk_mpp             )  CALL mpp_sum ( zdf_qiao_alloc )
+      IF( zdf_qiao_alloc > 0 )  CALL ctl_warn('zdf_qiao_alloc: allocation of arrays failed.')
+      !
+   END FUNCTION zdf_qiao_alloc
       
    !!======================================================================
