Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 7698)
@@ -237,4 +237,5 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zh, zt, zs, ztm, zn3, zn2, zn1, zn0, zn)
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -277,4 +278,5 @@
       CASE( np_seos )                !==  simplified EOS  ==!
          !
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zh, zt, zs, ztm, zn)
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -345,4 +347,5 @@
             END DO
             !
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, jsmp, jdof, zh, zt, zstemp, zs, ztm, zn3, zn2, zn1)
             DO jk = 1, jpkm1
                DO jj = 1, jpj
@@ -399,4 +402,5 @@
          ! Non-stochastic equation of state
          ELSE
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zh, zt, zs, ztm, zn3, zn2, zn1, zn0, zn)
             DO jk = 1, jpkm1
                DO jj = 1, jpj
@@ -441,4 +445,5 @@
       CASE( np_seos )                !==  simplified EOS  ==!
          !
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zh, zt, zs, ztm, zn)
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -493,5 +498,10 @@
       IF( nn_timing == 1 )   CALL timing_start('eos2d')
       !
-      prd(:,:) = 0._wp
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            prd(ji,jj) = 0._wp
+         END DO
+      END DO
       !
       SELECT CASE( neos )
@@ -499,4 +509,5 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
+!$OMP PARALLEL DO schedule(static) private(jj, ji, zh, zt, zs, zn3, zn2, zn1, zn0, zn)
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
@@ -538,4 +549,5 @@
       CASE( np_seos )                !==  simplified EOS  ==!
          !
+!$OMP PARALLEL DO schedule(static) private(jj, ji, zh, zt, zs, zn)
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
@@ -589,4 +601,5 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zh, zt, zs, ztm, zn3, zn2, zn1, zn0, zn)
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -646,4 +659,5 @@
       CASE( np_seos )                  !==  simplified EOS  ==!
          !
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zh, zt, zs, ztm, zn)
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -698,5 +712,12 @@
       IF( nn_timing == 1 ) CALL timing_start('rab_2d')
       !
-      pab(:,:,:) = 0._wp
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+         DO jk = 1, jpts
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  pab(ji,jj,jk) = 0._wp
+               END DO
+            END DO
+         END DO
       !
       SELECT CASE ( neos )
@@ -704,4 +725,5 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
+!$OMP PARALLEL DO schedule(static) private(jj, ji, zh, zt, zs, zn3, zn2, zn1, zn0, zn)
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
@@ -762,4 +784,5 @@
       CASE( np_seos )                  !==  simplified EOS  ==!
          !
+!$OMP PARALLEL DO schedule(static) private(jj, ji, zh, zt, zs, zn)
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
@@ -917,4 +940,5 @@
       IF( nn_timing == 1 ) CALL timing_start('bn2')
       !
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, 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
@@ -952,6 +976,6 @@
       !!                Rational approximation to TEOS10 algorithm (rms error on WOA13 values: 4.0e-5 degC)
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ctmp   ! Cons. Temp   [Celsius]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity     [psu]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ctmp   ! Cons. Temp [Celsius]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity   [psu]
       ! Leave result array automatic rather than making explicitly allocated
       REAL(wp), DIMENSION(jpi,jpj) ::   ptmp   ! potential temperature [Celsius]
@@ -969,4 +993,5 @@
       z1_T0   = 1._wp/40._wp
       !
+!$OMP PARALLEL DO schedule(static) private(jj, ji, zt, zs, ztm, zn,zd)
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -1024,4 +1049,6 @@
          !
          z1_S0 = 1._wp / 35.16504_wp
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj, ji, zs)
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -1031,14 +1058,39 @@
             END DO
          END DO
-         ptf(:,:) = ptf(:,:) * psal(:,:)
-         !
-         IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:)
+!$OMP DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ptf(ji,jj) = ptf(ji,jj) * psal(ji,jj)
+            END DO
+         END DO
+!$OMP END PARALLEL
+         !
+         IF( PRESENT( pdep ) ) THEN
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+           DO jj = 1, jpj
+              DO ji = 1, jpi
+                 ptf(ji,jj) = ptf(ji,jj) - 7.53e-4 * pdep(ji,jj)
+              END DO
+           END DO
+         END IF
          !
       CASE ( np_eos80 )                !==  PT,SP (UNESCO formulation)  ==!
          !
-         ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   &
-            &                     - 2.154996e-4_wp *       psal(:,:)   ) * psal(:,:)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         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(:,:)
+         IF( PRESENT( pdep ) ) THEN
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+           DO jj = 1, jpj
+              DO ji = 1, jpi
+                 ptf(ji,jj) = ptf(ji,jj) - 7.53e-4 * pdep(ji,jj)
+              END DO
+           END DO
+         END IF
          !
       CASE DEFAULT
@@ -1134,4 +1186,5 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zh, zt, zs, ztm, zn2, zn1, zn0, zn)
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -1197,4 +1250,5 @@
       CASE( np_seos )                !==  Vallis (2006) simplified EOS  ==!
          !
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zh, zt, zs, ztm, zn)
          DO jk = 1, jpkm1
             DO jj = 1, jpj
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 7698)
@@ -88,5 +88,5 @@
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
       !
-      INTEGER ::   jk   ! dummy loop index
+      INTEGER :: ji, jj, jk   ! dummy loop index
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zun, zvn, zwn
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdt, ztrds   ! 3D workspace
@@ -98,7 +98,14 @@
       !
       !                                          ! set time step
-      zun(:,:,:) = 0.0
-      zvn(:,:,:) = 0.0
-      zwn(:,:,:) = 0.0
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zun(ji,jj,jk) = 0.0
+               zvn(ji,jj,jk) = 0.0
+               zwn(ji,jj,jk) = 0.0
+            END DO
+         END DO
+      END DO
       !    
       IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000
@@ -110,25 +117,47 @@
       !                                         !==  effective transport  ==!
       IF( ln_wave .AND. ln_sdw )  THEN
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
          DO jk = 1, jpkm1                                                       ! eulerian transport + Stokes Drift
-            zun(:,:,jk) = e2u  (:,:) * e3u_n(:,:,jk) * ( un(:,:,jk) + usd(:,:,jk) )
-            zvn(:,:,jk) = e1v  (:,:) * e3v_n(:,:,jk) * ( vn(:,:,jk) + vsd(:,:,jk) )
-            zwn(:,:,jk) = e1e2t(:,:)                 * ( wn(:,:,jk) + wsd(:,:,jk) )
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zun(ji,jj,jk) = e2u(ji,jj) * e3u_n(ji,jj,jk) * ( un(ji,jj,jk) + usd(ji,jj,jk) )
+                  zvn(ji,jj,jk) = e1v(ji,jj) * e3v_n(ji,jj,jk) * ( vn(ji,jj,jk) + vsd(ji,jj,jk) )
+                  zwn(ji,jj,jk) = e1e2t(ji,jj) * ( wn(ji,jj,jk) + wsd(ji,jj,jk) )
+               END DO
+            END DO
          END DO
       ELSE
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
          DO jk = 1, jpkm1
