Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DIA/diaprod.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DIA/diaprod.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DIA/diaprod.F90	(revision 9176)
@@ -88,4 +88,5 @@
       IF( iom_use("ut") ) THEN
          z3d(:,:,:) = 0.e0 
+!$OMP PARALLEL DO
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -100,4 +101,5 @@
       IF( iom_use("vt") ) THEN
          z3d(:,:,:) = 0.e0 
+!$OMP PARALLEL DO
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -117,4 +119,5 @@
             END DO
          END DO
+!$OMP PARALLEL DO
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
@@ -129,4 +132,5 @@
       IF( iom_use("us") ) THEN
          z3d(:,:,:) = 0.e0 
+!$OMP PARALLEL DO
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -141,4 +145,5 @@
       IF( iom_use("vs") ) THEN
          z3d(:,:,:) = 0.e0 
+!$OMP PARALLEL DO
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -158,4 +163,5 @@
             END DO
          END DO
+!$OMP PARALLEL DO
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
@@ -170,4 +176,5 @@
       IF( iom_use("urhop") ) THEN
          z3d(:,:,:) = 0.e0 
+!$OMP PARALLEL DO
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -182,4 +189,5 @@
       IF( iom_use("vrhop") ) THEN
          z3d(:,:,:) = 0.e0 
+!$OMP PARALLEL DO
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -199,4 +207,5 @@
             END DO
          END DO
+!$OMP PARALLEL DO
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 9176)
@@ -131,12 +131,12 @@
       REAL(wp)                     ::   zztmp, zztmpx, zztmpy   ! 
       !!
-      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d      ! 2D workspace
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace
+      REAL(wp),  DIMENSION(jpi,jpj    ) :: z2d      ! 2D workspace
+      REAL(wp),  DIMENSION(jpi,jpj,jpk) :: z3d      ! 3D workspace
       !!----------------------------------------------------------------------
       ! 
       IF( nn_timing == 1 )   CALL timing_start('dia_wri')
       ! 
-      CALL wrk_alloc( jpi , jpj      , z2d )
-      CALL wrk_alloc( jpi , jpj, jpk , z3d )
+!     CALL wrk_alloc( jpi , jpj      , z2d )
+!     CALL wrk_alloc( jpi , jpj, jpk , z3d )
       !
       ! Output the initial state and forcings
@@ -176,4 +176,5 @@
       CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity
       IF ( iom_use("sbs") ) THEN
+!$OMP PARALLEL DO PRIVATE(jkbot)
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -187,4 +188,5 @@
       IF ( iom_use("taubot") ) THEN                ! bottom stress
          z2d(:,:) = 0._wp
+!$OMP PARALLEL DO PRIVATE(zztmpx, zztmpy)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -204,4 +206,5 @@
       CALL iom_put(  "ssu", un(:,:,1)         )    ! surface i-current
       IF ( iom_use("sbu") ) THEN
+!$OMP PARALLEL DO PRIVATE(jkbot)
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -221,4 +224,5 @@
       CALL iom_put(  "ssv", vn(:,:,1)         )    ! surface j-current
       IF ( iom_use("sbv") ) THEN
+!$OMP PARALLEL DO PRIVATE(jkbot)
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -239,4 +243,5 @@
          ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
          z2d(:,:) = rau0 * e12t(:,:)
+!$OMP PARALLEL DO
          DO jk = 1, jpk
             z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
@@ -258,4 +263,5 @@
 
       IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
+!$OMP PARALLEL DO PRIVATE(zztmp, zztmpx, zztmpy)
          DO jj = 2, jpjm1                                    ! sst gradient
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -269,5 +275,8 @@
          CALL lbc_lnk( z2d, 'T', 1. )
          CALL iom_put( "sstgrad2",  z2d               )    ! square of module of sst gradient
-         z2d(:,:) = SQRT( z2d(:,:) )
+!$OMP PARALLEL DO
+         DO jj = 1, jpj
+            z2d(:,jj) = SQRT( z2d(:,jj) )
+         ENDDO
          CALL iom_put( "sstgrad" ,  z2d               )    ! module of sst gradient
       ENDIF
@@ -276,4 +285,5 @@
       IF( iom_use("heatc") ) THEN
          z2d(:,:)  = 0._wp 
+!$OMP PARALLEL DO REDUCTION(+:z2d)
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -288,4 +298,5 @@
       IF( iom_use("saltc") ) THEN
          z2d(:,:)  = 0._wp 
+!$OMP PARALLEL DO REDUCTION(+:z2d)
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -300,4 +311,5 @@
       IF ( iom_use("eken") ) THEN
          rke(:,:,jk) = 0._wp                               !      kinetic energy 
+!$OMP PARALLEL DO PRIVATE(zztmp, zztmpx, zztmpy) 
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -325,4 +337,5 @@
       IF( iom_use("u_masstr") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
          z3d(:,:,jpk) = 0.e0
+!$OMP PARALLEL DO
          DO jk = 1, jpkm1
             z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * fse3u(:,:,jk) * umask(:,:,jk)
@@ -333,4 +346,5 @@
       IF( iom_use("u_heattr") ) THEN
          z2d(:,:) = 0.e0 
+!$OMP PARALLEL DO REDUCTION(+:z2d)
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -346,4 +360,5 @@
       IF( iom_use("u_salttr") ) THEN
          z2d(:,:) = 0.e0 
+!$OMP PARALLEL DO REDUCTION(+:z2d)
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -360,4 +375,5 @@
       IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
          z3d(:,:,jpk) = 0.e0
+!$OMP PARALLEL DO
          DO jk = 1, jpkm1
             z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) * vmask(:,:,jk)
@@ -368,4 +384,5 @@
       IF( iom_use("v_heattr") ) THEN
          z2d(:,:) = 0.e0 
+!$OMP PARALLEL DO REDUCTION(+:z2d)
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -381,4 +398,5 @@
       IF( iom_use("v_salttr") ) THEN
          z2d(:,:) = 0.e0 
+!$OMP PARALLEL DO REDUCTION(+:z2d)
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -392,6 +410,6 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi , jpj      , z2d )
-      CALL wrk_dealloc( jpi , jpj, jpk , z3d )
+!     CALL wrk_dealloc( jpi , jpj      , z2d )
+!     CALL wrk_dealloc( jpi , jpj, jpk , z3d )
       !
       IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 9176)
@@ -172,6 +172,7 @@
       fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1)
       fsdepw_b(:,:,1) = 0.0_wp
-
+!$OPM PARALLEL
       DO jk = 2, jpk
+!$OMP DO PRIVATE(zcoef)
          DO jj = 1,jpj
             DO ji = 1,jpi
@@ -189,5 +190,7 @@
             END DO
          END DO
+!$OMP END DO
       END DO
+!$OPM PARALLEL
 
       ! Before depth and Inverse of the local depth of the water column at u- and v- points
@@ -214,4 +217,5 @@
          ENDIF
          IF ( ln_vvl_zstar_at_eqtor ) THEN
+!$OMP PARALLEL DO
             DO jj = 1, jpj
                DO ji = 1, jpi
@@ -273,5 +277,5 @@
       !!----------------------------------------------------------------------
       REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zht, z_scale, zwu, zwv, zhdiv
+      REAL(wp),     DIMENSION(jpi, jpj  ) :: zht, z_scale, zwu, zwv, zhdiv
       !! * Arguments
       INTEGER, INTENT( in )                  :: kt                    ! time step
@@ -285,5 +289,5 @@
       !!----------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt')
-      CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
+!     CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
       CALL wrk_alloc( jpi, jpj, jpk, ze3t                     )
 
@@ -306,10 +310,17 @@
                                                        ! z_star coordinate and barotropic z-tilde part !
       !                                                ! --------------------------------------------- !
-
-      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
+!$OMP PARALLEL
+!$OMP DO
+      DO jj = 1, jpj 
+         z_scale(:,jj) = ( ssha(:,jj) - sshb(:,jj) ) * ssmask(:,jj) / ( ht_0(:,jj) + sshn(:,jj) + 1. - ssmask(:,jj) )
+      END DO
+!$OMP END DO
+!$OMP DO
       DO jk = 1, jpkm1
          ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0)
          fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
       END DO
+!$OMP END DO
+!$OMP END PARALLEL
 
       IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
@@ -323,9 +334,14 @@
          zhdiv(:,:) = 0.
          zht(:,:)   = 0.
+!$OMP PARALLEL DO REDUCTION(+:zhdiv, zht)
          DO jk = 1, jpkm1
             zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
             zht  (:,:) = zht  (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
          END DO
-         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
+!$OMP DO
+      DO jj = 1, jpj
+         zhdiv(:,jj) = zhdiv(:,jj) / ( zht(:,jj) + 1. - tmask_i(:,jj) )
+      END DO
+!$OMP END DO
 
          ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
@@ -333,4 +349,5 @@
          IF( ln_vvl_ztilde ) THEN
             IF( kt .GT. nit000 ) THEN
+!$OMP PARALLEL DO
                DO jk = 1, jpkm1
                   hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
@@ -347,8 +364,10 @@
          ! ----------------------------------
          IF( ln_vvl_ztilde ) THEN     ! z_tilde case
+!$OMP PARALLEL DO
             DO jk = 1, jpkm1
                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
             END DO
          ELSE                         ! layer case
+!$OMP PARALLEL DO
             DO jk = 1, jpkm1
                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
@@ -359,4 +378,5 @@
          ! ------------------
          IF( ln_vvl_ztilde ) THEN
+!$OMP PARALLEL DO
             DO jk = 1, jpk
                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
@@ -369,4 +389,6 @@
          zwv(:,:) = 0.0_wp
          ! a - first derivative: diffusive fluxes
+!$OMP PARALLEL 
+!$OMP DO
          DO jk = 1, jpkm1
             DO jj = 1, jpjm1
@@ -376,4 +398,12 @@
                   vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
                                   & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
+               END DO
+            END DO
+         END DO
+!$OMP END DO
+!$OMP DO REDUCTION(+:zwu, zwv)
+         DO jk = 1, jpkm1
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
                   zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
                   zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
@@ -381,5 +411,8 @@
             END DO
          END DO
+!$OMP END DO
+!$OMP END PARALLEL
          ! b - correction for last oceanic u-v points
+!$OMP PARALLEL DO
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -389,4 +422,5 @@
          END DO
          ! c - second derivative: divergence of diffusive fluxes
+!$OMP PARALLEL DO
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -413,12 +447,33 @@
          ENDIF
          CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp )
-         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
-
+!$OMP PARALLEL DO
+         DO jk = 1, jpk
+            tilde_e3t_a(:,:,jk) = tilde_e3t_b(:,:,jk) + z2dt * tmask(:,:,jk) * tilde_e3t_a(:,:,jk)
+         ENDDO
          ! Maximum deformation control
          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
          ze3t(:,:,jpk) = 0.0_wp
