Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 6520)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 6748)
@@ -112,4 +112,5 @@
          IF ( ln_loglayer.AND. .NOT.ln_linssh ) THEN ! "log layer" bottom friction coefficient
 
+!$OMP PARALLEL DO private(jj,ji,ikbt,ztmp)
             DO jj = 1, jpj
                DO ji = 1, jpi
@@ -123,4 +124,5 @@
 ! (ISF)
             IF ( ln_isfcav ) THEN
+!$OMP PARALLEL DO private(jj,ji,ikbt,ztmp)
                DO jj = 1, jpj
                   DO ji = 1, jpi
@@ -135,8 +137,11 @@
             !   
          ELSE
+!$OMP PARALLEL WORKSHARE
             zbfrt(:,:) = bfrcoef2d(:,:)
             ztfrt(:,:) = tfrcoef2d(:,:)
-         ENDIF
-
+!$OMP END PARALLEL WORKSHARE
+         ENDIF
+
+!$OMP PARALLEL DO private(jj,ji,ikbu,ikbv,zvu,zuv,zecu,zecv)
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
@@ -173,4 +178,5 @@
 
          IF( ln_isfcav ) THEN
+!$OMP PARALLEL DO private(jj,ji,ikbu,ikbv,zvu,zuv,zecu,zecv)
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
@@ -266,8 +272,10 @@
       CASE( 0 )
          IF(lwp) WRITE(numout,*) '      free-slip '
+!$OMP PARALLEL WORKSHARE
          bfrua(:,:) = 0.e0
          bfrva(:,:) = 0.e0
          tfrua(:,:) = 0.e0
          tfrva(:,:) = 0.e0
+!$OMP END PARALLEL WORKSHARE
          !
       CASE( 1 )
@@ -296,6 +304,8 @@
          ENDIF
          !
+!$OMP PARALLEL WORKSHARE
          bfrua(:,:) = - bfrcoef2d(:,:)
          bfrva(:,:) = - bfrcoef2d(:,:)
+!$OMP END PARALLEL WORKSHARE
          !
          IF ( ln_isfcav ) THEN
@@ -310,6 +320,8 @@
             ENDIF
             !
+!$OMP PARALLEL WORKSHARE
             tfrua(:,:) = - tfrcoef2d(:,:)
             tfrva(:,:) = - tfrcoef2d(:,:)
+!$OMP END PARALLEL WORKSHARE
          END IF
          !
@@ -371,4 +383,5 @@
          !
          IF( ln_loglayer.AND. ln_linssh ) THEN ! set "log layer" bottom friction once for all
+!$OMP PARALLEL DO private(jj,ji,ikbt,ztmp)
             DO jj = 1, jpj
                DO ji = 1, jpi
@@ -380,4 +393,5 @@
             END DO
             IF ( ln_isfcav ) THEN
+!$OMP PARALLEL DO private(jj,ji,ikbt,ztmp)
                DO jj = 1, jpj
                   DO ji = 1, jpi
@@ -419,4 +433,5 @@
       zmaxtfr = -1.e10_wp    ! initialise tracker for maximum of bottom friction coefficient
       !
+!$OMP PARALLEL DO private(jj,ji,ikbu,ikbv,zfru,zfrv,ictu,ictv,zminbfr,zmaxbfr,zmintfr,zmaxtfr)
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90	(revision 6520)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90	(revision 6748)
@@ -76,4 +76,5 @@
          zavm_evd(:,:,:) = avm(:,:,:)           ! set avm prior to evd application
          !
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
          DO jk = 1, jpkm1 
             DO jj = 2, jpj             ! no vector opt.
@@ -97,4 +98,5 @@
          !
       CASE DEFAULT         ! enhance vertical eddy diffusivity only (if rn2<-1.e-12) 
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
          DO jk = 1, jpkm1
 !!!         WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd   ! agissant sur T SEUL! 
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 6520)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 6748)
@@ -96,6 +96,8 @@
 
       ! w-level of the mixing and mixed layers
+!$OMP PARALLEL WORKSHARE
       nmln(:,:)  = nlb10               ! Initialization to the number of w ocean point
       hmlp(:,:)  = 0._wp               ! here hmlp used as a dummy variable, integrating vertically N^2