-            zun(:,:,jk) = e2u  (:,:) * e3u_n(:,:,jk) * un(:,:,jk)               ! eulerian transport only
-            zvn(:,:,jk) = e1v  (:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
-            zwn(:,:,jk) = e1e2t(:,:)                 * wn(:,:,jk)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zun(ji,jj,jk) = e2u  (ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)    ! eulerian transport only
+                  zvn(ji,jj,jk) = e1v  (ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
+                  zwn(ji,jj,jk) = e1e2t(ji,jj)                   * wn(ji,jj,jk)
+               END DO
+            END DO
          END DO
       ENDIF
       !
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN                                ! add z-tilde and/or vvl corrections
-         zun(:,:,:) = zun(:,:,:) + un_td(:,:,:)
-         zvn(:,:,:) = zvn(:,:,:) + vn_td(:,:,:)
-      ENDIF
-      !
-      zun(:,:,jpk) = 0._wp                                                      ! no transport trough the bottom
-      zvn(:,:,jpk) = 0._wp
-      zwn(:,:,jpk) = 0._wp
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zun(ji,jj,jk) = zun(ji,jj,jk) + un_td(ji,jj,jk)
+                  zvn(ji,jj,jk) = zvn(ji,jj,jk) + vn_td(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ENDIF
+      !
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zun(ji,jj,jpk) = 0._wp                                              ! no transport trough the bottom
+            zvn(ji,jj,jpk) = 0._wp
+            zwn(ji,jj,jpk) = 0._wp
+         END DO
+      END DO
       !
       IF( ln_ldfeiv .AND. .NOT. ln_traldf_triad )   &
@@ -147,6 +176,13 @@
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
          CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdt(ji,jj,jk) = tsa(ji,jj,jk,jp_tem)
+                  ztrds(ji,jj,jk) = tsa(ji,jj,jk,jp_sal)
+               END DO
+            END DO
+         END DO
       ENDIF
       !
@@ -169,7 +205,12 @@
       !
       IF( l_trdtra )   THEN                      ! save the advective trends for further diagnostics
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
          DO jk = 1, jpkm1
-            ztrdt(:,:,jk) = tsa(:,:,jk,jp_tem) - ztrdt(:,:,jk)
-            ztrds(:,:,jk) = tsa(:,:,jk,jp_sal) - ztrds(:,:,jk)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdt(ji,jj,jk) = tsa(ji,jj,jk,jp_tem) - ztrdt(ji,jj,jk)
+                  ztrds(ji,jj,jk) = tsa(ji,jj,jk,jp_sal) - ztrds(ji,jj,jk)
+               END DO
+            END DO
          END DO
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_totad, ztrdt )
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90	(revision 7698)
@@ -113,16 +113,48 @@
       IF( l_trd .OR. l_hst )  THEN
          CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
-         ztrdx(:,:,:) = 0._wp   ;    ztrdy(:,:,:) = 0._wp   ;   ztrdz(:,:,:) = 0._wp
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdx(ji,jj,jk) = 0._wp
+                  ztrdy(ji,jj,jk) = 0._wp
+                  ztrdz(ji,jj,jk) = 0._wp
+               END DO
+            END DO
+         END DO
       ENDIF
       !
       IF( l_ptr ) THEN  
          CALL wrk_alloc( jpi, jpj, jpk, zptry )
-         zptry(:,:,:) = 0._wp
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zptry(ji,jj,jk) = 0._wp
+               END DO
+            END DO
+         END DO
       ENDIF
       !                          ! surface & bottom value : flux set to zero one for all
-      zwz(:,:, 1 ) = 0._wp            
-      zwx(:,:,jpk) = 0._wp   ;   zwy(:,:,jpk) = 0._wp    ;    zwz(:,:,jpk) = 0._wp
-      !
-      zwi(:,:,:) = 0._wp        
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj, ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zwz(ji,jj, 1 ) = 0._wp
+            zwx(ji,jj,jpk) = 0._wp
+            zwy(ji,jj,jpk) = 0._wp
+            zwz(ji,jj,jpk) = 0._wp
+         END DO
+      END DO
+!$OMP END DO NOWAIT
+!$OMP DO schedule(static) private(jk, jj, ji)
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zwi(ji,jj,jk) = 0._wp
+            END DO
+         END DO
+      END DO
+!$OMP END PARALLEL
       !
       DO jn = 1, kjpt            !==  loop over the tracers  ==!
@@ -130,4 +162,6 @@
          !        !==  upstream advection with initial mass fluxes & intermediate update  ==!
          !                    !* upstream tracer flux in the i and j direction 
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jk, jj, ji, zfp_vj, zfm_vj, zfp_ui,zfm_ui)
          DO jk = 1, jpkm1
             DO jj = 1, jpjm1
@@ -143,5 +177,7 @@
             END DO
          END DO
+!$OMP END DO NOWAIT
          !                    !* upstream tracer flux in the k direction *!
+!$OMP DO schedule(static) private(jk, jj, ji, zfp_wk, zfm_wk)
          DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask)
             DO jj = 1, jpj
@@ -153,6 +189,8 @@
             END DO
          END DO
+!$OMP END PARALLEL
          IF( ln_linssh ) THEN    ! top ocean value (only in linear free surface as zwz has been w-masked)
             IF( ln_isfcav ) THEN             ! top of the ice-shelf cavities and at the ocean surface
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
                DO jj = 1, jpj
                   DO ji = 1, jpi
@@ -161,8 +199,14 @@
                END DO   
             ELSE                             ! no cavities: only at the ocean surface
-               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zwz(ji,jj,1) = pwn(ji,jj,1) * ptb(ji,jj,1,jn)
+                  END DO
+               END DO
             ENDIF
          ENDIF
          !               
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, ztra)
          DO jk = 1, jpkm1     !* trend and after field with monotonic scheme
             DO jj = 2, jpjm1
@@ -181,8 +225,26 @@
          !                
          IF( l_trd .OR. l_hst )  THEN             ! trend diagnostics (contribution of upstream fluxes)
-            ztrdx(:,:,:) = zwx(:,:,:)   ;   ztrdy(:,:,:) = zwy(:,:,:)   ;   ztrdz(:,:,:) = zwz(:,:,:)
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdx(ji,jj,jk) = zwx(ji,jj,jk)
+                     ztrdy(ji,jj,jk) = zwy(ji,jj,jk)
+                     ztrdz(ji,jj,jk) = zwz(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
          END IF
          !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( l_ptr )  zptry(:,:,:) = zwy(:,:,:) 
+         IF( l_ptr ) THEN
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zptry(ji,jj,jk) = zwy(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         END IF
          !
          !        !==  anti-diffusive flux : high order minus low order  ==!
@@ -191,4 +253,5 @@
          !
          CASE(  2  )                   !- 2nd order centered
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
             DO jk = 1, jpkm1
                DO jj = 1, jpjm1
@@ -201,6 +264,13 @@
             !
          CASE(  4  )                   !- 4th order centered
-            zltu(:,:,jpk) = 0._wp            ! Bottom value : flux set to zero
-            zltv(:,:,jpk) = 0._wp
+!$OMP PARALLEL 
+!$OMP DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zltu(ji,jj,jpk) = 0._wp            ! Bottom value : flux set to zero
+                  zltv(ji,jj,jpk) = 0._wp
+               END DO
+            END DO
+!$OMP DO schedule(static) private(jk, jj, ji)
             DO jk = 1, jpkm1                 ! Laplacian
                DO jj = 1, jpjm1                    ! 1st derivative (gradient)
@@ -217,6 +287,8 @@
                END DO
             END DO
+!$OMP END PARALLEL
             CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
             !
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zC2t_u, zC2t_v)
             DO jk = 1, jpkm1                 ! Horizontal advective fluxes
                DO jj = 1, jpjm1
@@ -232,6 +304,13 @@
             !
          CASE(  41 )                   !- 4th order centered       ==>>   !!gm coding attempt   need to be tested
-            ztu(:,:,jpk) = 0._wp             ! Bottom value : flux set to zero
-            ztv(:,:,jpk) = 0._wp
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztu(ji,jj,jpk) = 0._wp             ! Bottom value : flux set to zero
+                  ztv(ji,jj,jpk) = 0._wp
+               END DO
+            END DO
+!$OMP DO schedule(static) private(jk, jj, ji)
             DO jk = 1, jpkm1                 ! 1st derivative (gradient)
                DO jj = 1, jpjm1
@@ -242,6 +321,8 @@
                END DO
             END DO
+!$OMP END PARALLEL
             CALL lbc_lnk( ztu, 'U', -1. )   ;    CALL lbc_lnk( ztv, 'V', -1. )   ! Lateral boundary cond. (unchanged sgn)
             !
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zC2t_u, zC2t_v, zC4t_u, zC4t_v)
             DO jk = 1, jpkm1                 ! Horizontal advective fluxes
                DO jj = 2, jpjm1