-         DO jk = 1, jpkm1
-            ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
-         END DO
+!$OMP PARALLEL DO
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+!dir$ IVDEP
+               DO ji = 1, jpi
+                  ze3t(ji, jj ,jk) = tilde_e3t_a(ji, jj ,jk) / e3t_0(ji, jj ,jk) * tmask(ji, jj ,jk) * tmask_i(ji, jj)
+               ENDDO
+             ENDDO
+         END DO
+!$OMP PARALLEL
+!$OMP DO
+         DO jk = 1, jpkm1
+                  ze3t(:, : ,jk) = tilde_e3t_a(:, : ,jk) / e3t_0(:, : ,jk) * tmask_i(:, :)
+         END DO
+!$OMP END DO
+!$OMP DO
+         DO jk = 1, jpkm1 
+                  ze3t(:, : ,jk) = ze3t(:, : ,jk) * tmask(:, : ,jk) 
+         END DO
+!$OMP END DO
+!$OMP END PARALLEL
+
          z_tmax = MAXVAL( ze3t(:,:,:) )
          IF( lk_mpp )   CALL mpp_max( z_tmax )                 ! max over the global domain
@@ -448,10 +503,13 @@
          ! - ML - end test
          ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
-         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
-         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
-
+!$OMP PARALLEL DO
+         DO jk =1, jpk
+            tilde_e3t_a(:,:,jk) = MIN( tilde_e3t_a(:,:,jk),   rn_zdef_max * e3t_0(:,:,jk) )
+            tilde_e3t_a(:,:,jk) = MAX( tilde_e3t_a(:,:,jk), - rn_zdef_max * e3t_0(:,:,jk) )
+         ENDDO
          !
          ! "tilda" change in the after scale factor
          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!$OMP PARALLEL DO
          DO jk = 1, jpkm1
             dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
@@ -464,8 +522,10 @@
          !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
          zht(:,:) = 0.
+!$OMP PARALLEL DO REDUCTION(+:zht)
          DO jk = 1, jpkm1
             zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
          z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
+!$OMP PARALLEL DO
          DO jk = 1, jpkm1
             dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
@@ -476,4 +536,5 @@
       IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
       !                                           ! ---baroclinic part--------- !
+!$OMP PARALLEL DO
          DO jk = 1, jpkm1
             fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
@@ -491,4 +552,5 @@
          !
          zht(:,:) = 0.0_wp
+!$OMP PARALLEL DO REDUCTION(+:zht)
          DO jk = 1, jpkm1
             zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
@@ -499,4 +561,5 @@
          !
          zht(:,:) = 0.0_wp
+!$OMP PARALLEL DO REDUCTION(+:zht)
          DO jk = 1, jpkm1
             zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk)
@@ -507,4 +570,5 @@
          !
          zht(:,:) = 0.0_wp
+!$OMP PARALLEL DO REDUCTION(+:zht)
          DO jk = 1, jpkm1
             zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk)
@@ -540,4 +604,5 @@
       hu_a(:,:) = 0._wp                        ! Ocean depth at U-points
       hv_a(:,:) = 0._wp                        ! Ocean depth at V-points
+
       DO jk = 1, jpkm1
          hu_a(:,:) = hu_a(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)
@@ -545,8 +610,11 @@
       END DO
       !                                        ! Inverse of the local depth
-      hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
-      hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
-
-      CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
+!$OMP PARALLEL DO
+      DO jj = 1, jpj
+         hur_a(:,jj) = 1._wp / ( hu_a(:,jj) + 1._wp - umask_i(:,jj) ) * umask_i(:,jj)
+         hvr_a(:,jj) = 1._wp / ( hv_a(:,jj) + 1._wp - vmask_i(:,jj) ) * vmask_i(:,jj)
+      ENDDO
+
+!     CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
       CALL wrk_dealloc( jpi, jpj, jpk, ze3t                     )
 
@@ -603,6 +671,9 @@
             tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
          ELSE
-            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
+!$OMP PARALLEL DO
+            DO jk = 1,jpk
+               tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
             &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
+            ENDDO
          ENDIF
          tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
@@ -636,6 +707,7 @@
       fsdepw_n(:,:,1) = 0.0_wp
       fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
-
+!$OMP PARALLEL SHARED(jk)
       DO jk = 2, jpk
+!$OMP DO PRIVATE(zcoef)
          DO jj = 1,jpj
             DO ji = 1,jpi
@@ -649,5 +721,7 @@
             END DO
          END DO
+!$OMP END DO
       END DO
+!$OMP END PARALLEL
 
       ! Local depth and Inverse of the local depth of the water column at u- and v- points
@@ -663,4 +737,5 @@
       ! --------------------------------------------
       ht(:,:) = 0.
+!$OMP PARALLEL DO REDUCTION(+:ht)
       DO jk = 1, jpkm1
          ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
@@ -705,4 +780,5 @@
          !               ! ------------------------------------- !
          ! horizontal surface weighted interpolation
+!$OMP PARALLEL DO
          DO jk = 1, jpk
             DO jj = 1, jpjm1
@@ -718,9 +794,13 @@
          ! boundary conditions
          CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
-         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
+!$OMP PARALLEL DO
+         DO jk = 1, jpk
+            pe3_out(:,:,jk) = pe3_out(:,:,jk) + e3u_0(:,:,jk)
+         ENDDO
          !               ! ------------------------------------- !
       CASE( 'V' )        ! interpolation from T-point to V-point !
          !               ! ------------------------------------- !
          ! horizontal surface weighted interpolation
+!$OMP PARALLEL DO
          DO jk = 1, jpk
             DO jj = 1, jpjm1
@@ -736,9 +816,13 @@
          ! boundary conditions
          CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
-         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
+!$OMP PARALLEL DO
+         DO jk = 1, jpk
+            pe3_out(:,:,jk) = pe3_out(:,:,jk) + e3v_0(:,:,jk)
+         ENDDO
          !               ! ------------------------------------- !
       CASE( 'F' )        ! interpolation from U-point to F-point !
          !               ! ------------------------------------- !
          ! horizontal surface weighted interpolation
+!$OMP PARALLEL DO
          DO jk = 1, jpk
             DO jj = 1, jpjm1
@@ -754,5 +838,8 @@
          ! boundary conditions
          CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )
-         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
+!$OMP PARALLEL DO
+         DO jk = 1, jpk
+            pe3_out(:,:,jk) = pe3_out(:,:,jk) + e3f_0(:,:,jk)
+         ENDDO
          !               ! ------------------------------------- !
       CASE( 'W' )        ! interpolation from T-point to W-point !
@@ -761,4 +848,5 @@
          pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
          ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
+!$OMP PARALLEL DO
          DO jk = 2, jpk
             pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )   &
@@ -771,4 +859,5 @@
          pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
          ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
+!$OMP PARALLEL DO
          DO jk = 2, jpk
             pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )   &
@@ -781,4 +870,5 @@
          pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
          ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
+!$OMP PARALLEL DO
          DO jk = 2, jpk
             pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )   &
@@ -857,4 +947,5 @@
                IF(lwp) write(numout,*) 'Compute scale factor from sshn'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
+!$OMP PARALLEL DO
                DO jk=1,jpk
                   fse3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 9176)
@@ -105,6 +105,6 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~   NOT optimal for auto-tasking case'
       ENDIF
-
-      !                                                ! ===============
+      !
+!$OMP PARALLEL DO                                       ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab
          !                                             ! ===============
@@ -287,5 +287,5 @@
       ENDIF
 
-      !                                                ! ===============
+!$OMP PARALLEL DO                                      ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab
          !                                             ! ===============
Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90	(revision 9176)
@@ -77,5 +77,5 @@
       INTEGER  ::   ji, jj, jk                  ! dummy loop indices
       REAL(wp) ::   zua, zva, zbt, ze2u, ze2v   ! temporary scalar
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zcu, zcv
+      REAL(wp),      DIMENSION(jpi,jpj  ) :: zcu, zcv
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zuf, zut, zlu, zlv
       !!----------------------------------------------------------------------
@@ -83,5 +83,5 @@
       IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_bilap')
       !
-      CALL wrk_alloc( jpi, jpj,      zcu, zcv           )
+!     CALL wrk_alloc( jpi, jpj,      zcu, zcv           )
       CALL wrk_alloc( jpi, jpj, jpk, zuf, zut, zlu, zlv ) 
       !
@@ -102,36 +102,44 @@
       zlv(:,:,:) = 0._wp
 
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
+         !
          ! Laplacian
          ! ---------
 
          IF( ln_sco .OR. ln_zps ) THEN   ! s-coordinate or z-coordinate with partial steps
-            zuf(:,:,jk) = rotb(:,:,jk) * fse3f(:,:,jk)
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zlu(ji,jj,jk) = - ( zuf(ji,jj,jk) - zuf(ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   &
-                     &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj,jk) ) / e1u(ji,jj)
-   
-                  zlv(ji,jj,jk) = + ( zuf(ji,jj,jk) - zuf(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   &
-                     &         + ( hdivb(ji,jj+1,jk) - hdivb(ji,jj,jk) ) / e2v(ji,jj)
-               END DO
+!$OMP PARALLEL DO
+            !                                                ! ===============
+            DO jk = 1, jpkm1                                 ! Horizontal slab
+               !                                             ! ===============
+               ! Laplacian
+               ! ---------
+                  zuf(:,:,jk) = rotb(:,:,jk) * fse3f(:,:,jk)
+                  DO jj = 2, jpjm1
+                     DO ji = fs_2, fs_jpim1   ! vector opt.
+                        zlu(ji,jj,jk) = - ( zuf(ji,jj,jk) - zuf(ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   &
+                           &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj,jk) ) / e1u(ji,jj)
+   
+                        zlv(ji,jj,jk) = + ( zuf(ji,jj,jk) - zuf(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   &
+                           &         + ( hdivb(ji,jj+1,jk) - hdivb(ji,jj,jk) ) / e2v(ji,jj)
+                     END DO
+                  END DO
             END DO
          ELSE                            ! z-coordinate - full step
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zlu(ji,jj,jk) = - ( rotb (ji  ,jj,jk) - rotb (ji,jj-1,jk) ) / e2u(ji,jj)   &
-                     &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj  ,jk) ) / e1u(ji,jj)
-   
-                  zlv(ji,jj,jk) = + ( rotb (ji,jj  ,jk) - rotb (ji-1,jj,jk) ) / e1v(ji,jj)   &
-                     &         + ( hdivb(ji,jj+1,jk) - hdivb(ji  ,jj,jk) ) / e2v(ji,jj)
+!$OMP PARALLEL DO
+           !                                                ! ===============
+            DO jk = 1, jpkm1                                 ! Horizontal slab
+               !                                             ! ===============
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     zlu(ji,jj,jk) = - ( rotb (ji  ,jj,jk) - rotb (ji,jj-1,jk) ) / e2u(ji,jj)   &
+                        &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj  ,jk) ) / e1u(ji,jj)
+   
+                     zlv(ji,jj,jk) = + ( rotb (ji,jj  ,jk) - rotb (ji-1,jj,jk) ) / e1v(ji,jj)   &
+                        &         + ( hdivb(ji,jj+1,jk) - hdivb(ji  ,jj,jk) ) / e2v(ji,jj)
+                  END DO  
                END DO  
-            END DO  
+            END DO
          ENDIF
