Index: trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 6497)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 7698)
@@ -171,4 +171,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
+      INTEGER             ::   jk, jj, ji  
       !!----------------------------------------------------------------------
       !
@@ -179,8 +180,15 @@
       !
       IF( kt /= nit000 ) THEN   ! restore before value to compute tke
-         avt (:,:,:) = avt_k (:,:,:) 
-         avm (:,:,:) = avm_k (:,:,:) 
-         avmu(:,:,:) = avmu_k(:,:,:) 
-         avmv(:,:,:) = avmv_k(:,:,:) 
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  avt (ji,jj,jk) = avt_k (ji,jj,jk) 
+                  avm (ji,jj,jk) = avm_k (ji,jj,jk) 
+                  avmu(ji,jj,jk) = avmu_k(ji,jj,jk) 
+                  avmv(ji,jj,jk) = avmv_k(ji,jj,jk) 
+               END DO
+            END DO
+         END DO
       ENDIF 
       !
@@ -189,8 +197,15 @@
       CALL tke_avn      ! now avt, avm, avmu, avmv
       !
-      avt_k (:,:,:) = avt (:,:,:) 
-      avm_k (:,:,:) = avm (:,:,:) 
-      avmu_k(:,:,:) = avmu(:,:,:) 
-      avmv_k(:,:,:) = avmv(:,:,:) 
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  avt_k (ji,jj,jk) = avt (ji,jj,jk) 
+                  avm_k (ji,jj,jk) = avm (ji,jj,jk) 
+                  avmu_k(ji,jj,jk) = avmu(ji,jj,jk) 
+                  avmv_k(ji,jj,jk) = avmv(ji,jj,jk) 
+               END DO
+            END DO
+         END DO
       !
 #if defined key_agrif
@@ -253,4 +268,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 +275,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.
@@ -293,12 +310,29 @@
          !
          !                        !* total energy produce by LC : cumulative sum over jk
-         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * gdepw_n(:,:,1) * e3w_n(:,:,1)
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj, ji)
+         DO jj =1, jpj
+            DO ji=1, jpi
+               zpelc(ji,jj,1) =  MAX( rn2b(ji,jj,1), 0._wp ) * gdepw_n(ji,jj,1) * e3w_n(ji,jj,1)
+            END DO
+         END DO
          DO jk = 2, jpk
-            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw_n(:,:,jk) * e3w_n(:,:,jk)
+!$OMP 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
          zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
-         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
+!$OMP DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  imlc(ji,jj) = mbkt(ji,jj) + 1       ! Initialization to the number of w ocean point (=2 over land)
+               END DO
+            END DO
          DO jk = jpkm1, 2, -1
+!$OMP DO schedule(static) private(jj, ji, zus)
             DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us 
                DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1)
@@ -309,4 +343,5 @@
          END DO
          !                               ! finite LC depth
+!$OMP DO schedule(static) private(jj, ji)
          DO jj = 1, jpj 
             DO ji = 1, jpi
@@ -315,4 +350,5 @@
          END DO
          zcof = 0.016 / SQRT( zrhoa * zcdrag )
+!$OMP 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
@@ -328,4 +364,5 @@
             END DO
          END DO
+!$OMP END PARALLEL
          !
       ENDIF
@@ -338,4 +375,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 +394,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 +411,6 @@
       ENDIF
       !         
+!$OMP PARALLEL
+!$OMP 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 +446,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 DO schedule(static) private(jj, ji)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
@@ -411,4 +453,5 @@
          END DO
       END DO
+!$OMP 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 +460,5 @@
       END DO
       DO jk = 3, jpkm1
+!$OMP DO schedule(static) private(jj, ji)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
@@ -423,4 +467,5 @@
          END DO
       END DO
+!$OMP 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 +474,5 @@
       END DO
       DO jk = jpk-2, 2, -1
+!$OMP DO schedule(static) private(jj, ji)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
@@ -435,4 +481,5 @@
          END DO
       END DO
+!$OMP DO schedule(static) private(jk,jj, ji)
       DO jk = 2, jpkm1                             ! set the minimum value of tke
          DO jj = 2, jpjm1
@@ -442,4 +489,5 @@
          END DO
       END DO
+!$OMP END PARALLEL
 
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -450,4 +498,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 +508,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 +517,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 +596,16 @@
       !
       ! initialisation of interior minimum value (avoid a 2d loop with mikt)
-      zmxlm(:,:,:)  = rmxl_min    
-      zmxld(:,:,:)  = rmxl_min
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zmxlm(ji,jj,jk)  = rmxl_min    
+               zmxld(ji,jj,jk)  = rmxl_min
+            END DO
+         END DO
+      END DO
       !
       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 +615,14 @@
          END DO
       ELSE 
-         zmxlm(:,:,1) = rn_mxl0
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zmxlm(ji,jj,1) = rn_mxl0
+            END DO
+         END DO
       ENDIF
       !
+!$OMP PARALLEL
+!$OMP 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 +636,12 @@
       !                     !* Physical limits for the mixing length
       !