@@ -264,4 +345,5 @@
          !
          CASE(  2  )                   !- 2nd order centered
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
             DO jk = 2, jpkm1    
                DO jj = 2, jpjm1
@@ -275,4 +357,5 @@
          CASE(  4  )                   !- 4th order COMPACT
             CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )   ! zwt = COMPACT interpolation of T at w-point
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
             DO jk = 2, jpkm1
                DO jj = 2, jpjm1
@@ -285,5 +368,10 @@
          END SELECT
          IF( ln_linssh ) THEN    ! top ocean value: high order = upstream  ==>>  zwz=0
-            zwz(:,:,1) = 0._wp   ! only ocean surface as interior zwz values have been w-masked
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zwz(ji,jj,1) = 0._wp   ! only ocean surface as interior zwz values have been w-masked
+               END DO
+            END DO
          ENDIF
          !
@@ -297,4 +385,5 @@
          !        !==  final trend with corrected fluxes  ==!
          !
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -309,7 +398,14 @@
          !
          IF( l_trd .OR. l_hst ) THEN     ! trend diagnostics (contribution of upstream fluxes)
-            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed
-            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
-            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdx(ji,jj,jk) = ztrdx(ji,jj,jk) + zwx(ji,jj,jk)  ! <<< Add to previously computed
+                     ztrdy(ji,jj,jk) = ztrdy(ji,jj,jk) + zwy(ji,jj,jk)  ! <<< Add to previously computed
+                     ztrdz(ji,jj,jk) = ztrdz(ji,jj,jk) + zwz(ji,jj,jk)  ! <<< Add to previously computed
+                  END DO
+               END DO
+            END DO
          ENDIF
             !
@@ -325,5 +421,12 @@
          !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
          IF( l_ptr ) THEN  
-            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zptry(ji,jj,jk) = zptry(ji,jj,jk) + zwy(ji,jj,jk)  ! <<< Add to previously computed
+                  END DO
+               END DO
+            END DO
             CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
          ENDIF
@@ -662,5 +765,4 @@
       zbig  = 1.e+40_wp
       zrtrn = 1.e-15_wp
-      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
 
       ! Search local extrema
@@ -672,4 +774,15 @@
          &        paft * tmask + zbig * ( 1._wp - tmask )  )
 
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jk, jj, ji)
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zbetup(ji,jj,jk) = 0._wp
+               zbetdo(ji,jj,jk) = 0._wp
+            END DO
+         END DO
+      END DO
+!$OMP DO schedule(static) private(jk, jj, ji, ikm1, zup, zdo, zpos, zneg, zbt)
       DO jk = 1, jpkm1
          ikm1 = MAX(jk-1,1)
@@ -706,8 +819,10 @@
          END DO
       END DO
+!$OMP END PARALLEL
       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 schedule(static) private(jk, jj, ji, za, zb, zc, zav, zbv, zcv, zau, zbu, zcu)
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mle.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mle.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mle.F90	(revision 7698)
@@ -327,4 +327,5 @@
             IF( ierr /= 0 )   CALL ctl_stop( 'tra_adv_mle_init: failed to allocate arrays' )
             z1_t2 = 1._wp / ( rn_time * rn_time )
+!$OMP PARALLEL DO schedule(static) private(jj, ji, zfu, zfv)
             DO jj = 2, jpj                           ! "coriolis+ time^-1" at u- & v-points
                DO ji = fs_2, jpi   ! vector opt.
@@ -347,5 +348,10 @@
          !
          z1_t2 = 1._wp / ( rn_time * rn_time )
-         r1_ft(:,:) = 1._wp / SQRT(  ff_t(:,:) * ff_t(:,:) + z1_t2  )
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               r1_ft(ji,jj) = 1._wp / SQRT(  ff_t(ji,jj) * ff_t(ji,jj) + z1_t2 )
+            END DO
+         END DO
          !
       ENDIF
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90	(revision 7698)
@@ -108,15 +108,35 @@
          !
          ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr )
-         xind(:,:,:) = 1._wp              ! set equal to 1 where up-stream is not needed
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  xind(ji,jj,jk) = 1._wp              ! set equal to 1 where up-stream is not needed
+               END DO
+            END DO
+         END DO
          !
          IF( ld_msc_ups ) THEN            ! define the upstream indicator (if asked)
             ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
-            upsmsk(:,:) = 0._wp                             ! not upstream by default
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  upsmsk(ji,jj) = 0._wp                             ! not upstream by default
+               END DO
+            END DO
             !
+!$OMP DO schedule(static) private(jk,jj,ji)
             DO jk = 1, jpkm1
-               xind(:,:,jk) = 1._wp                              &                 ! =>1 where up-stream is not needed
-                  &         - MAX ( rnfmsk(:,:) * rnfmsk_z(jk),  &                 ! =>0 near runoff mouths (& closed sea outflows)
-                  &                 upsmsk(:,:)                ) * tmask(:,:,jk)   ! =>0 in some user defined area
-            END DO
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     xind(ji,jj,jk) = 1._wp                              &                   ! =>1 where up-stream is not needed
+                        &         - MAX ( rnfmsk(ji,jj) * rnfmsk_z(jk),  &                   ! =>0 near runoff mouths (& closed sea outflows)
+                        &                 upsmsk(ji,jj)                ) * tmask(ji,jj,jk)   ! =>0 in some user defined area
+                  END DO
+               END DO
+            END DO
+!$OMP END DO NOWAIT
+!$OMP END PARALLEL
          ENDIF 
          !
@@ -136,6 +156,13 @@
          !
          !                                !-- first guess of the slopes
-         zwx(:,:,jpk) = 0._wp                   ! bottom values
-         zwy(:,:,jpk) = 0._wp  
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zwx(ji,jj,jpk) = 0._wp           ! bottom values
+               zwy(ji,jj,jpk) = 0._wp
+            END DO
+         END DO
+!$OMP DO schedule(static) private(jk, jj, ji)
          DO jk = 1, jpkm1                       ! interior values
             DO jj = 1, jpjm1      
@@ -146,9 +173,18 @@
            END DO
          END DO