-      END DO
       CALL lbc_lnk( zlu, 'U', -1. )   ;   CALL lbc_lnk( zlv, 'V', -1. )   ! Boundary conditions
-
-         
+!$OMP PARALLE DO PRIVATE(zcu, zcv, zbt)
       DO jk = 1, jpkm1
    
@@ -145,6 +153,10 @@
          
          ! Contravariant "laplacian"
-         zcu(:,:) = e1u(:,:) * zlu(:,:,jk)
-         zcv(:,:) = e2v(:,:) * zlv(:,:,jk)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zcu(ji,jj) = e1u(ji,jj) * zlu(ji,jj,jk)
+               zcv(ji,jj) = e2v(ji,jj) * zlv(ji,jj,jk)
+            END DO
+         END DO
          
          ! Laplacian curl ( * e3f if s-coordinates or z-coordinate with partial steps)
@@ -180,5 +192,5 @@
       CALL lbc_lnk( zuf, 'F', 1. )
       CALL lbc_lnk( zut, 'T', 1. )
-
+!OMP PARALLEL DO PRIVATE(ze2u, ze2v, zua, zva)
       DO jk = 1, jpkm1      
    
@@ -205,5 +217,5 @@
       END DO                                           !   End of slab
       !                                                ! ===============
-      CALL wrk_dealloc( jpi, jpj,      zcu, zcv           )
+!     CALL wrk_dealloc( jpi, jpj,      zcu, zcv           )
       CALL wrk_dealloc( jpi, jpj, jpk, zuf, zut, zlu, zlv ) 
       !
Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 9176)
@@ -222,10 +222,10 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   zx1, zy1, zfact2, zx2, zy2   ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz
+      REAL(wp), DIMENSION(jpi, jpj) :: zwx, zwy, zwz
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('vor_ene')
       !
-      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz ) 
+!     CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz ) 
       !
       IF( kt == nit000 ) THEN
@@ -237,5 +237,5 @@
       zfact2 = 0.5 * 0.5      ! Local constant initialization
 
-!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
+!$OMP PARALLEL DO PRIVATE( zwx, zwy, zwz, zy1, zy2, zx1, zx2 )
       !                                                ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab
@@ -292,5 +292,5 @@
       END DO                                           !   End of slab
       !                                                ! ===============
-      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz ) 
+!     CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz ) 
       !
       IF( nn_timing == 1 )  CALL timing_stop('vor_ene')
@@ -350,5 +350,5 @@
       zfact2 = 0.5 * 0.5
 
-!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww )
+!!!!$OMP PARALLEL DO PRIVATE( zwx, zwy, zwz, zww, zy1, zy2, zx1, zx2, zua, zva, zcua, zcva)
       !                                                ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab
@@ -466,5 +466,5 @@
       zfact1 = 0.5 * 0.25      ! Local constant initialization
 
-!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
+!!!!$OMP PARALLEL DO PRIVATE( zwx, zwy, zwz, zuav, zvau )
       !                                                ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab
@@ -599,4 +599,5 @@
 
          IF( ln_dynvor_een_old ) THEN ! original formulation
+!$OMP PARALLEL DO PRIVATE(ze3)
             DO jk = 1, jpk
                DO jj = 1, jpjm1
@@ -609,4 +610,5 @@
             END DO
          ELSE ! new formulation from NEMO 3.6
+!$OMP PARALLEL DO PRIVATE(ze3, zmsk)
             DO jk = 1, jpk
                DO jj = 1, jpjm1
@@ -628,5 +630,5 @@
 
       
-!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse )
+!!!!!$OMP PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse, zua, zva )
       !                                                ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab
@@ -637,7 +639,7 @@
          SELECT CASE( kvor )      ! vorticity considered
          CASE ( 1 )                                                ! planetary vorticity (Coriolis)
-            zwz(:,:) = ff(:,:)      * ze3f(:,:,jk)
+               zwz(:,:) = ff(:,:)      * ze3f(:,:,jk)
          CASE ( 2 )                                                ! relative  vorticity
-            zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk)
+               zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk)
          CASE ( 3 )                                                ! metric term
             DO jj = 1, jpjm1
@@ -650,5 +652,5 @@
             CALL lbc_lnk( zwz, 'F', 1. )
         CASE ( 4 )                                                ! total (relative + planetary vorticity)
-            zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk)
+              zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk)
          CASE ( 5 )                                                ! total (coriolis + metric)
             DO jj = 1, jpjm1
Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 9176)
@@ -105,5 +105,5 @@
       INTEGER  ::   ii0, ii1, iku   ! temporary integer
       INTEGER  ::   ij0, ij1, ikv   ! temporary integer
-      REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_16, zcofw ! local scalars
+      REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_16, zcofw, zcofwa ! local scalars
       REAL(wp) ::   zci, zfi, zau, zbu, zai, zbi   !   -      -
       REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      -
@@ -131,4 +131,5 @@
          zwz(:,:,:) = 0._wp
          !
+!$OMP PARALLEL DO
          DO jk = 1, jpk             !==   i- & j-gradient of density   ==!
             DO jj = 1, jpjm1
@@ -139,5 +140,7 @@
             END DO
          END DO
+
          IF( ln_zps ) THEN                           ! partial steps correction at the bottom ocean level
+!$OMP PARALLEL DO
             DO jj = 1, jpjm1
                DO ji = 1, jpim1
@@ -148,4 +151,5 @@
          ENDIF
          IF( ln_zps .AND. ln_isfcav ) THEN           ! partial steps correction at the bottom ocean level
+!$OMP PARALLEL DO
             DO jj = 1, jpjm1
                DO ji = 1, jpim1
@@ -158,4 +162,5 @@
          !==   Local vertical density gradient at T-point   == !   (evaluated from N^2)
          ! interior value
+!$OMP PARALLEL DO
          DO jk = 2, jpkm1
             !                                ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
@@ -171,4 +176,5 @@
          IF ( ln_isfcav ) THEN
             ! if isf need to overwrite the interior value at at the first ocean point
+!$OMP PARALLEL DO
             DO jj = 1, jpjm1
                DO ji = 1, jpim1
@@ -186,4 +192,5 @@
          !
          IF ( ln_isfcav ) THEN
+!$OMP PARALLEL DO
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -195,4 +202,5 @@
             ENDDO
          ELSE
+!$OMP PARALLEL DO
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -202,4 +210,5 @@
             ENDDO
          END IF
+!$OMP PARALLEL DO PRIVATE(zau, zav, zbu, zbv, zbu, zbv, zfi, zfj, zdepu, zdepv)
          DO jk = 2, jpkm1                            !* Slopes at u and v points
             DO jj = 2, jpjm1
@@ -243,7 +252,9 @@
             END DO
          END DO
+!$OMP END PARALLEL DO
          CALL lbc_lnk( zwz, 'U', -1. )   ;   CALL lbc_lnk( zww, 'V', -1. )      ! lateral boundary conditions
          !
          !                                            !* horizontal Shapiro filter
+!$OMP PARALLEL DO
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1, MAX(1, jpj-3)                        ! rows jj=2 and =jpjm1 only
@@ -287,9 +298,10 @@
             END DO
          END DO