+!$OMP END PARALLEL WORKSHARE
       zN2_c = grav * rho_c * r1_rau0   ! convert density criteria into N^2 criteria
       DO jk = nlb10, jpkm1
@@ -110,5 +112,8 @@
       !
       ! w-level of the turbocline and mixing layer (iom_use)
+!$OMP PARALLEL WORKSHARE
       imld(:,:) = mbkt(:,:) + 1        ! Initialization to the number of w ocean point
+!$OMP END PARALLEL WORKSHARE
+
       DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10 
          DO jj = 1, jpj
@@ -119,4 +124,5 @@
       END DO
       ! depth of the mixing and mixed layers
+!$OMP PARALLEL DO schedule(static) private(jj, ji, iiki, iikn)
       DO jj = 1, jpj
          DO ji = 1, jpi
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 6520)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 6748)
@@ -179,8 +179,10 @@
       !
       IF( kt /= nit000 ) THEN   ! restore before value to compute tke
+!$OMP PARALLEL WORKSHARE
          avt (:,:,:) = avt_k (:,:,:) 
          avm (:,:,:) = avm_k (:,:,:) 
          avmu(:,:,:) = avmu_k(:,:,:) 
          avmv(:,:,:) = avmv_k(:,:,:) 
+!$OMP END PARALLEL WORKSHARE
       ENDIF 
       !
@@ -189,8 +191,10 @@
       CALL tke_avn      ! now avt, avm, avmu, avmv
       !
+!$OMP PARALLEL WORKSHARE
       avt_k (:,:,:) = avt (:,:,:) 
       avm_k (:,:,:) = avm (:,:,:) 
       avmu_k(:,:,:) = avmu(:,:,:) 
       avmv_k(:,:,:) = avmv(:,:,:) 
+!$OMP END PARALLEL WORKSHARE
       !
 #if defined key_agrif
@@ -253,4 +257,5 @@
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF ( ln_isfcav ) THEN
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
          DO jj = 2, jpjm1            ! en(mikt(ji,jj))   = rn_emin
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -259,4 +264,5 @@
          END DO
       END IF
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
       DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -295,5 +301,10 @@
          zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * gdepw_n(:,:,1) * e3w_n(:,:,1)
          DO jk = 2, jpk
-            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw_n(:,:,jk) * e3w_n(:,:,jk)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj =1, jpj
+               DO ji=1, jpi
+                  zpelc(ji,jj,jk)  = zpelc(ji,jj,jk-1) + MAX( rn2b(ji,jj,jk), 0._wp ) * gdepw_n(ji,jj,jk) * e3w_n(ji,jj,jk)
+               END DO
+            END DO
          END DO
          !                        !* finite Langmuir Circulation depth
@@ -309,4 +320,5 @@
          END DO
          !                               ! finite LC depth
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
          DO jj = 1, jpj 
             DO ji = 1, jpi
@@ -315,4 +327,5 @@
          END DO
          zcof = 0.016 / SQRT( zrhoa * zcdrag )
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zus, zind, zwlc)
          DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en
             DO jj = 2, jpjm1
@@ -338,4 +351,5 @@
       !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
       !
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
       DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form)
          DO jj = 1, jpjm1
@@ -356,4 +370,5 @@
          ! Note that zesh2 is also computed in the next loop.
          ! We decided to compute it twice to keep code readability and avoid an IF case in the DO loops
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zesh2, zri)
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
@@ -372,4 +387,5 @@
       ENDIF
       !         
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zcof, zzd_up, zzd_lw, zesh2)
       DO jk = 2, jpkm1           !* Matrix and right hand side in en
          DO jj = 2, jpjm1
@@ -405,4 +421,5 @@
       !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
       DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
@@ -411,4 +428,5 @@
          END DO
       END DO
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
       DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -417,4 +435,5 @@
       END DO
       DO jk = 3, jpkm1
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
@@ -423,4 +442,5 @@
          END DO
       END DO
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
       DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -429,4 +449,5 @@
       END DO
       DO jk = jpk-2, 2, -1
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
@@ -435,4 +456,5 @@
          END DO
       END DO
+!$OMP PARALLEL DO schedule(static) private(jk,jj, ji)
       DO jk = 2, jpkm1                             ! set the minimum value of tke
          DO jj = 2, jpjm1
@@ -450,4 +472,5 @@
       
       IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