+!$OMP END DO NOWAIT
+!$OMP END PARALLEL
          CALL lbc_lnk( zwx, 'U', -1. )          ! lateral boundary conditions   (changed sign)
          CALL lbc_lnk( zwy, 'V', -1. )
          !                                !-- Slopes of tracer
-         zslpx(:,:,jpk) = 0._wp                 ! bottom values
-         zslpy(:,:,jpk) = 0._wp
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zslpx(ji,jj,jpk) = 0._wp                 ! bottom values
+               zslpy(ji,jj,jpk) = 0._wp
+            END DO
+         END DO
+!$OMP DO schedule(static) private(jk, jj, ji)
          DO jk = 1, jpkm1                       ! interior values
             DO jj = 2, jpj
@@ -162,4 +198,5 @@
          END DO
          !
+!$OMP DO schedule(static) private(jk, jj, ji)
          DO jk = 1, jpkm1                 !-- Slopes limitation
             DO jj = 2, jpj
@@ -175,4 +212,5 @@
          END DO
          !
+!$OMP DO schedule(static) private(jk, jj, ji, z0u, zalpha, zu, zv, zzwx, zzwy, z0v)
          DO jk = 1, jpkm1                 !-- MUSCL horizontal advective fluxes
             DO jj = 2, jpjm1
@@ -195,6 +233,9 @@
             END DO
          END DO
+!$OMP END DO NOWAIT
+!$OMP END PARALLEL
          CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )   ! lateral boundary conditions   (changed sign)
          !
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
          DO jk = 1, jpkm1                 !-- Tracer advective trend
             DO jj = 2, jpjm1      
@@ -219,11 +260,29 @@
          !
          !                                !-- first guess of the slopes
-         zwx(:,:, 1 ) = 0._wp                   ! surface & bottom boundary conditions
-         zwx(:,:,jpk) = 0._wp
+!$OMP PARALLEL 
+!$OMP DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zwx(ji,jj, 1 ) = 0._wp           ! surface & bottom boundary conditions
+               zwx(ji,jj,jpk) = 0._wp
+           END DO
+         END DO
+!$OMP DO schedule(static) private(jk, jj, ji)
          DO jk = 2, jpkm1                       ! interior values
-            zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) )
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zwx(ji,jj,jk) = tmask(ji,jj,jk) * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
+              END DO
+            END DO
          END DO
          !                                !-- Slopes of tracer
-         zslpx(:,:,1) = 0._wp                   ! surface values
+!$OMP END DO NOWAIT
+!$OMP DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zslpx(ji,jj,1) = 0._wp                   ! surface values
+           END DO
+         END DO
+!$OMP DO schedule(static) private(jk, jj, ji)
          DO jk = 2, jpkm1                       ! interior value
             DO jj = 1, jpj
@@ -234,4 +293,5 @@
             END DO
          END DO
+!$OMP DO schedule(static) private(jk, jj, ji)
          DO jk = 2, jpkm1                 !-- Slopes limitation
             DO jj = 1, jpj                      ! interior values
@@ -243,4 +303,5 @@
             END DO
          END DO
+!$OMP DO schedule(static) private(jk, jj, ji, z0w, zalpha, zw, zzwx, zzwy)
          DO jk = 1, jpk-2                 !-- vertical advective flux
             DO jj = 2, jpjm1      
@@ -255,6 +316,9 @@
             END DO
          END DO
+!$OMP END DO NOWAIT
+!$OMP END PARALLEL
          IF( ln_linssh ) THEN                   ! top values, linear free surface only
             IF( ln_isfcav ) THEN                      ! ice-shelf cavities (top of the ocean)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
                DO jj = 1, jpj
                   DO ji = 1, jpi
@@ -263,8 +327,14 @@
                END DO   
             ELSE                                      ! no cavities: only at the ocean surface
-               zwx(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zwx(ji,jj,1) = pwn(ji,jj,1) * ptb(ji,jj,1,jn)
+                  END DO
+               END DO
             ENDIF
          ENDIF
          !
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
          DO jk = 1, jpkm1                 !-- vertical advective trend
             DO jj = 2, jpjm1      
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 7698)
@@ -76,5 +76,5 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
-      INTEGER  ::   ji, jj    ! dummy loop indices
+      INTEGER  ::   ji, jj, jk    ! dummy loop indices
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdt
       !!----------------------------------------------------------------------
@@ -84,7 +84,15 @@
       IF( l_trdtra )   THEN         ! Save the input temperature trend
          CALL wrk_alloc( jpi,jpj,jpk,   ztrdt )
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdt(ji,jj,jk) = tsa(ji,jj,jk,jp_tem)
+               END DO
+            END DO
+         END DO
       ENDIF
       !                             !  Add the geothermal trend on temperature
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
@@ -96,5 +104,12 @@
       !
       IF( l_trdtra ) THEN        ! Send the trend for diagnostics
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdt(ji,jj,jk) = tsa(ji,jj,jk,jp_tem) - ztrdt(ji,jj,jk)
+               END DO
+            END DO
+         END DO
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbc, ztrdt )
          CALL wrk_dealloc( jpi,jpj,jpk,   ztrdt )
@@ -162,5 +177,10 @@
          CASE ( 1 )                          !* constant flux
             IF(lwp) WRITE(numout,*) '      *** constant heat flux  =   ', rn_geoflx_cst
-            qgh_trd0(:,:) = r1_rau0_rcp * rn_geoflx_cst
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  qgh_trd0(ji,jj) = r1_rau0_rcp * rn_geoflx_cst
+               END DO
+            END DO
             !
          CASE ( 2 )                          !* variable geothermal heat flux : read the geothermal fluxes in mW/m2
@@ -179,5 +199,10 @@
 
             CALL fld_read( nit000, 1, sf_qgh )                         ! Read qgh data
-            qgh_trd0(:,:) = r1_rau0_rcp * sf_qgh(1)%fnow(:,:,1) * 1.e-3 ! conversion in W/m2
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  qgh_trd0(ji,jj) = r1_rau0_rcp * sf_qgh(1)%fnow(ji,jj,1) * 1.e-3 ! conversion in W/m2
+               END DO
+            END DO
             !
          CASE DEFAULT
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 7698)
@@ -105,4 +105,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step
+      INTEGER  ::   ji, jj, jk        ! dummy loop indices
       !
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
@@ -113,6 +114,13 @@
       IF( l_trdtra )   THEN                         !* Save the input trends
          CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdt(ji,jj,jk) = tsa(ji,jj,jk,jp_tem)
+                  ztrds(ji,jj,jk) = tsa(ji,jj,jk,jp_sal)
+               END DO
+            END DO
+         END DO
       ENDIF
 
@@ -146,6 +154,13 @@
 
       IF( l_trdtra )   THEN                      ! send the trends for further diagnostics
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdt(ji,jj,jk) = tsa(ji,jj,jk,jp_tem) - ztrdt(ji,jj,jk)
+                  ztrds(ji,jj,jk) = tsa(ji,jj,jk,jp_sal) - ztrds(ji,jj,jk)
+               END DO
+            END DO
+         END DO
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrdt )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrds )
@@ -195,4 +210,5 @@
       DO jn = 1, kjpt                                     ! tracer loop
          !                                                ! ===========
+!$OMP PARALLEL DO schedule(static) private(jj,ji,ik)
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -202,4 +218,5 @@
          END DO
          !               