-      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value 
-      zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value
+!$OMP DO schedule(static) private(jj,ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zmxld(ji,jj, 1 ) = zmxlm(ji,jj,1)   ! surface set to the minimum value 
+            zmxld(ji,jj,jpk) = rmxl_min       ! last level  set to the minimum value
+         END DO
+      END DO
+!$OMP END PARALLEL
       !
       SELECT CASE ( nn_mxl )
@@ -578,4 +650,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 +664,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
@@ -602,5 +676,7 @@
          !
       CASE ( 2 )           ! |dk[xml]| bounded by e3t :
+!$OMP PARALLEL
          DO jk = 2, jpkm1         ! from the surface to the bottom :
+!$OMP DO schedule(static) private(jj, ji)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -610,4 +686,5 @@
          END DO
          DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
+!$OMP DO schedule(static) private(jj, ji, zemxl)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -618,7 +695,10 @@
             END DO
          END DO
+!$OMP END PARALLEL
          !
       CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
+!$OMP PARALLEL
          DO jk = 2, jpkm1         ! from the surface to the bottom : lup
+!$OMP DO schedule(static) private(jj, ji)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -628,4 +708,5 @@
          END DO
          DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
+!$OMP DO schedule(static) private(jj, ji)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -634,4 +715,5 @@
             END DO
          END DO
+!$OMP DO schedule(static) private(jk, jj, ji, zemlm, zemlp)
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
@@ -644,10 +726,18 @@
             END DO
          END DO
+!$OMP END PARALLEL
          !
       END SELECT
       !
 # if defined key_c1d
-      e_dis(:,:,:) = zmxld(:,:,:)      ! c1d configuration : save mixing and dissipation turbulent length scales
-      e_mix(:,:,:) = zmxlm(:,:,:)
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               e_dis(ji,jj,jk) = zmxld(ji,jj,jk)      ! c1d configuration : save mixing and dissipation turbulent length scales
+               e_mix(ji,jj,jk) = zmxlm(ji,jj,jk)
+            END DO
+         END DO
+      END DO
 # endif
 
@@ -655,4 +745,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 +759,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 +771,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
@@ -798,17 +891,42 @@
          SELECT CASE( nn_htau )             ! Choice of the depth of penetration
          CASE( 0 )                                 ! constant depth penetration (here 10 meters)
-            htau(:,:) = 10._wp
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  htau(ji,jj) = 10._wp
+               END DO
+            END DO
          CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
-            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )            
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  htau(ji,jj) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   )            
+               END DO
+            END DO
          END SELECT
       ENDIF
       !                               !* set vertical eddy coef. to the background value
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jk,jj,ji)
       DO jk = 1, jpk
-         avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
-         avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
-         avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
-         avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
-      END DO
-      dissl(:,:,:) = 1.e-12_wp
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               avt (ji,jj,jk) = avtb(jk) * wmask (ji,jj,jk)
+               avm (ji,jj,jk) = avmb(jk) * wmask (ji,jj,jk)
+               avmu(ji,jj,jk) = avmb(jk) * wumask(ji,jj,jk)
+               avmv(ji,jj,jk) = avmb(jk) * wvmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+!$OMP END DO NOWAIT
+!$OMP DO schedule(static) private(jk,jj,ji)
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               dissl(ji,jj,jk) = 1.e-12_wp
+            END DO
+         END DO
+      END DO
+!$OMP END PARALLEL
       !                              
       CALL tke_rst( nit000, 'READ' )  !* read or initialize all required files
@@ -830,5 +948,5 @@
      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
      !
-     INTEGER ::   jit, jk   ! dummy loop indices
+     INTEGER ::   jit, jk, jj, ji   ! dummy loop indices
      INTEGER ::   id1, id2, id3, id4, id5, id6   ! local integers
      !!----------------------------------------------------------------------
@@ -857,22 +975,51 @@
            ELSE                                     ! No TKE array found: initialisation
               IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
-              en (:,:,:) = rn_emin * tmask(:,:,:)
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+              DO jk = 1, jpk
+                 DO jj = 1, jpj
+                    DO ji = 1, jpi
+                       en (ji,jj,jk) = rn_emin * tmask(ji,jj,jk)
+                    END DO
+                 END DO
+              END DO
               CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
               !
-              avt_k (:,:,:) = avt (:,:,:)
-              avm_k (:,:,:) = avm (:,:,:)
-              avmu_k(:,:,:) = avmu(:,:,:)
-              avmv_k(:,:,:) = avmv(:,:,:)
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+              DO jk = 1, jpk
+                 DO jj = 1, jpj
+                    DO ji = 1, jpi
+                       avt_k (ji,jj,jk) = avt (ji,jj,jk)
+                       avm_k (ji,jj,jk) = avm (ji,jj,jk)
+                       avmu_k(ji,jj,jk) = avmu(ji,jj,jk)
+                       avmv_k(ji,jj,jk) = avmv(ji,jj,jk)
+                    END DO
+                 END DO
+              END DO
               !
               DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO
            ENDIF
         ELSE                                   !* Start from rest
-           en(:,:,:) = rn_emin * tmask(:,:,:)
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jk,jj,ji)
+           DO jk = 1, jpk
+              DO jj = 1, jpj
+                 DO ji = 1, jpi
+                    en(ji,jj,jk) = rn_emin * tmask(ji,jj,jk)
+                 END DO
+              END DO
+           END DO
+!$OMP END DO NOWAIT
+!$OMP DO schedule(static) private(jk)
            DO jk = 1, jpk                           ! set the Kz to the background value
-              avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
-              avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
-              avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
-              avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
+              DO jj = 1, jpj
+                 DO ji = 1, jpi
+                    avt (ji,jj,jk) = avtb(jk) * wmask (ji,jj,jk)
+                    avm (ji,jj,jk) = avmb(jk) * wmask (ji,jj,jk)
+                    avmu(ji,jj,jk) = avmb(jk) * wumask(ji,jj,jk)
+                    avmv(ji,jj,jk) = avmb(jk) * wvmask(ji,jj,jk)
+                 END DO
+              END DO
            END DO
+!$OMP END PARALLEL
         ENDIF
         !