-
+!$OMP END PARALLEL DO
 
          ! II.  slopes at w point           | wslpi = mij( d/di( prd ) / d/dz( prd )
          ! ===========================      | wslpj = mij( d/dj( prd ) / d/dz( prd )
          !
+!$OMP PARALLEL DO PRIVATE(zbw, zci, zcj, zai, zaj, zbi, zbj, zfk, zck)
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
@@ -329,16 +341,18 @@
             END DO
          END DO
+!$OMP END PARALLEL DO
          CALL lbc_lnk( zwz, 'T', -1. )   ;    CALL lbc_lnk( zww, 'T', -1. )      ! lateral boundary conditions
          !
          !                                           !* horizontal Shapiro filter
+!$OMP PARALLEL DO PRIVATE(zcofwa, zcofw, zck)
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1, MAX(1, jpj-3)                        ! rows jj=2 and =jpjm1 only
                DO ji = 2, jpim1
-                  zcofw = tmask(ji,jj,jk) * z1_16
+                  zcofwa = tmask(ji,jj,jk) * z1_16
                   wslpi(ji,jj,jk) = (          zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
                        &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     &
                        &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     &
                        &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   &
-                       &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * zcofw
+                       &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * zcofwa
 
                   wslpj(ji,jj,jk) = (          zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
@@ -346,5 +360,5 @@
                        &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     &
                        &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   &
-                       &                + 4.*  zww(ji  ,jj  ,jk)                         ) * zcofw
+                       &                + 4.*  zww(ji  ,jj  ,jk)                         ) * zcofwa
                END DO
             END DO
@@ -422,4 +436,5 @@
          ! set the slope of diffusion to the slope of s-surfaces 
          !      ( c a u t i o n : minus sign as fsdep has positive value ) 
+!$OMP PARALLEL DO
          DO jj = 2, jpjm1 
             DO ji = fs_2, fs_jpim1   ! vector opt. 
@@ -430,5 +445,5 @@
             END DO 
          END DO 
-
+!$OMP PARALLEL DO
          DO jk = 2, jpk 
             DO jj = 2, jpjm1 
@@ -746,4 +761,5 @@
       !
       !                                            !==   surface mixed layer mask   !
+!$OMP PARALLEL DO PRIVATE(ik)
       DO jk = 1, jpk                               ! =1 inside the mixed layer, =0 otherwise
          DO jj = 1, jpj
@@ -770,4 +786,6 @@
       !-----------------------------------------------------------------------
       !
+!$OMP PARALLEL DO PRIVATE(iku, ikv, zbu, zbv, zau, zav, ik, ikm1, zbw, &
+!$OMP&                    zci, zcj, zai, zaj, zbi, zbj)
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
@@ -872,4 +890,5 @@
             ! set the slope of diffusion to the slope of s-surfaces
             !      ( c a u t i o n : minus sign as fsdep has positive value )
+!$OMP PARALLEL DO
             DO jk = 1, jpk
                DO jj = 2, jpjm1
Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/SOL/solpcg.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/SOL/solpcg.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/SOL/solpcg.F90	(revision 9176)
@@ -94,4 +94,5 @@
       REAL(wp), DIMENSION(2) ::   zsum
       REAL(wp), POINTER, DIMENSION(:,:) ::   zgcr
+      REAL(wp),     DIMENSION(jpi, jpj) ::   tmp1, tmp2
       !!----------------------------------------------------------------------
       !
@@ -109,4 +110,5 @@
       ! gcr   = gcb-a.gcx
       ! gcdes = gcr
+!$OMP PARALLEL DO PRIVATE(zgcad)
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -120,11 +122,21 @@
          END DO
       END DO
+!$OMP END PARALLEL DO 
 
       ! rnorme = (gcr,gcr)
-      rnorme = glob_sum(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  )
-
+      tmp1 = 0.
+!$OMP PARALLEL DO 
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            tmp1(ji, jj) = gcr(ji, jj) * gcdmat(ji, jj) * gcr(ji, jj)
+         END DO
+      END DO      
+!$OMP END PARALLEL DO
+      rnorme = glob_sum(  tmp1(:,:)  )
       CALL lbc_lnk( gcdes, c_solver_pt, 1. )   ! lateral boundary condition
 
       ! gccd = matrix . gcdes
+      gccd = 0.
+!$OMP PARALLEL DO 
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -134,11 +146,19 @@
          END DO
       END DO 
-
+!$OMP END PARALLEL DO
       ! alph = (gcr,gcr)/(gcdes,gccd)
-      radd = glob_sum(  gcdes(:,:) * gcdmat(:,:) * gccd(:,:)  )
+!$OMP PARALLEL DO 
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            tmp1(ji, jj) = gcdes(ji, jj) * gcdmat(ji, jj) * gccd(ji, jj)
+         END DO
+      END DO      
+!$OMP END PARALLEL DO
+      radd = glob_sum(  tmp1  )
       alph = rnorme /radd
 
       ! gcx = gcx + alph * gcdes
       ! gcr = gcr - alph * gccd
+!$OMP PARALLEL DO
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -147,5 +167,5 @@
          END DO
       END DO
-
+!$OMP END PARALLEL DO
       ! Algorithm wtih Eijkhout rearrangement
       ! -------------------------------------
@@ -158,4 +178,5 @@
 
          ! zgcr = matrix . gcr
+!$OMP PARALLEL DO
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -168,8 +189,30 @@
          ! rnorme = (gcr,gcr)
          rr = rnorme
-
-         ! zgcad = (zgcr,gcr) 
-         zsum(1) = glob_sum(gcr(:,:) * gcdmat(:,:) * gcr(:,:))
-         zsum(2) = glob_sum(gcr(:,:) * gcdmat(:,:) * zgcr(:,:) * bmask(:,:))
+         
+         ! zgcad = (zgcr,gcr)
+      tmp2 = 0.
+!$OMP PARALLEL
+!$OMP DO 
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            tmp2(ji, jj) = gcr(ji, jj) * gcdmat(ji, jj)
+            tmp1(ji, jj) = tmp2(ji, jj) * gcr(ji, jj)
+         END DO
+      END DO      
+!$OMP END DO
+!$OMP DO 
+!DIR$ IVDEP
+      DO jj = 1, jpj
+!DIR$ IVDEP
+         DO ji = 1, jpi
+            tmp2(ji, jj) = tmp2(ji, jj) * zgcr(ji, jj) * bmask(ji, jj)
+         END DO
+      END DO      
+!$OMP END DO
+!$OMP END PARALLEL
+ 
+!        zsum(1) = glob_sum(gcr(:,:) * gcdmat(:,:) * gcr(:,:))
+!        zsum(2) = glob_sum(gcr(:,:) * gcdmat(:,:) * zgcr(:,:) * bmask(:,:))
+         zsum = glob_asum_2d(tmp1, tmp2)
 
          !!RB we should gather the 2 glob_sum
@@ -190,4 +233,5 @@
          ! gcx = gcx + alph * gcdes
          ! gcr = gcr - alph * gccd
+!$OMP PARALLEL DO
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 9176)
@@ -231,4 +231,5 @@
       CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
+!$OMP PARALLEL DO PRIVATE(zh, zt, zs, ztm, zn3, zn2, zn1, zn0, zn)
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -271,4 +272,5 @@
       CASE( 1 )                !==  simplified EOS  ==!
          !
+!$OMP PARALLEL DO PRIVATE(zt, zs, zh, ztm, zn)
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -393,4 +395,5 @@
          ! Non-stochastic equation of state
          ELSE
+!$OMP PARALLEL DO PRIVATE(zh, zt, zs, ztm, zn3, zn2, zn1, zn0, zn0, zn)
             DO jk = 1, jpkm1
                DO jj = 1, jpj
@@ -435,4 +438,5 @@
       CASE( 1 )                !==  simplified EOS  ==!
          !
+!$OMP PARALLEL DO PRIVATE(zt, zs, zh, ztm, zn)
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -493,4 +497,5 @@
       CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
+!$OMP PARALLEL DO PRIVATE(zh, zt, zs, zn3, zn2, zn1, zn0, zn)
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
@@ -532,4 +537,5 @@
       CASE( 1 )                !==  simplified EOS  ==!
          !
+!$OMP PARALLEL DO PRIVATE(zt, zs, zh, zn)
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
@@ -583,4 +589,5 @@
       CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
+!$OMP PARALLEL DO PRIVATE(zh, zt, zs, ztm, zn3, zn2, zn1, zn0, zn)
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -640,4 +647,5 @@
       CASE( 1 )                  !==  simplified EOS  ==!
          !
+!$OMP PARALLEL DO PRIVATE(zt, zs, zh, ztm, zn)
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -697,4 +705,5 @@
       CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
+!$OMP PARALLEL DO PRIVATE(zh, zt, zs, zn3, zn2, zn1, zn0, zn)
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
@@ -755,4 +764,5 @@
       CASE( 1 )                  !==  simplified EOS  ==!
          !
+!$OMP PARALLEL DO PRIVATE(zt, zs, zh, zn)
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
@@ -910,4 +920,5 @@
       IF( nn_timing == 1 ) CALL timing_start('bn2')
       !
+!$OMP PARALLEL DO PRIVATE(zrw, zaw, zbw)
       DO jk = 2, jpkm1           ! interior points only (2=< jk =< jpkm1 )
          DO jj = 1, jpj          ! surface and bottom value set to zero one for all in istate.F90
@@ -962,4 +973,5 @@
       z1_T0   = 1._wp/40._wp
       !
+!$OMP PARALLEL DO PRIVATE(zt, zs,ztm, zn, zd)
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -1016,4 +1028,5 @@
       CASE ( -1, 1 )                !==  CT,SA (TEOS-10 formulation) ==!
          !
+!$OMP PARALLEL DO PRIVATE(zs)
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -1023,5 +1036,10 @@
             END DO
          END DO
-         ptf(:,:) = ptf(:,:) * psal(:,:)
+!$OMP PARALLEL DO
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ptf(ji,jj) = ptf(ji,jj) * psal(ji,jj)
+            END DO
+         END DO
          !
          IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:)
@@ -1029,6 +1047,11 @@
       CASE ( 0 )                     !==  PT,SP (UNESCO formulation)  ==!
          !
-         ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   &
-            &                     - 2.154996e-4_wp *       psal(:,:)   ) * psal(:,:)
+!$OMP PARALLEL DO
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ptf(ji,jj) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(ji,jj) )   &
+            &                     - 2.154996e-4_wp *       psal(ji,jj)   ) * psal(ji,jj)
+            END DO
+         END DO
             !
          IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:)
@@ -1125,4 +1148,5 @@
       CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
+!$OMP PARALLEL DO PRIVATE(zh, zt, zs, ztm, zn2, zn1, zn0, zn)
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -1188,4 +1212,5 @@
       CASE( 1 )                !==  Vallis (2006) simplified EOS  ==!
          !
+!$OMP PARALLEL DO PRIVATE(zt, zs, zh, ztm, zn)
          DO jk = 1, jpkm1
             DO jj = 1, jpj
Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 9176)
@@ -120,4 +120,5 @@
          ! --------------------------------------------------------------------
          ! upstream tracer flux in the i and j direction
+!$OMP PARALLEL DO PRIVATE(zfp_ui, zfm_ui, zfp_vj, zfm_vj)
          DO jk = 1, jpkm1
             DO jj = 1, jpjm1
@@ -133,7 +134,8 @@
             END DO
          END DO
-
+!$OMP END PARALLEL DO 
          ! upstream tracer flux in the k direction
          ! Interior value
+!$OMP PARALLEL DO PRIVATE(zfp_wk, zfm_wk)
          DO jk = 2, jpkm1
             DO jj = 1, jpj
@@ -145,4 +147,5 @@
             END DO
          END DO
+!$OMP END PARALLEL DO
          ! Surface value
          IF( lk_vvl ) THEN   
@@ -158,9 +161,11 @@
          ELSE                
             IF ( ln_isfcav ) THEN
+!$OMP PARALLEL DO 
                DO jj = 1, jpj
                   DO ji = 1, jpi
                      zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface 
                   END DO
-               END DO   
+               END DO
+!$OMP END PARALLEL DO   
             ELSE
                zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface
@@ -169,4 +174,5 @@
 
          ! total advective trend
+!$OMP PARALLEL DO PRIVATE(z2dtt, ztra)
          DO jk = 1, jpkm1
             z2dtt = p2dt(jk)
@@ -183,4 +189,5 @@
             END DO
          END DO
+!$OMP END PARALLEL DO
          !                             ! Lateral boundary conditions on zwi  (unchanged sign)
          CALL lbc_lnk( zwi, 'T', 1. )  
@@ -200,4 +207,5 @@
          ! --------------------------------------------------
          ! antidiffusive flux on i and j
+!$OMP PARALLEL DO
          DO jk = 1, jpkm1
             DO jj = 1, jpjm1
@@ -208,7 +216,8 @@
             END DO
          END DO
-      
+!$OMP END PARALLEL DO
          ! antidiffusive flux on k
          ! Interior value
+!$OMP PARALLEL DO
          DO jk = 2, jpkm1                    
             DO jj = 1, jpj
@@ -218,4 +227,5 @@
             END DO
          END DO
+!$OMP END PARALLEL DO
          ! surface value
          IF ( ln_isfcav ) THEN
@@ -238,4 +248,5 @@
          ! 5. final trend with corrected fluxes
          ! ------------------------------------
+!$OMP PARALLEL DO PRIVATE(zbtr, ztra)
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -251,5 +262,5 @@
             END DO
          END DO
-
+!$OMP END PARALLEL DO
          !                                 ! trend diagnostics (contribution of upstream fluxes)
          IF( l_trd )  THEN 
@@ -356,4 +367,5 @@
          ! --------------------------------------------------------------------
          ! upstream tracer flux in the i and j direction
+!$OMP PARALLEL DO PRIVATE(zfp_ui, zfm_ui, zfp_vj, zfm_vj)
          DO jk = 1, jpkm1
             DO jj = 1, jpjm1
@@ -369,7 +381,8 @@
             END DO
          END DO
-
+!$OMP END PARALLEL DO
          ! upstream tracer flux in the k direction
          ! Interior value
+!$OMP PARALLEL DO PRIVATE(zfp_wk, zfm_wk)
          DO jk = 2, jpkm1
             DO jj = 1, jpj
@@ -381,4 +394,5 @@
             END DO
          END DO
+!$OMP END PARALLEL DO
          ! Surface value
          IF( lk_vvl ) THEN
@@ -394,4 +408,5 @@
          ELSE
             IF ( ln_isfcav ) THEN
+!$OMP PARALLEL DO 
                DO jj = 1, jpj
                   DO ji = 1, jpi
@@ -399,4 +414,5 @@
                   END DO
                END DO
+!$OMP END PARALLEL DO
             ELSE
                zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)                                               ! linear free surface + no isf
@@ -405,4 +421,5 @@
 
          ! total advective trend
+!$OMP PARALLEL DO PRIVATE(z2dtt, ztra)
          DO jk = 1, jpkm1
             z2dtt = p2dt(jk)
@@ -419,4 +436,5 @@
             END DO
          END DO