+!$OMP PARALLEL DO schedule(static) private(jj,ji,ik)
          DO jj = 2, jpjm1                                    ! Compute the trend
             DO ji = 2, jpim1
@@ -357,4 +374,5 @@
       ENDIF
       !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
+!$OMP PARALLEL DO schedule(static) private(jj,ji,ik)
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -374,4 +392,5 @@
       IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
          !                                !-------------------!
+!$OMP PARALLEL DO schedule(static) private(jj,ji,za,zb,zgdrho,zsign)
          DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
             DO ji = 1, fs_jpim1   ! vector opt.
@@ -406,4 +425,5 @@
          !
          CASE( 1 )                                   != use of upper velocity
+!$OMP PARALLEL DO schedule(static) private(jj,ji,za,zb,zgdrho,zsign,zsigna)
             DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
                DO ji = 1, fs_jpim1   ! vector opt.
@@ -437,4 +457,5 @@
          CASE( 2 )                                 != bbl velocity = F( delta rho )
             zgbbl = grav * rn_gambbl
+!$OMP PARALLEL DO schedule(static) private(jj,ji,iid,iis,ikud,ikus,za,zb,zgdrho,ijd,ijs,ikvd,ikvs)
             DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
                DO ji = 1, fs_jpim1   ! vector opt.
@@ -533,4 +554,5 @@
 
       !                             !* vertical index of  "deep" bottom u- and v-points
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
       DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
          DO ji = 1, jpim1
@@ -547,4 +569,5 @@
       !                                 !* sign of grad(H) at u- and v-points
       mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
       DO jj = 1, jpjm1
          DO ji = 1, jpim1
@@ -554,4 +577,5 @@
       END DO
       !
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
       DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
          DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
@@ -563,6 +587,11 @@
       !
       !                             !* masked diffusive flux coefficients
-      ahu_bbl_0(:,:) = rn_ahtbbl * e2_e1u(:,:) * e3u_bbl_0(:,:) * umask(:,:,1)
-      ahv_bbl_0(:,:) = rn_ahtbbl * e1_e2v(:,:) * e3v_bbl_0(:,:) * vmask(:,:,1)
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ahu_bbl_0(ji,jj) = rn_ahtbbl * e2_e1u(ji,jj) * e3u_bbl_0(ji,jj) * umask(ji,jj,1)
+            ahv_bbl_0(ji,jj) = rn_ahtbbl * e1_e2v(ji,jj) * e3v_bbl_0(ji,jj) * vmask(ji,jj,1)
+         END DO
+      END DO
 
       !
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 7698)
@@ -102,5 +102,14 @@
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
          CALL wrk_alloc( jpi,jpj,jpk,jpts,   ztrdts ) 
-         ztrdts(:,:,:,:) = tsa(:,:,:,:) 
+         DO jn = 1, jpts
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdts(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) 
+                  END DO
+               END DO
+            END DO
+         END DO
       ENDIF
       !                           !==  input T-S data at kt  ==!
@@ -111,4 +120,5 @@
       CASE( 0 )                        !*  newtonian damping throughout the water column  *!
          DO jn = 1, jpts
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
             DO jk = 1, jpkm1
                DO jj = 2, jpjm1
@@ -121,4 +131,5 @@
          !
       CASE ( 1 )                       !*  no damping in the turbocline (avt > 5 cm2/s)  *!
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -135,4 +146,5 @@
          !
       CASE ( 2 )                       !*  no damping in the mixed layer   *!
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -151,5 +163,14 @@
       !
       IF( l_trdtra )   THEN       ! trend diagnostic
-         ztrdts(:,:,:,:) = tsa(:,:,:,:) - ztrdts(:,:,:,:)
+         DO jn = 1, jpts
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdts(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) - ztrdts(ji,jj,jk,jn)
+                  END DO
+               END DO
+            END DO
+         END DO
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_dmp, ztrdts(:,:,:,jp_tem) )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_dmp, ztrdts(:,:,:,jp_sal) )
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 7698)
@@ -57,4 +57,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      INTEGER ::   jk, jj, ji         ! dummy loop indices
       !!
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
@@ -65,6 +66,13 @@
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
          CALL wrk_alloc( jpi,jpj,jpk,   ztrdt, ztrds ) 
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdt(ji,jj,jk) = tsa(ji,jj,jk,jp_tem)
+                  ztrds(ji,jj,jk) = tsa(ji,jj,jk,jp_sal)
+               END DO
+            END DO
+         END DO
       ENDIF
       !
@@ -81,6 +89,13 @@
       !
       IF( l_trdtra )   THEN                    !* save the horizontal diffusive trends for further diagnostics
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdt(ji,jj,jk) = tsa(ji,jj,jk,jp_tem) - ztrdt(ji,jj,jk)
+                  ztrds(ji,jj,jk) = tsa(ji,jj,jk,jp_sal) - ztrds(ji,jj,jk)
+               END DO
+            END DO
+         END DO
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrdt )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_ldf, ztrds )
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 7698)
@@ -125,6 +125,13 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
          !
-         akz     (:,:,:) = 0._wp      
-         ah_wslp2(:,:,:) = 0._wp
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  akz     (ji,jj,jk) = 0._wp
+                  ah_wslp2(ji,jj,jk) = 0._wp
+               END DO
+            END DO
+         END DO
       ENDIF
       !   
@@ -151,4 +158,5 @@
       IF( kpass == 1 ) THEN                  !==  first pass only  ==!
          !
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zmsku, zmskv, zahu_w, zahv_w)
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
@@ -172,4 +180,5 @@
          !
          IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
             DO jk = 2, jpkm1
                DO jj = 2, jpjm1
@@ -185,4 +194,5 @@
             !
             IF( ln_traldf_blp ) THEN                ! bilaplacian operator
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
                DO jk = 2, jpkm1
                   DO jj = 1, jpjm1
@@ -194,4 +204,5 @@
                END DO
             ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, ze3w_2, zcoef0)
                DO jk = 2, jpkm1
                   DO jj = 1, jpjm1
@@ -206,5 +217,12 @@
            !
          ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
-            akz(:,:,:) = ah_wslp2(:,:,:)      
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+           DO jk = 1, jpk
+              DO jj = 1, jpj
+                 DO ji = 1, jpi
+                    akz(ji,jj,jk) = ah_wslp2(ji,jj,jk)
+                 END DO
+              END DO
+           END DO
          ENDIF
       ENDIF
@@ -218,9 +236,16 @@
          !!----------------------------------------------------------------------
 !!gm : bug.... why (x,:,:)?   (1,jpj,:) and (jpi,1,:) should be sufficient....
-         zdit (1,:,:) = 0._wp     ;     zdit (jpi,:,:) = 0._wp
-         zdjt (1,:,:) = 0._wp     ;     zdjt (jpi,:,:) = 0._wp
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jk, jj)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               zdit (1,jj,jk) = 0._wp     ;     zdit (jpi,jj,jk) = 0._wp
+               zdjt (1,jj,jk) = 0._wp     ;     zdjt (jpi,jj,jk) = 0._wp
+            END DO
+         END DO
          !!end
 
          ! Horizontal tracer gradient 
+!$OMP DO schedule(static) private(jk, jj, ji)
          DO jk = 1, jpkm1
             DO jj = 1, jpjm1
@@ -231,5 +256,7 @@
             END DO
          END DO
+!$OMP END PARALLEL
          IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
             DO jj = 1, jpjm1              ! bottom correction (partial bottom cell)
                DO ji = 1, fs_jpim1   ! vector opt.