@@ -459,4 +482,5 @@
          END DO
       ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
+!$OMP PARALLEL DO schedule(static) private(jj, ji, jk)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -467,4 +491,5 @@
          END DO
       ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, ztx2, zty2, ztau, zdif)
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
@@ -545,8 +570,11 @@
       !
       ! initialisation of interior minimum value (avoid a 2d loop with mikt)
+!$OMP PARALLEL WORKSHARE
       zmxlm(:,:,:)  = rmxl_min    
       zmxld(:,:,:)  = rmxl_min
+!$OMP END PARALLEL WORKSHARE
       !
       IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g)
+!$OMP PARALLEL DO schedule(static) private(jj, ji, zraug)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
@@ -556,7 +584,10 @@
          END DO
       ELSE 
+!$OMP PARALLEL WORKSHARE
          zmxlm(:,:,1) = rn_mxl0
+!$OMP END PARALLEL WORKSHARE
       ENDIF
       !
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zrn2)
       DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2)
          DO jj = 2, jpjm1
@@ -570,6 +601,8 @@
       !                     !* Physical limits for the mixing length
       !
+!$OMP PARALLEL WORKSHARE
       zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value 
       zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value
+!$OMP END PARALLEL WORKSHARE
       !
       SELECT CASE ( nn_mxl )
@@ -578,4 +611,5 @@
       ! where wmask = 0 set zmxlm == e3w_n
       CASE ( 0 )           ! bounded by the distance to surface and bottom
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zemxl)
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
@@ -591,4 +625,5 @@
          !
       CASE ( 1 )           ! bounded by the vertical scale factor
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zemxl)
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
@@ -603,4 +638,5 @@
       CASE ( 2 )           ! |dk[xml]| bounded by e3t :
          DO jk = 2, jpkm1         ! from the surface to the bottom :
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -610,4 +646,5 @@
          END DO
          DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
+!$OMP PARALLEL DO schedule(static) private(jj, ji, zemxl)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -621,4 +658,5 @@
       CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
          DO jk = 2, jpkm1         ! from the surface to the bottom : lup
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -628,4 +666,5 @@
          END DO
          DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -634,4 +673,5 @@
             END DO
          END DO
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zemlm, zemlp)
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
@@ -648,6 +688,8 @@
       !
 # if defined key_c1d
+!$OMP PARALLEL WORKSHARE
       e_dis(:,:,:) = zmxld(:,:,:)      ! c1d configuration : save mixing and dissipation turbulent length scales
       e_mix(:,:,:) = zmxlm(:,:,:)
+!$OMP END PARALLEL WORKSHARE
 # endif
 
@@ -655,4 +697,5 @@
       !                     !  Vertical eddy viscosity and diffusivity  (avmu, avmv, avt)
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zsqen, zav)
       DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points
          DO jj = 2, jpjm1
@@ -668,4 +711,5 @@
       CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
       !
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
       DO jk = 2, jpkm1            !* vertical eddy viscosity at wu- and wv-points
          DO jj = 2, jpjm1
@@ -679,4 +723,5 @@
       !
       IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
@@ -804,4 +849,5 @@
       ENDIF
       !                               !* set vertical eddy coef. to the background value
+!$OMP PARALLEL DO schedule(static) private(jk)
       DO jk = 1, jpk
          avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
@@ -857,16 +903,23 @@
            ELSE                                     ! No TKE array found: initialisation
               IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
+!$OMP PARALLEL WORKSHARE
               en (:,:,:) = rn_emin * tmask(:,:,:)
+!$OMP END PARALLEL WORKSHARE
               CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
               !
+!$OMP PARALLEL WORKSHARE
               avt_k (:,:,:) = avt (:,:,:)
               avm_k (:,:,:) = avm (:,:,:)
               avmu_k(:,:,:) = avmu(:,:,:)
               avmv_k(:,:,:) = avmv(:,:,:)
+!$OMP END PARALLEL WORKSHARE
               !
               DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO
            ENDIF
         ELSE                                   !* Start from rest
+!$OMP PARALLEL WORKSHARE
            en(:,:,:) = rn_emin * tmask(:,:,:)
+!$OMP END PARALLEL WORKSHARE
+!$OMP PARALLEL DO schedule(static) private(jk)
            DO jk = 1, jpk                           ! set the Kz to the background value
               avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