+!$OMP END PARALLEL DO
          !                             ! Lateral boundary conditions on zwi  (unchanged sign)
          CALL lbc_lnk( zwi, 'T', 1. )  
@@ -437,4 +455,5 @@
          ! antidiffusive flux on i and j
          !
+!$OMP PARALLEL DO 
          DO jk = 1, jpkm1
             !
@@ -463,5 +482,5 @@
             END DO
          END DO
-      
+!$OMP END PARALLEL DO      
          ! antidiffusive flux on k
          zwz(:,:,1) = 0._wp        ! Surface value
@@ -489,4 +508,5 @@
               jta = MOD(jta,3) + 1
             ENDIF
+!$OMP PARALLEL DO
             DO jk = 2, jpkm1          ! Interior value
                DO jj = 2, jpjm1
@@ -497,7 +517,7 @@
                END DO
             END DO
-
+!$OMP END PARALLEL DO
             jtaken = MOD( jtaken + 1 , 2 )
-
+!$OMP PARALLEL DO PRIVATE (zbtr, ztra)
             DO jk = 2, jpkm1          ! Interior value
                DO jj = 2, jpjm1
@@ -510,7 +530,7 @@
                END DO
             END DO
-
-         END DO
-
+!$OMP END PARALLEL DO
+         END DO
+!$OMP PARALLEL DO 
          DO jk = 2, jpkm1          ! Anti-diffusive vertical flux using average flux from the sub-timestepping
             DO jj = 2, jpjm1
@@ -520,4 +540,5 @@
             END DO
          END DO
+!$OMP END PARALLEL DO
          CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions
          CALL lbc_lnk( zwz, 'W',  1. )
@@ -530,4 +551,5 @@
          ! 5. final trend with corrected fluxes
          ! ------------------------------------
+!$OMP PARALLEL DO PRIVATE(zbtr, ztra)
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -543,5 +565,5 @@
             END DO
          END DO
-
+!$OMP END PARALLEL DO
          !                                 ! trend diagnostics (contribution of upstream fluxes)
          IF( l_trd )  THEN 
@@ -612,4 +634,5 @@
          &        paft * tmask + zbig * ( 1._wp - tmask )  )
 
+!$OMP PARALLEL DO PRIVATE(ikm1, z2dtt, zup, zdo, zpos, zneg, zbt)
       DO jk = 1, jpkm1
          ikm1 = MAX(jk-1,1)
@@ -647,8 +670,10 @@
          END DO
       END DO
+!$OMP END PARALLEL DO
       CALL lbc_lnk( zbetup, 'T', 1. )   ;   CALL lbc_lnk( zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
 
       ! 3. monotonic flux in the i & j direction (paa & pbb)
       ! ----------------------------------------
+!$OMP PARALLEL DO PRIVATE(zau, zbu, zcu, zav, zbv, zcv, za, zb, zc)
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1
@@ -673,4 +698,5 @@
          END DO
       END DO
+!$OMP END PARALLEL DO
       CALL lbc_lnk( paa, 'U', -1. )   ;   CALL lbc_lnk( pbb, 'V', -1. )   ! lateral boundary condition (changed sign)
       !
Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 9176)
@@ -89,4 +89,5 @@
       !
       !                             !  Add the geothermal heat flux trend on temperature
+!$OMP PARALLEL DO PRIVATE(ik, zqgh_trd)
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 9176)
@@ -108,5 +108,5 @@
       REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      -
       REAL(wp) ::  zcoef0, zbtr, ztra            !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
+      REAL(wp), DIMENSION(jpi,jpj  ) ::  z2d
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdkt, zdk1t, zdit, zdjt, ztfw 
       !!----------------------------------------------------------------------
@@ -114,5 +114,5 @@
       IF( nn_timing == 1 )  CALL timing_start('tra_ldf_iso')
       !
-      CALL wrk_alloc( jpi, jpj,      z2d ) 
+!     CALL wrk_alloc( jpi, jpj,      z2d ) 
       CALL wrk_alloc( jpi, jpj, jpk, zdit, zdjt, ztfw, zdkt, zdk1t ) 
       !
@@ -137,4 +137,5 @@
 
          ! Horizontal tracer gradient 
+!$OMP PARALLEL DO
          DO jk = 1, jpkm1
             DO jj = 1, jpjm1
@@ -145,7 +146,9 @@
             END DO
          END DO
+!$OMP END PARALLEL DO
 
          ! partial cell correction
          IF( ln_zps ) THEN      ! partial steps correction at the last ocean level 
+!$OMP PARALLEL DO
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
@@ -157,4 +160,5 @@
          ENDIF
          IF( ln_zps .AND. ln_isfcav ) THEN      ! partial steps correction at the first wet level beneath a cavity
+!$OMP PARALLEL DO
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
@@ -173,4 +177,5 @@
          ! 
          ! interior value 
+!$OMP PARALLEL DO
          DO jk = 2, jpkm1               
             DO jj = 1, jpj
@@ -182,8 +187,10 @@
             END DO
          END DO
+!$OMP END PARALLEL DO
          ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
          zdk1t(:,:,1) = ( ptb(:,:,1,jn  ) - ptb(:,:,2,jn) ) * wmask(:,:,2)
          zdkt (:,:,1) = zdk1t(:,:,1)
          IF ( ln_isfcav ) THEN
+!$OMP PARALLEL DO PRIVATE(ikt)
             DO jj = 1, jpj
                DO ji = 1, jpi   ! vector opt.
@@ -193,8 +200,10 @@
                END DO
             END DO
+!$OMP END PARALLEL DO
          END IF
 
          ! 2. Horizontal fluxes
          ! --------------------   
+!$OMP PARALLEL DO PRIVATE(zabe1, zabe2, zmsku, zmskv, zcof1, zcof2, zbtr, ztra)
          DO jk = 1, jpkm1
             DO jj = 1 , jpjm1
@@ -233,4 +242,5 @@
          END DO                                        !   End of slab  
          !                                             ! ===============
+!$OMP END PARALLEL DO
          !
          ! "Poleward" diffusive heat or salt transports (T-S case only)
@@ -245,4 +255,5 @@
            IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN
                z2d(:,:) = 0._wp 
+!$OMP PARALLEL DO REDUCTION(+:z2d)
                DO jk = 1, jpkm1
                   DO jj = 2, jpjm1
@@ -252,9 +263,12 @@
                   END DO
                END DO
+
                z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043)
+
                CALL lbc_lnk( z2d, 'U', -1. )
                CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction
                !
                z2d(:,:) = 0._wp 
+!$OMP PARALLEL DO REDUCTION(+:z2d)
                DO jk = 1, jpkm1
                   DO jj = 2, jpjm1
@@ -264,4 +278,5 @@
                   END DO
                END DO
+
                z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043)
                CALL lbc_lnk( z2d, 'V', -1. )
@@ -286,4 +301,5 @@
          
          ! interior (2=<jk=<jpk-1)
+!$OMP PARALLEL DO PRIVATE(zcoef0, zmsku, zmskv, zcoef3, zcoef4 )
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
@@ -306,8 +322,9 @@
             END DO
          END DO
-         
+!$OMP END PARALLEL DO         
          
          ! I.5 Divergence of vertical fluxes added to the general tracer trend
          ! -------------------------------------------------------------------
+!$OMP PARALLEL DO PRIVATE(zbtr, ztra)
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -319,8 +336,9 @@
             END DO
          END DO
+!$OMP END PARALLEL DO
          !
       END DO
       !
-      CALL wrk_dealloc( jpi, jpj, z2d ) 
+!     CALL wrk_dealloc( jpi, jpj, z2d ) 
       CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt, ztfw, zdkt, zdk1t ) 
       !
Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 9176)
@@ -153,4 +153,5 @@
      ! trends computation
       IF( l_trdtra ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
+!$OMP PARALLEL DO PRIVATE(zfact)
          DO jk = 1, jpkm1
             zfact = 1._wp / r2dtra(jk)             
@@ -304,4 +305,6 @@
       !
       DO jn = 1, kjpt      
+!$OMP PARALLEL DO PRIVATE(zfact1, zfact2, ze3t_b, ze3t_n, ze3t_a, ztc_b, ztc_n, ztc_a, ze3t_d, ztc_d, &
+!$OMP&                    ze3t_f, ztc_f, ze3t_d)
          DO jk = 1, jpkm1
             zfact1 = atfp * p2dt(jk)
Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 9176)
@@ -129,4 +129,5 @@
             ! isoneutral diffusion: add the contribution 
             IF( ln_traldf_grif    ) THEN     ! Griffies isoneutral diff
+!$OMP PARALLEL DO
                DO jk = 2, jpkm1
                   DO jj = 2, jpjm1
@@ -137,4 +138,5 @@
                END DO
             ELSE IF( l_traldf_rot ) THEN     ! standard isoneutral diff
+!$OMP PARALLEL DO
                DO jk = 2, jpkm1
                   DO jj = 2, jpjm1
@@ -149,4 +151,5 @@
 #endif
             ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked)
+!$OMP PARALLEL DO PRIVATE(ze3ta, ze3tn)
             DO jk = 1, jpkm1
                DO jj = 2, jpjm1
@@ -187,5 +190,8 @@
                END DO
             END DO
+
+!$OMP PARALLEL
             DO jk = 2, jpkm1
+!$OMP DO
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1
@@ -193,9 +199,12 @@
                   END DO
                END DO
-            END DO
+!$OMP END DO
+            END DO
+!$OMP END PARALLEL
             !
          END IF 
          !         
          ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
+!$OMP PARALLEL DO PRIVATE(ze3tb, ze3tn)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
@@ -206,5 +215,8 @@
             END DO
          END DO
+
+!$OMP PARALLEL
          DO jk = 2, jpkm1
+!$OMP DO PRIVATE(ze3tb, ze3tn, zrhs)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
@@ -215,7 +227,10 @@
                END DO
             END DO
-         END DO
-
+!$OMP END DO
+         END DO
+!$OMP END PARALLEL
          ! third recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer)
+!$OMP PARALLEL 
+!$OMP DO 
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
@@ -223,5 +238,7 @@
             END DO
          END DO
+
          DO jk = jpk-2, 1, -1
+!$OMP DO
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
@@ -230,5 +247,7 @@
                END DO
             END DO
-         END DO
+!$OMP END DO
+         END DO
+!$OMP END PARALLEL
          !                                            ! ================= !
       END DO                                          !  end tracer loop  !
Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 9176)
@@ -127,4 +127,5 @@
             zwt(:,:, 1 ) = 0._wp   ;   zws(:,:, 1 ) = 0._wp            ! vertical diffusive fluxes
             zwt(:,:,jpk) = 0._wp   ;   zws(:,:,jpk) = 0._wp