@@ -239,4 +266,5 @@
             END DO
             IF( ln_isfcav ) THEN      ! first wet level beneath a cavity
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
                DO jj = 1, jpjm1
                   DO ji = 1, fs_jpim1   ! vector opt.
@@ -252,12 +280,31 @@
          !!----------------------------------------------------------------------
          !
+!$OMP PARALLEL
          DO jk = 1, jpkm1                                 ! Horizontal slab
             !
             !                             !== Vertical tracer gradient
-            zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * wmask(:,:,jk+1)     ! level jk+1
-            !
-            IF( jk == 1 ) THEN   ;   zdkt(:,:) = zdk1t(:,:)                          ! surface: zdkt(jk=1)=zdkt(jk=2)
-            ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk)
+!$OMP DO schedule(static) private(jj, ji)
+            DO jj = 1 , jpj
+               DO ji = 1, jpi
+                  zdk1t(ji,jj) = ( ptb(ji,jj,jk,jn) - ptb(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1)     ! level jk+1
+               END DO
+            END DO
+            !
+            IF( jk == 1 ) THEN   
+!$OMP DO schedule(static) private(jj, ji)
+               DO jj = 1 , jpj
+                  DO ji = 1, jpi
+                     zdkt(ji,jj) = zdk1t(ji,jj)                          ! surface: zdkt(jk=1)=zdkt(jk=2)
+                  END DO
+               END DO
+            ELSE  
+!$OMP DO schedule(static) private(jj, ji)
+               DO jj = 1 , jpj
+                  DO ji = 1, jpi
+                     zdkt(ji,jj) = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * wmask(ji,jj,jk)
+                  END DO
+               END DO
             ENDIF
+!$OMP DO schedule(static) private(jj, ji, zmsku, zmskv, zabe1, zabe2, zcof1, zcof2)
             DO jj = 1 , jpjm1            !==  Horizontal fluxes
                DO ji = 1, fs_jpim1   ! vector opt.
@@ -283,4 +330,5 @@
             END DO
             !
+!$OMP DO schedule(static) private(jj, ji)
             DO jj = 2 , jpjm1          !== horizontal divergence and add to pta
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -296,11 +344,22 @@
          !!----------------------------------------------------------------------
          !
-         ztfw(1,:,:) = 0._wp     ;     ztfw(jpi,:,:) = 0._wp
+!$OMP DO schedule(static) private(jk, jj)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               ztfw(1,jj,jk) = 0._wp     ;     ztfw(jpi,jj,jk) = 0._wp
+            END DO
+         END DO
          !
          ! Vertical fluxes
          ! ---------------
          !                          ! Surface and bottom vertical fluxes set to zero
-         ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp
+!$OMP DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ztfw(ji,jj, 1 ) = 0._wp      ;      ztfw(ji,jj,jpk) = 0._wp
+            END DO
+         END DO
          
+!$OMP DO schedule(static) private(jk, jj, ji, zmsku, zmskv, zahu_w, zahv_w, zcoef3, zcoef4)
          DO jk = 2, jpkm1           ! interior (2=<jk=<jpk-1)
             DO jj = 2, jpjm1
@@ -327,6 +386,8 @@
             END DO
          END DO
+!$OMP END PARALLEL
          !                                !==  add the vertical 33 flux  ==!
          IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
             DO jk = 2, jpkm1       
                DO jj = 1, jpjm1
@@ -342,4 +403,5 @@
             SELECT CASE( kpass )
             CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
                DO jk = 2, jpkm1 
                   DO jj = 1, jpjm1
@@ -352,4 +414,5 @@
                END DO 
             CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp.
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
                DO jk = 2, jpkm1 
                   DO jj = 1, jpjm1
@@ -364,4 +427,5 @@
          ENDIF
          !         
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
          DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==!
             DO jj = 2, jpjm1
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 7698)
@@ -121,6 +121,13 @@
       IF( l_trdtra )   THEN                    
          CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
-         ztrdt(:,:,jk) = 0._wp
-         ztrds(:,:,jk) = 0._wp
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdt(ji,jj,jk) = 0._wp 
+                  ztrds(ji,jj,jk) = 0._wp
+               END DO
+            END DO
+         END DO
          IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend 
             CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdfp, ztrdt )
@@ -129,7 +136,12 @@
          ! total trend for the non-time-filtered variables. 
             zfact = 1.0 / rdt
-         DO jk = 1, jpkm1
-            ztrdt(:,:,jk) = ( tsa(:,:,jk,jp_tem) - tsn(:,:,jk,jp_tem) ) * zfact 
-            ztrds(:,:,jk) = ( tsa(:,:,jk,jp_sal) - tsn(:,:,jk,jp_sal) ) * zfact 
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdt(ji,jj,jk) = ( tsa(ji,jj,jk,jp_tem) - tsn(ji,jj,jk,jp_tem) ) * zfact
+                  ztrds(ji,jj,jk) = ( tsa(ji,jj,jk,jp_sal) - tsn(ji,jj,jk,jp_sal) ) * zfact
+               END DO
+            END DO
          END DO
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_tot, ztrdt )
@@ -137,12 +149,24 @@
          ! Store now fields before applying the Asselin filter 
          ! in order to calculate Asselin filter trend later.
-         ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
-         ztrds(:,:,:) = tsn(:,:,:,jp_sal)
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdt(ji,jj,jk) = tsn(ji,jj,jk,jp_tem)
+                  ztrds(ji,jj,jk) = tsn(ji,jj,jk,jp_sal)
+               END DO
+            END DO
+         END DO
       ENDIF
 
       IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
          DO jn = 1, jpts
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
             DO jk = 1, jpkm1
-               tsn(:,:,jk,jn) = tsa(:,:,jk,jn)    
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     tsn(ji,jj,jk,jn) = tsa(ji,jj,jk,jn)    
+                  END DO
+               END DO
             END DO
          END DO
@@ -163,8 +187,13 @@
       !
       IF( l_trdtra ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
-         DO jk = 1, jpkm1
-            zfact = 1._wp / r2dt             
-            ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
-            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
+!$OMP PARALLEL DO schedule(static) private(jk, zfact)
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zfact = 1._wp / r2dt             
+                  ztrdt(ji,jj,jk) = ( tsb(ji,jj,jk,jp_tem) - ztrdt(ji,jj,jk) ) * zfact
+                  ztrds(ji,jj,jk) = ( tsb(ji,jj,jk,jp_sal) - ztrds(ji,jj,jk) ) * zfact
+               END DO
+            END DO
          END DO
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
@@ -214,4 +243,5 @@
       DO jn = 1, kjpt
          !
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,ztn,ztd)
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -280,4 +310,5 @@
       !
       DO jn = 1, kjpt      
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,zfact1,zfact2,ze3t_b,ze3t_n,ze3t_a,ze3t_d,ze3t_f,ztc_b,ztc_n,ztc_a,ztc_d,ztc_f)
          DO jk = 1, jpkm1
             zfact1 = atfp * p2dt
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 7698)
@@ -128,5 +128,12 @@
       IF( l_trdtra ) THEN      ! trends diagnostic: save the input temperature trend
          CALL wrk_alloc( jpi,jpj,jpk,   ztrdt ) 
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdt(ji,jj,jk) = tsa(ji,jj,jk,jp_tem)
+                  END DO
+               END DO
+            END DO
       ENDIF
       !
@@ -142,9 +149,23 @@
          ELSE                                           ! No restart or restart not found: Euler forward time stepping
             z1_2 = 1._wp
-            qsr_hc_b(:,:,:) = 0._wp
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     qsr_hc_b(ji,jj,jk) = 0._wp
+                  END DO
+               END DO
+            END DO
          ENDIF
       ELSE                             !==  Swap of qsr heat content  ==!
          z1_2 = 0.5_wp
-         qsr_hc_b(:,:,:) = qsr_hc(:,:,:)
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     qsr_hc_b(ji,jj,jk) = qsr_hc(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
       ENDIF
       !
@@ -155,6 +176,11 @@
       CASE( np_BIO )                   !==  bio-model fluxes  ==!
          !
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
          DO jk = 1, nksr
-            qsr_hc(:,:,jk) = r1_rau0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  qsr_hc(ji,jj,jk) = r1_rau0_rcp * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) )
+               END DO
+             END DO
          END DO
          !
@@ -166,4 +192,5 @@
          IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll
             CALL fld_read( kt, 1, sf_chl )         ! Read Chl data and provides it at the current time step
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,zchl,zCtot,zze,zpsi,zlogc,zlogc2,zlogc3,zCb,zCmax,zpsimax,zdelpsi,zCze)
             DO jk = 1, nksr + 1
                DO jj = 2, jpjm1                       ! Separation in R-G-B depending of the surface Chl
@@ -190,10 +217,17 @@
             END DO
          ELSE                                !* constant chrlorophyll
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
            DO jk = 1, nksr + 1
-              zchl3d(:,:,jk) = 0.05 
+              DO jj = 1, jpj
+                 DO ji = 1, jpi
+                    zchl3d(ji,jj,jk) = 0.05
+                 ENDDO
+              ENDDO
             ENDDO
          ENDIF
          !
          zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj,ji)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
@@ -205,6 +239,8 @@
             END DO
          END DO
+!$OMP END DO NOWAIT
          !
          DO jk = 2, nksr+1                   !* interior equi-partition in R-G-B depending of vertical profile of Chl
+!$OMP DO schedule(static) private(jj,ji,zchl,irgb)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
@@ -217,4 +253,5 @@
             END DO
 
+!$OMP DO schedule(static) private(jj,ji,zc0,zc1,zc2,zc3)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
@@ -232,4 +269,5 @@
          END DO
          !
+!$OMP DO schedule(static) private(jk,jj,ji)
          DO jk = 1, nksr                     !* now qsr induced heat content
             DO jj = 2, jpjm1
@@ -239,4 +277,5 @@
             END DO
          END DO
+!$OMP END PARALLEL
          !
          CALL wrk_dealloc( jpi,jpj,        zekb, zekg, zekr        ) 
@@ -247,4 +286,5 @@
          zz0 =        rn_abs   * r1_rau0_rcp      ! surface equi-partition in 2-bands
          zz1 = ( 1. - rn_abs ) * r1_rau0_rcp
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,zc0,zc1)
          DO jk = 1, nksr                          ! solar heat absorbed at T-point in the top 400m 
             DO jj = 2, jpjm1
@@ -260,4 +300,5 @@
       !
       !                          !-----------------------------!
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
       DO jk = 1, nksr            !  update to the temp. trend  !
          DO jj = 2, jpjm1        !-----------------------------!
@@ -270,4 +311,5 @@
       !
       IF( ln_qsr_ice ) THEN      ! sea-ice: store the 1st ocean level attenuation coefficient
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
          DO jj = 2, jpjm1 
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -284,8 +326,20 @@
          CALL wrk_alloc( jpi,jpj,jpk,   zetot )
          !
-         zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj,ji)
+         DO jj = 1, jpj 
+            DO ji = 1, jpi   ! vector opt.
+               zetot(ji,jj,nksr+1:jpk) = 0._wp     ! below ~400m set to zero
+            END DO
+         END DO
          DO jk = nksr, 1, -1
-            zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) / r1_rau0_rcp
+!$OMP DO schedule(static) private(jj,ji)
+            DO jj = 1, jpj 
+               DO ji = 1, jpi   ! vector opt.
+                  zetot(ji,jj,jk) = zetot(ji,jj,jk+1) + qsr_hc(ji,jj,jk) / r1_rau0_rcp
+               END DO
+            END DO
          END DO         
+!$OMP END PARALLEL
          CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation
          !
@@ -299,5 +353,12 @@
       !
       IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdt(ji,jj,jk) = tsa(ji,jj,jk,jp_tem) - ztrdt(ji,jj,jk)
+               END DO
+            END DO
+         END DO
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt )
          CALL wrk_dealloc( jpi,jpj,jpk,   ztrdt ) 
@@ -426,5 +487,12 @@
       END SELECT
       !
-      qsr_hc(:,:,:) = 0._wp     ! now qsr heat content set to zero where it will not be computed
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               qsr_hc(ji,jj,jk) = 0._wp     ! now qsr heat content set to zero where it will not be computed
+            END DO
+         END DO
+      END DO
       !
       ! 1st ocean level attenuation coefficient (used in sbcssm)
@@ -432,5 +500,10 @@
          CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev'  , fraqsr_1lev  )
       ELSE
-         fraqsr_1lev(:,:) = 1._wp   ! default : no penetration
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               fraqsr_1lev(ji,jj) = 1._wp   ! default : no penetration
+            END DO
+         END DO
       ENDIF
       !
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 7698)
@@ -88,12 +88,24 @@
       IF( l_trdtra ) THEN                    !* Save ta and sa trends
          CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdt(ji,jj,jk) = tsa(ji,jj,jk,jp_tem)
+                  ztrds(ji,jj,jk) = tsa(ji,jj,jk,jp_sal)
+               END DO
+            END DO
+         END DO
       ENDIF
       !
 !!gm  This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist)
       IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
-         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns
-         qsr(:,:) = 0._wp                     ! qsr set to zero
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               qns(ji,jj) = qns(ji,jj) + qsr(ji,jj)      ! total heat flux in qns
+               qsr(ji,jj) = 0._wp                     ! qsr set to zero
+            END DO
+         END DO
       ENDIF
 
@@ -111,12 +123,27 @@
          ELSE                                   ! No restart or restart not found: Euler forward time stepping
             zfact = 1._wp
-            sbc_tsc(:,:,:) = 0._wp
-            sbc_tsc_b(:,:,:) = 0._wp
+            DO jn = 1, jpts
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     sbc_tsc(ji,jj,jn) = 0._wp
+                     sbc_tsc_b(ji,jj,jn) = 0._wp
+                  END DO
+               END DO
+            END DO
          ENDIF
       ELSE                                !* other time-steps: swap of forcing fields
          zfact = 0.5_wp
-         sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:)
+         DO jn = 1, jpts
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  sbc_tsc_b(ji,jj,jn) = sbc_tsc(ji,jj,jn)
+               END DO
+            END DO
+         END DO
       ENDIF
       !                             !==  Now sbc tracer content fields  ==!
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
       DO jj = 2, jpj
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -126,4 +153,5 @@
       END DO
       IF( ln_linssh ) THEN                !* linear free surface  
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
          DO jj = 2, jpj                         !==>> add concentration/dilution effect due to constant volume cell
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -137,4 +165,5 @@
       !
       DO jn = 1, jpts               !==  update tracer trend  ==!
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
          DO jj = 2, jpj
             DO ji = fs_2, fs_jpim1   ! vector opt.  