+!$OMP PARALLEL DO
             DO jk = 2, jpk
                zwt(:,:,jk) =   avt(:,:,jk) * ( tsa(:,:,jk-1,jp_tem) - tsa(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk)
@@ -133,4 +134,5 @@
             !
             ztrdt(:,:,jpk) = 0._wp   ;   ztrds(:,:,jpk) = 0._wp
+!$OMP PARALLEL DO
             DO jk = 1, jpkm1
                ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / fse3t(:,:,jk)
@@ -142,5 +144,8 @@
             !
          CASE DEFAULT                 ! other trends: mask and send T & S trends to trd_tra_mng
-            ztrds(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
+!$OMP PARALLEL DO
+           DO jk = 1, jpk
+              ztrds(:,:,jk) = ptrd(:,:,jk) * tmask(:,:,jk)
+           ENDDO
             CALL trd_tra_mng( trdt, ztrds, ktrd, kt )  
          END SELECT
@@ -200,4 +205,5 @@
       ptrd(:,:,jpk) = 0._wp
       !
+!$OMP PARALLEL DO
       DO jk = 1, jpkm1         ! advective trend
          DO jj = 2, jpjm1
Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90	(revision 9176)
@@ -77,4 +77,5 @@
          zavm_evd(:,:,:) = avm(:,:,:)           ! set avm prior to evd application
          !
+!$OMP PARALLEL DO
          DO jk = 1, jpkm1 
             DO jj = 2, jpj             ! no vector opt.
@@ -103,4 +104,5 @@
          !
       CASE DEFAULT         ! enhance vertical eddy diffusivity only (if rn2<-1.e-12) 
+!$OMP PARALLEL DO
          DO jk = 1, jpkm1
 !!!         WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd   ! agissant sur T SEUL! 
Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 9176)
@@ -227,6 +227,6 @@
       REAL(wp) ::   zbbrau, zesh2                   ! temporary scalars
       REAL(wp) ::   zfact1, zfact2, zfact3          !    -         -
-      REAL(wp) ::   ztx2  , zty2  , zcof            !    -         -
-      REAL(wp) ::   ztau  , zdif                    !    -         -
+      REAL(wp) ::   ztx2  , zty2  , zcof, zcofa     !    -         -
+      REAL(wp) ::   ztau  , zdif, zdifa             !    -         -
       REAL(wp) ::   zus   , zwlc  , zind            !    -         -
       REAL(wp) ::   zzd_up, zzd_lw                  !    -         -
@@ -253,4 +253,5 @@
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF ( ln_isfcav ) THEN
+!$OMP PARALLEL DO
          DO jj = 2, jpjm1            ! en(mikt(ji,jj))   = rn_emin
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -259,4 +260,5 @@
          END DO
       END IF
+!$OMP PARALLEL DO
       DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -296,11 +298,20 @@
          !                        !* total energy produce by LC : cumulative sum over jk
          zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * fsdepw(:,:,1) * fse3w(:,:,1)
+!$OMP PARALLEL
          DO jk = 2, jpk
-            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
-         END DO
+!$OMP DO
+            DO jj = 1, jpj
+                  zpelc(:,jj,jk)  = zpelc(:,jj,jk-1) + MAX( rn2b(:,jj,jk), 0._wp ) * fsdepw(:,jj,jk) * fse3w(:,jj,jk)
+            END DO
+!$OMP END DO
+         END DO
+!$OMP END PARALLEL
          !                        !* finite Langmuir Circulation depth
          zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
+         zcofa = 0.016 / SQRT( zrhoa * zcdrag )
          imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
+!$OMP PARALLEL SHARED(imlc)
          DO jk = jpkm1, 2, -1
+!$OMP DO PRIVATE(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,6 +320,8 @@
                END DO
             END DO
+!$OMP END DO
          END DO
          !                               ! finite LC depth
+!$OMP DO
          DO jj = 1, jpj 
             DO ji = 1, jpi
@@ -316,12 +329,11 @@
             END DO
          END DO
-         zcof = 0.016 / SQRT( zrhoa * zcdrag )
-!CDIR NOVERRCHK
+!$OMP END DO
+!         zcof = 0.016 / SQRT( zrhoa * zcdrag )
+!$OMP DO PRIVATE(zus, zind, zwlc)
          DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en
-!CDIR NOVERRCHK
-            DO jj = 2, jpjm1
-!CDIR NOVERRCHK
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zus  = zcofa * SQRT( taum(ji,jj) )           ! Stokes drift
                   !                                           ! vertical velocity due to LC
                   zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) )
@@ -333,4 +345,6 @@
             END DO
          END DO
+!$OMP END DO
+!$OMP END PARALLEL
          !
       ENDIF
@@ -343,4 +357,5 @@
       !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
       !
+!$OMP PARALLEL DO
       DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form)
          DO jj = 1, jpj                 ! here avmu, avmv used as workspace
@@ -358,4 +373,5 @@
       END DO
       !
+!$OMP PARALLEL DO PRIVATE(zcof, zzd_up, zzd_lw, zesh2)
       DO jk = 2, jpkm1           !* Matrix and right hand side in en
          DO jj = 2, jpjm1
@@ -390,5 +406,7 @@
       END DO
       !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
+!$OMP PARALLEL
       DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
+!$OMP DO 
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
@@ -396,7 +414,10 @@
             END DO
          END DO
-      END DO
+!$OMP END DO
+      END DO
+!$OMP END PARALLEL
       !
       ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+!$OMP PARALLEL DO
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -404,5 +425,7 @@
          END DO
       END DO
+!$OMP PARALLEL
       DO jk = 3, jpkm1
+!$OMP DO
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
@@ -410,7 +433,10 @@
             END DO
          END DO
-      END DO
+!$OMP END DO
+      END DO
+!$OMP END PARALLEL
       !
       ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+!$OMP PARALLEL DO
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -418,5 +444,7 @@
          END DO
       END DO
+!$OMP PARALLEL
       DO jk = jpk-2, 2, -1
+!$OMP DO
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
@@ -424,5 +452,8 @@
             END DO
          END DO
-      END DO
+!$OMP END DO
+      END DO
+!$OMP END PARALLEL
+!$OMP PARALLEL DO
       DO jk = 2, jpkm1                             ! set the minimum value of tke
          DO jj = 2, jpjm1
@@ -440,4 +471,5 @@
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF( nn_htau == 2 ) THEN           !* mixed-layer depth dependant length scale
+!$OMP PARALLEL DO
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -452,4 +484,5 @@
       !
       IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
+!$OMP PARALLEL DO
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
@@ -461,4 +494,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 PRIVATE(jk)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -469,15 +503,13 @@
          END DO
       ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
-!CDIR NOVERRCHK
+!$OMP PARALLEL DO PRIVATE(ztx2, zty2, ztau, zdif, zdifa)
          DO jk = 2, jpkm1
-!CDIR NOVERRCHK
-            DO jj = 2, jpjm1
-!CDIR NOVERRCHK
+            DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
                   zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
                   ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress 
-                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean 
-                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
+                  zdifa = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean 
+                  zdif = rhftau_scl * MAX( 0._wp, zdifa + rhftau_add )  ! apply some modifications...
                   en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
                      &                        * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
@@ -487,4 +519,5 @@
       ELSEIF( nn_etau == 4 ) THEN       !* column integral independant of htau (rn_efr must be scaled up)
          IF( nn_htau == 2 ) THEN        ! efr dependant on time-varying htau 
+!$OMP PARALLEL DO
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -493,4 +526,5 @@
             END DO
          ENDIF
+!$OMP PARALLEL DO
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
@@ -504,4 +538,5 @@
       CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
       !
+!$OMP PARALLEL DO
       DO jk = 2, jpkm1                             ! TKE budget: near-inertial waves term  
          DO jj = 2, jpjm1  
@@ -580,4 +615,5 @@
       !
       IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g)
+!$OMP PARALLEL DO PRIVATE(zraug)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
@@ -590,9 +626,7 @@
       ENDIF
       !
-!CDIR NOVERRCHK
+!$OMP PARALLEL DO PRIVATE(zrn2)
       DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2)
-!CDIR NOVERRCHK
          DO jj = 2, jpjm1
-!CDIR NOVERRCHK
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zrn2 = MAX( rn2(ji,jj,jk), rsmall )
@@ -611,4 +645,5 @@
       ! where wmask = 0 set zmxlm == fse3w
       CASE ( 0 )           ! bounded by the distance to surface and bottom
+!$OMP PARALLEL DO PRIVATE(zemxl)
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
@@ -624,4 +659,5 @@
          !
       CASE ( 1 )           ! bounded by the vertical scale factor
+!$OMP PARALLEL DO PRIVATE(zemxl)
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
@@ -635,5 +671,7 @@
          !
       CASE ( 2 )           ! |dk[xml]| bounded by e3t :
+!$OMP PARALLEL
          DO jk = 2, jpkm1         ! from the surface to the bottom :
+!$OMP DO
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -643,4 +681,5 @@
          END DO
          DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
+!$OMP DO PRIVATE(zemxl)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -651,7 +690,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
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -661,4 +703,5 @@
          END DO
          DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
+!$OMP DO
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -667,9 +710,7 @@
             END DO
          END DO
-!CDIR NOVERRCHK
+!$OMP DO PRIVATE(zemlm, zemlp)
          DO jk = 2, jpkm1
-!CDIR NOVERRCHK
-            DO jj = 2, jpjm1
-!CDIR NOVERRCHK
+            DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
@@ -680,4 +721,5 @@
             END DO
          END DO
+!$OMP END PARALLEL
          !
       END SELECT
@@ -691,9 +733,6 @@
       !                     !  Vertical eddy viscosity and diffusivity  (avmu, avmv, avt)
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-!CDIR NOVERRCHK
       DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points
-!CDIR NOVERRCHK
          DO jj = 2, jpjm1
-!CDIR NOVERRCHK
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zsqen = SQRT( en(ji,jj,jk) )
@@ -894,4 +933,5 @@
       ENDIF
       !                               !* set vertical eddy coef. to the background value
+!$OMP PARALLEL DO
       DO jk = 1, jpk
          avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
@@ -959,4 +999,5 @@
         ELSE                                   !* Start from rest
            en(:,:,:) = rn_emin * tmask(:,:,:)
+!$OMP PARALLEL DO
            DO jk = 1, jpk                           ! set the Kz to the background value
               avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90	(revision 9176)
@@ -111,10 +111,10 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   ztpc         ! scalar workspace
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zkz
+      REAL(wp), DIMENSION(jpi,jpj) ::   zkz
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('zdf_tmx')
       !
-      CALL wrk_alloc( jpi,jpj, zkz )
+!     CALL wrk_alloc( jpi,jpj, zkz )
 
       !                          ! ----------------------- !