@@ -218,4 +247,5 @@
       !
       IF( ln_iscpl .AND. ln_hsb) THEN         ! input of heat and salt due to river runoff 
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zdep)
          DO jk = 1,jpk
             DO jj = 2, jpj 
@@ -232,6 +262,13 @@
 
       IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdt(ji,jj,jk) = tsa(ji,jj,jk,jp_tem) - ztrdt(ji,jj,jk)
+                  ztrds(ji,jj,jk) = tsa(ji,jj,jk,jp_sal) - ztrds(ji,jj,jk)
+               END DO  
+            END DO  
+         END DO
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds )
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 7698)
@@ -58,5 +58,5 @@
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
       !
-      INTEGER  ::   jk                   ! Dummy loop indices
+      INTEGER  ::   jk, jj, ji           ! Dummy loop indices
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdt, ztrds   ! 3D workspace
       !!---------------------------------------------------------------------
@@ -72,6 +72,13 @@
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
          CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdt(ji,jj,jk) = tsa(ji,jj,jk,jp_tem)
+                  ztrds(ji,jj,jk) = tsa(ji,jj,jk,jp_sal)
+               END DO
+            END DO
+         END DO
       ENDIF
       !
@@ -84,11 +91,23 @@
       ! JMM avoid negative salinities near river outlet ! Ugly fix
       ! JMM : restore negative salinities to small salinities:
-      WHERE( tsa(:,:,:,jp_sal) < 0._wp )   tsa(:,:,:,jp_sal) = 0.1_wp
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( tsa(ji,jj,jk,jp_sal) < 0._wp )   tsa(ji,jj,jk,jp_sal) = 0.1_wp
+            END DO
+         END DO
+      END DO
 !!gm
 
       IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
          DO jk = 1, jpkm1
-            ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem) - tsb(:,:,jk,jp_tem) ) / r2dt ) - ztrdt(:,:,jk)
-            ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal) - tsb(:,:,jk,jp_sal) ) / r2dt ) - ztrds(:,:,jk)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdt(ji,jj,jk) = ( ( tsa(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) / r2dt ) - ztrdt(ji,jj,jk)
+                  ztrds(ji,jj,jk) = ( ( tsa(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) / r2dt ) - ztrds(ji,jj,jk)
+               END DO
+            END DO
          END DO
 !!gm this should be moved in trdtra.F90 and done on all trends
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 7698)
@@ -106,11 +106,29 @@
             !
             ! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
-            IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN   ;   zwt(:,:,2:jpk) = avt  (:,:,2:jpk)
-            ELSE                                            ;   zwt(:,:,2:jpk) = fsavs(:,:,2:jpk)
+            IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zwt(ji,jj,2:jpk) = avt  (ji,jj,2:jpk)
+                  END DO
+               END DO
+            ELSE                                            
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zwt(ji,jj,2:jpk) = fsavs(ji,jj,2:jpk)
+                  END DO
+               END DO
             ENDIF
-            zwt(:,:,1) = 0._wp
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zwt(ji,jj,1) = 0._wp
+               END DO
+            END DO
             !
             IF( l_ldfslp ) THEN            ! isoneutral diffusion: add the contribution 
                IF( ln_traldf_msc  ) THEN     ! MSC iso-neutral operator 
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
                   DO jk = 2, jpkm1
                      DO jj = 2, jpjm1
@@ -121,4 +139,5 @@
                   END DO
                ELSE                          ! standard or triad iso-neutral operator
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
                   DO jk = 2, jpkm1
                      DO jj = 2, jpjm1
@@ -132,4 +151,6 @@
             !
             ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked)
+!$OMP PARALLEL 
+!$OMP DO schedule(static) private(jk, jj, ji)
             DO jk = 1, jpkm1
                DO jj = 2, jpjm1
@@ -162,4 +183,5 @@
             !   used as a work space array: its value is modified.
             !
+!$OMP DO schedule(static) private(jj, ji)
             DO jj = 2, jpjm1        !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
                DO ji = fs_2, fs_jpim1            ! done one for all passive tracers (so included in the IF instruction)
@@ -167,5 +189,7 @@
                END DO
             END DO
+!$OMP END DO NOWAIT 
             DO jk = 2, jpkm1
+!$OMP DO schedule(static) private(jj, ji)
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1
@@ -174,7 +198,10 @@
                END DO
             END DO
+!$OMP END PARALLEL 
             !
          ENDIF 
          !         
+!$OMP PARALLEL 
+!$OMP DO schedule(static) private(jj, ji)
          DO jj = 2, jpjm1           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
             DO ji = fs_2, fs_jpim1
@@ -183,4 +210,5 @@
          END DO
          DO jk = 2, jpkm1
+!$OMP DO schedule(static) private(jj, ji, zrhs)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
@@ -191,4 +219,5 @@
          END DO
          !
+!$OMP DO schedule(static) private(jj, ji)
          DO jj = 2, jpjm1           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer)
             DO ji = fs_2, fs_jpim1
@@ -197,4 +226,5 @@
          END DO
          DO jk = jpk-2, 1, -1
+!$OMP DO schedule(static) private(jj, ji)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
@@ -204,4 +234,5 @@
             END DO
          END DO
+!$OMP END PARALLEL 
          !                                            ! ================= !
       END DO                                          !  end tracer loop  !
Index: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 7646)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 7698)
@@ -101,9 +101,24 @@
       IF( nn_timing == 1 )   CALL timing_start( 'zps_hde')
       !
-      pgtu(:,:,:)=0._wp   ;   zti (:,:,:)=0._wp   ;   zhi (:,:  )=0._wp
-      pgtv(:,:,:)=0._wp   ;   ztj (:,:,:)=0._wp   ;   zhj (:,:  )=0._wp
+      DO jn = 1, kjpt
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               pgtu(ji,jj,jn)=0._wp   ;   zti (ji,jj,jn)=0._wp
+               pgtv(ji,jj,jn)=0._wp   ;   ztj (ji,jj,jn)=0._wp
+            END DO
+         END DO
+      END DO
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1
+            zhi (ji,jj  )=0._wp
+            zhj (ji,jj  )=0._wp
+         END DO
+       END DO
       !
       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
          !
+!$OMP PARALLEL DO schedule(static) private(jj,ji,iku,ikv,ze3wu,ze3wv,zmaxu,zmaxv)
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
@@ -150,6 +165,14 @@
       !                
       IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
-         pgru(:,:) = 0._wp
-         pgrv(:,:) = 0._wp                ! depth of the partial step level
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj,ji)
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               pgru(ji,jj) = 0._wp
+               pgrv(ji,jj) = 0._wp                ! depth of the partial step level
+            END DO
+         END DO
+!$OMP END DO NOWAIT
+!$OMP DO schedule(static) private(jj,ji,iku,ikv,ze3wu,ze3wv)
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
@@ -166,8 +189,11 @@
             END DO
          END DO
+!$OMP END DO NOWAIT
+!$OMP END PARALLEL
          !
          CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj 
          CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj 
          !
+!$OMP PARALLEL DO schedule(static) private(jj,ji,iku,ikv,ze3wu,ze3wv)
          DO jj = 1, jpjm1                 ! Gradient of density at the last level 
             DO ji = 1, jpim1