@@ -128,5 +128,4 @@
          zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk) * wmask(:,:,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 ji = 1, jpi
@@ -134,5 +133,4 @@
          END DO
       END DO
-
       DO jk = 2, jpkm1     !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s
          DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
@@ -142,5 +140,4 @@
          END DO
       END DO
-
       IF( kt == nit000 ) THEN       !* check at first time-step: diagnose the energy consumed by zav_tide
          ztpc = 0.e0
@@ -166,4 +163,5 @@
       !                          !   Update  mixing coefs  !                          
       !                          ! ----------------------- !
+!$OMP PARALLEL DO
       DO jk = 2, jpkm1              !* update momentum & tracer diffusivity with tidal mixing
          DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
@@ -174,5 +172,5 @@
          END DO
       END DO
-      
+!$OMP PARALLEL DO
       DO jk = 2, jpkm1              !* update momentum & tracer diffusivity with tidal mixing
          DO jj = 2, jpjm1
@@ -190,5 +188,5 @@
       IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_tide , clinfo1=' tmx - av_tide: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk)
       !
-      CALL wrk_dealloc( jpi,jpj, zkz )
+!     CALL wrk_dealloc( jpi,jpj, zkz )
       !
       IF( nn_timing == 1 )  CALL timing_stop('zdf_tmx')
@@ -222,6 +220,6 @@
       INTEGER  ::   ji, jj, jk    ! dummy loop indices
       REAL(wp) ::   zcoef, ztpc   ! temporary scalar
-      REAL(wp), DIMENSION(:,:)  , POINTER ::   zkz                        ! 2D workspace
-      REAL(wp), DIMENSION(:,:)  , POINTER ::   zsum1 , zsum2 , zsum       !  -      -
+      REAL(wp), DIMENSION(jpi, jpj)  ::   zkz                        ! 2D workspace
+      REAL(wp), DIMENSION(jpi, jpj)  ::   zsum1 , zsum2 , zsum       !  -      -
       REAL(wp), DIMENSION(:,:,:), POINTER ::   zempba_3d_1, zempba_3d_2   ! 3D workspace
       REAL(wp), DIMENSION(:,:,:), POINTER ::   zempba_3d  , zdn2dz        !  -      -
@@ -231,5 +229,5 @@
       IF( nn_timing == 1 )  CALL timing_start('tmx_itf')
       !
-      CALL wrk_alloc( jpi,jpj, zkz, zsum1 , zsum2 , zsum )
+!     CALL wrk_alloc( jpi,jpj, zkz, zsum1 , zsum2 , zsum )
       CALL wrk_alloc( jpi,jpj,jpk, zempba_3d_1, zempba_3d_2, zempba_3d, zdn2dz, zavt_itf )
 
@@ -237,7 +235,7 @@
       zempba_3d_1(:,:,jpk) = 0.e0
       zempba_3d_2(:,:,jpk) = 0.e0
+!$OMP PARALLEL DO
       DO jk = 1, jpkm1             
          zdn2dz     (:,:,jk) = rn2(:,:,jk) - rn2(:,:,jk+1)           ! Vertical profile of dN2/dz
-!CDIR NOVERRCHK
          zempba_3d_1(:,:,jk) = SQRT(  MAX( 0.e0, rn2(:,:,jk) )  )    !    -        -    of N
          zempba_3d_2(:,:,jk) =        MAX( 0.e0, rn2(:,:,jk) )       !    -        -    of N^2
@@ -257,5 +255,4 @@
          END DO
       END DO
-
       DO jk= 1, jpk
          DO jj = 1, jpj
@@ -313,4 +310,5 @@
 
       !                             ! Update pav with the ITF mixing coefficient
+!$OMP PARALLEL DO
       DO jk = 2, jpkm1
          pav(:,:,jk) = pav     (:,:,jk) * ( 1.e0 - mask_itf(:,:) )   &
@@ -318,5 +316,5 @@
       END DO
       !
-      CALL wrk_dealloc( jpi,jpj, zkz, zsum1 , zsum2 , zsum )
+!     CALL wrk_dealloc( jpi,jpj, zkz, zsum1 , zsum2 , zsum )
       CALL wrk_dealloc( jpi,jpj,jpk, zempba_3d_1, zempba_3d_2, zempba_3d, zdn2dz, zavt_itf )
       !
Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90	(revision 9176)
@@ -27,4 +27,5 @@
    PUBLIC   DDPDD      ! also used in closea module
    PUBLIC   glob_min, glob_max
+   PUBLIC   glob_asum_2d
 #if defined key_nosignedzero
    PUBLIC SIGN
@@ -189,20 +190,23 @@
    END FUNCTION glob_sum_1d
 
-   FUNCTION glob_sum_2d( ptab )
-      !!----------------------------------------------------------------------
-      !!                  ***  FUNCTION  glob_sum_2d ***
+   FUNCTION sum_2d_ref( ptab )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION  sum_2d_ref ***
       !!
       !! ** Purpose : perform a sum in calling DDPDD routine
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab
-      REAL(wp)                             ::   glob_sum_2d   ! global masked sum
+      COMPLEX(wp)                          ::   sum_2d_ref   ! global masked sum
       !!
       COMPLEX(wp)::   ctmp
       REAL(wp)   ::   ztmp
+!$    COMPLEX(wp)::   comp
       INTEGER    ::   ji, jj   ! dummy loop indices
       !!-----------------------------------------------------------------------
       !
-      ztmp = 0.e0
       ctmp = CMPLX( 0.e0, 0.e0, wp )
+!$    comp = CMPLX( 0.e0, 0.e0, wp )
+!$OMP PARALLEL FIRSTPRIVATE(ctmp) PRIVATE(ztmp) SHARED(comp)
+!$OMP DO
       DO jj = 1, jpj
          DO ji =1, jpi
@@ -211,4 +215,29 @@
          END DO
       END DO
+!$OMP ENDDO
+!$OMP CRITICAL
+!$    CALL DDPDD( ctmp, comp ) 
+!$OMP END CRITICAL
+!$OMP END PARALLEL
+!$    ctmp = comp
+      sum_2d_ref = ctmp
+      !
+   END FUNCTION sum_2d_ref
+
+   FUNCTION glob_sum_2d( ptab )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_sum_2d ***
+      !!
+      !! ** Purpose : perform a sum in calling DDPDD routine
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab
+      REAL(wp)                             ::   glob_sum_2d   ! global masked sum
+      !!
+      COMPLEX(wp)::   ctmp
+      REAL(wp)   ::   ztmp
+      INTEGER    ::   ji, jj   ! dummy loop indices
+      !!-----------------------------------------------------------------------
+      !
+      ctmp = sum_2d_ref(ptab)
       IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain
       glob_sum_2d = REAL(ctmp,wp)
@@ -228,4 +257,5 @@
       COMPLEX(wp)::   ctmp
       REAL(wp)   ::   ztmp
+!$    COMPLEX(wp)::   comp
       INTEGER    ::   ji, jj, jk   ! dummy loop indices
       INTEGER    ::   ijpk ! local variables: size of ptab
@@ -234,6 +264,8 @@
       ijpk = SIZE(ptab,3)
       !
-      ztmp = 0.e0
       ctmp = CMPLX( 0.e0, 0.e0, wp )
+!$    comp = CMPLX( 0.e0, 0.e0, wp )
+!$OMP PARALLEL FIRSTPRIVATE(ctmp) PRIVATE(ztmp) SHARED(comp)
+!$OMP DO
       DO jk = 1, ijpk
          DO jj = 1, jpj
@@ -244,4 +276,10 @@
          END DO
       END DO
+!$OMP ENDDO
+!$OMP CRITICAL
+!$    CALL DDPDD( ctmp, comp ) 
+!$OMP END CRITICAL
+!$OMP END PARALLEL
+!$    ctmp = comp
       IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain
       glob_sum_3d = REAL(ctmp,wp)
@@ -261,4 +299,5 @@
       COMPLEX(wp)::   ctmp
       REAL(wp)   ::   ztmp
+!$    COMPLEX(wp)::   comp
       INTEGER    ::   ji, jj   ! dummy loop indices
       !!-----------------------------------------------------------------------
@@ -266,4 +305,7 @@
       ztmp = 0.e0
       ctmp = CMPLX( 0.e0, 0.e0, wp )
+!$    comp = CMPLX( 0.e0, 0.e0, wp )
+!$OMP PARALLEL FIRSTPRIVATE(ctmp) PRIVATE(ztmp) SHARED(comp)
+!$OMP DO
       DO jj = 1, jpj
          DO ji =1, jpi
@@ -274,4 +316,10 @@
          END DO
       END DO
+!$OMP ENDDO
+!$OMP CRITICAL
+!$    CALL DDPDD( ctmp, comp ) 
+!$OMP END CRITICAL
+!$OMP END PARALLEL
+!$    ctmp = comp
       IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain
       glob_sum_2d_a = REAL(ctmp,wp)
@@ -279,4 +327,22 @@
    END FUNCTION glob_sum_2d_a
 
+   FUNCTION glob_asum_2d( ptab1, ptab2 )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_sum_2d_a ***
+      !!
+      !! ** Purpose : perform a sum on two 2D arrays in calling DDPDD routine
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab1, ptab2
+      REAL(wp),             DIMENSION(2)   ::   glob_asum_2d   ! global masked sum
+      !!
+      COMPLEX(wp),          DIMENSION(2)   ::   ctmp
+      !!-----------------------------------------------------------------------
+      !
+      ctmp(1) =  sum_2d_ref(ptab1)
+      ctmp(2) =  sum_2d_ref(ptab2)
+      IF( lk_mpp )   CALL mpp_sum( ctmp, 2 )   ! sum over the global domain
+      glob_asum_2d = REAL(ctmp,wp)
+      !
+   END FUNCTION glob_asum_2d
 
    FUNCTION glob_sum_3d_a( ptab1, ptab2 )
@@ -291,4 +357,5 @@
       COMPLEX(wp)::   ctmp
       REAL(wp)   ::   ztmp
+!$    COMPLEX(wp)::   comp
       INTEGER    ::   ji, jj, jk   ! dummy loop indices
       INTEGER    ::   ijpk ! local variables: size of ptab
@@ -299,4 +366,7 @@
       ztmp = 0.e0
       ctmp = CMPLX( 0.e0, 0.e0, wp )
+!$    comp = CMPLX( 0.e0, 0.e0, wp )
+!$OMP PARALLEL FIRSTPRIVATE(ctmp) PRIVATE(ztmp) SHARED(comp)
+!$OMP DO
       DO jk = 1, ijpk
          DO jj = 1, jpj
@@ -309,4 +379,10 @@
          END DO    
       END DO
+!$OMP ENDDO
+!$OMP CRITICAL
+!$    CALL DDPDD( ctmp, comp ) 
+!$OMP END CRITICAL
+!$OMP END PARALLEL
+!$    ctmp = comp
       IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain
       glob_sum_3d_a = REAL(ctmp,wp)
@@ -317,4 +393,26 @@
 
    ! --- MIN ---
+   FUNCTION glob_min_2d_ref( ptab ) 
+      !!-----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_min_2D  ***
+      !!
+      !! ** Purpose : perform a masked min on the inner global domain of a 2D array
+      !!-----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::   ptab          ! input 2D array
+      REAL(wp)                             ::   glob_min_2d_ref   ! global masked min
+      INTEGER                              ::   jj, ji        ! local index
+      !!-----------------------------------------------------------------------
+      !
+      glob_min_2d_ref = 1.e32
+!$OMP PARALLEL DO REDUCTION(MIN:glob_min_2d_ref)
+      DO jj = 1, jpj
+         DO ji =1, jpi
+            glob_min_2d_ref = MIN(glob_min_2d_ref, ptab(ji,jj)*tmask_i(ji,jj) )
+         ENDDO
+      ENDDO
+!$OMP END PARALLEL DO
+      !
+   END FUNCTION glob_min_2d_ref
+
    FUNCTION glob_min_2d( ptab ) 
       !!-----------------------------------------------------------------------
@@ -323,12 +421,37 @@
       !! ** Purpose : perform a masked min on the inner global domain of a 2D array
       !!-----------------------------------------------------------------------
-      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab          ! input 2D array
+      REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::   ptab          ! input 2D array
       REAL(wp)                             ::   glob_min_2d   ! global masked min
-      !!-----------------------------------------------------------------------
-      !
-      glob_min_2d = MINVAL( ptab(:,:)*tmask_i(:,:) )
+      INTEGER                              ::   jj, ji        ! local index
+      !!-----------------------------------------------------------------------
+      !
+      glob_min_2d = glob_min_2d_ref (ptab)
       IF( lk_mpp )   CALL mpp_min( glob_min_2d )
       !
    END FUNCTION glob_min_2d
+
+   FUNCTION glob_min_3d_ref( ptab ) 
+      !!-----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_min_3D  ***
+      !!
+      !! ** Purpose : perform a masked min on the inner global domain of a 3D array
+      !!-----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab          ! input 3D array
+      REAL(wp)                               ::   glob_min_3d_ref   ! global masked min
+      !!
+      INTEGER :: jk
+      INTEGER :: ijpk ! local variable: size of the 3d dimension of ptab
+      !!-----------------------------------------------------------------------
+      !
+      ijpk = SIZE(ptab,3)
+      !
+      glob_min_3d_ref = 1.e32
+!$OMP PARALLEL DO REDUCTION(MIN:glob_min_3d_ref)
+      DO jk = 1, ijpk
+         glob_min_3d_ref = MIN( glob_min_3d_ref, MINVAL( ptab(:,:,jk)*tmask_i(:,:) ) )
+      END DO
+!$OMP END PARALLEL DO
+      !
+   END FUNCTION glob_min_3d_ref
  
    FUNCTION glob_min_3d( ptab ) 
@@ -340,15 +463,7 @@
       REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab          ! input 3D array
       REAL(wp)                               ::   glob_min_3d   ! global masked min
-      !!
-      INTEGER :: jk
-      INTEGER :: ijpk ! local variable: size of the 3d dimension of ptab
-      !!-----------------------------------------------------------------------
-      !
-      ijpk = SIZE(ptab,3)
-      !
-      glob_min_3d = MINVAL( ptab(:,:,1)*tmask_i(:,:) )
-      DO jk = 2, ijpk
-         glob_min_3d = MIN( glob_min_3d, MINVAL( ptab(:,:,jk)*tmask_i(:,:) ) )
-      END DO
+      !!-----------------------------------------------------------------------
+      !
+      glob_min_3d = glob_min_3d_ref(ptab) 
       IF( lk_mpp )   CALL mpp_min( glob_min_3d )
       !
@@ -366,6 +481,6 @@
       !!-----------------------------------------------------------------------
       !             
-      glob_min_2d_a(1) = MINVAL( ptab1(:,:)*tmask_i(:,:) )
-      glob_min_2d_a(2) = MINVAL( ptab2(:,:)*tmask_i(:,:) )
+      glob_min_2d_a(1) = glob_min_2d_ref( ptab1 )
+      glob_min_2d_a(2) = glob_min_2d_ref( ptab2 )
       IF( lk_mpp )   CALL mpp_min( glob_min_2d_a, 2 )
       !
@@ -381,4 +496,58 @@
       REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2    ! input 3D array
       REAL(wp)            , DIMENSION(2)     ::   glob_min_3d_a   ! global masked min
+      !!-----------------------------------------------------------------------
+      !
+      glob_min_3d_a(1) = glob_min_3d_ref( ptab1 )
+      glob_min_3d_a(2) = glob_min_3d_ref( ptab2 )
+      IF( lk_mpp )   CALL mpp_min( glob_min_3d_a, 2 )
+      !
+   END FUNCTION glob_min_3d_a
+
+   ! --- MAX ---
+   FUNCTION glob_max_2d_ref( ptab ) 
+      !!-----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_max_2D  ***
+      !!
+      !! ** Purpose : perform a masked max on the inner global domain of a 2D array
+      !!-----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab          ! input 2D array
+      REAL(wp)                             ::   glob_max_2d_ref   ! global masked max
+      INTEGER                              ::   jj, ji        ! local index
+      !!-----------------------------------------------------------------------
+      !
+      glob_max_2d_ref = -1.e32
+!$OMP PARALLEL DO REDUCTION(MAX:glob_max_2d_ref)
+      DO jj = 1, jpj
+         DO ji =1, jpi
+            glob_max_2d_ref = MAX(glob_max_2d_ref, ptab(ji,jj)*tmask_i(ji,jj) )
+         ENDDO
+      ENDDO
+!$OMP END PARALLEL DO
+      !
+   END FUNCTION glob_max_2d_ref
+
+   FUNCTION glob_max_2d( ptab ) 
+      !!-----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_max_2D  ***
+      !!
+      !! ** Purpose : perform a masked max on the inner global domain of a 2D array
+      !!-----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab          ! input 2D array
+      REAL(wp)                             ::   glob_max_2d   ! global masked max
+      !!-----------------------------------------------------------------------
+      !
+      glob_max_2d = glob_max_2d_ref( ptab )
+      IF( lk_mpp )   CALL mpp_max( glob_max_2d )
+      !
+   END FUNCTION glob_max_2d
+
+   FUNCTION glob_max_3d_ref( ptab ) 
+      !!-----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_max_3D  ***
+      !!
+      !! ** Purpose : perform a masked max on the inner global domain of a 3D array
+      !!-----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab          ! input 3D array
+      REAL(wp)                               ::   glob_max_3d_ref   ! global masked max
       !!
       INTEGER :: jk
@@ -386,31 +555,14 @@
       !!-----------------------------------------------------------------------
       !
-      ijpk = SIZE(ptab1,3)
-      !
-      glob_min_3d_a(1) = MINVAL( ptab1(:,:,1)*tmask_i(:,:) )
-      glob_min_3d_a(2) = MINVAL( ptab2(:,:,1)*tmask_i(:,:) )
-      DO jk = 2, ijpk
-         glob_min_3d_a(1) = MIN( glob_min_3d_a(1), MINVAL( ptab1(:,:,jk)*tmask_i(:,:) ) )
-         glob_min_3d_a(2) = MIN( glob_min_3d_a(2), MINVAL( ptab2(:,:,jk)*tmask_i(:,:) ) )
+      ijpk = SIZE(ptab,3)
+      !
+      glob_max_3d_ref = -1e32
+!$OMP PARALLEL DO REDUCTION(MAX:glob_max_3d_ref)
+      DO jk = 1, ijpk
+         glob_max_3d_ref = MAX( glob_max_3d_ref, MAXVAL( ptab(:,:,jk)*tmask_i(:,:) ) )
       END DO
-      IF( lk_mpp )   CALL mpp_min( glob_min_3d_a, 2 )
-      !
-   END FUNCTION glob_min_3d_a
-
-   ! --- MAX ---
-   FUNCTION glob_max_2d( ptab ) 
-      !!-----------------------------------------------------------------------
-      !!                  ***  FUNCTION  glob_max_2D  ***
-      !!
-      !! ** Purpose : perform a masked max on the inner global domain of a 2D array
-      !!-----------------------------------------------------------------------
-      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab          ! input 2D array
-      REAL(wp)                             ::   glob_max_2d   ! global masked max
-      !!-----------------------------------------------------------------------
-      !
-      glob_max_2d = MAXVAL( ptab(:,:)*tmask_i(:,:) )
-      IF( lk_mpp )   CALL mpp_max( glob_max_2d )
-      !
-   END FUNCTION glob_max_2d
+!$OMP END PARALLEL DO 
+      !
+   END FUNCTION glob_max_3d_ref
  
    FUNCTION glob_max_3d( ptab ) 
@@ -422,15 +574,7 @@
       REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab          ! input 3D array
       REAL(wp)                               ::   glob_max_3d   ! global masked max
-      !!
-      INTEGER :: jk
-      INTEGER :: ijpk ! local variable: size of the 3d dimension of ptab
-      !!-----------------------------------------------------------------------
-      !
-      ijpk = SIZE(ptab,3)
-      !
-      glob_max_3d = MAXVAL( ptab(:,:,1)*tmask_i(:,:) )
-      DO jk = 2, ijpk
-         glob_max_3d = MAX( glob_max_3d, MAXVAL( ptab(:,:,jk)*tmask_i(:,:) ) )
-      END DO
+      !!-----------------------------------------------------------------------
+      !
+      glob_max_3d = glob_max_3d_ref( ptab )
       IF( lk_mpp )   CALL mpp_max( glob_max_3d )
       !
@@ -448,6 +592,6 @@
       !!-----------------------------------------------------------------------
       !             
-      glob_max_2d_a(1) = MAXVAL( ptab1(:,:)*tmask_i(:,:) )
-      glob_max_2d_a(2) = MAXVAL( ptab2(:,:)*tmask_i(:,:) )
+      glob_max_2d_a(1) = glob_max_2d_ref( ptab1 )
+      glob_max_2d_a(2) = glob_max_2d_ref( ptab2 )
       IF( lk_mpp )   CALL mpp_max( glob_max_2d_a, 2 )
       !
@@ -463,17 +607,8 @@
       REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2    ! input 3D array
       REAL(wp)            , DIMENSION(2)     ::   glob_max_3d_a   ! global masked max
-      !!
-      INTEGER :: jk
-      INTEGER :: ijpk ! local variable: size of the 3d dimension of ptab
-      !!-----------------------------------------------------------------------
-      !
-      ijpk = SIZE(ptab1,3)
-      !
-      glob_max_3d_a(1) = MAXVAL( ptab1(:,:,1)*tmask_i(:,:) )
-      glob_max_3d_a(2) = MAXVAL( ptab2(:,:,1)*tmask_i(:,:) )
-      DO jk = 2, ijpk
-         glob_max_3d_a(1) = MAX( glob_max_3d_a(1), MAXVAL( ptab1(:,:,jk)*tmask_i(:,:) ) )
-         glob_max_3d_a(2) = MAX( glob_max_3d_a(2), MAXVAL( ptab2(:,:,jk)*tmask_i(:,:) ) )
-      END DO
+      !!-----------------------------------------------------------------------
+      !
+      glob_max_3d_a(1) = glob_max_3d_ref( ptab1 )
+      glob_max_3d_a(2) = glob_max_3d_ref( ptab2 )
       IF( lk_mpp )   CALL mpp_max( glob_max_3d_a, 2 )
       !
Index: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/stpctl.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/stpctl.F90	(revision 9175)
+++ branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/stpctl.F90	(revision 9176)
@@ -156,4 +156,5 @@
          !
          IF(lwp) WRITE(numsol,9200) kt, niter, res, SQRT(epsr)/eps       ! Solver
+         IF(lwp) call flush(numsol)
          !
          IF( kindic < 0 .AND. zsmin > 0.e0 .AND. zumax <= 20.e0 ) THEN   ! create a abort file if problem found 
