Index: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/depth_e3.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DOM/depth_e3.F90	(revision 7698)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DOM/depth_e3.F90	(revision 7753)
@@ -150,25 +150,13 @@
       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdept_3d, pdepw_3d   ! depth = SUM( e3 )     [m]
       !
-      INTEGER  ::   jk, jj, ji           ! dummy loop indices
+      INTEGER  ::   jk           ! dummy loop indices
       !!----------------------------------------------------------------------      
       !
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            pdepw_3d(ji,jj,1) = 0.0_wp
-            pdept_3d(ji,jj,1) = 0.5_wp * pe3w_3d(ji,jj,1)
-         END DO
+      pdepw_3d(:,:,1) = 0.0_wp
+      pdept_3d(:,:,1) = 0.5_wp * pe3w_3d(:,:,1)
+      DO jk = 2, jpk
+         pdepw_3d(:,:,jk) = pdepw_3d(:,:,jk-1) + pe3t_3d(:,:,jk-1) 
+         pdept_3d(:,:,jk) = pdept_3d(:,:,jk-1) + pe3w_3d(:,:,jk  ) 
       END DO
-      DO jk = 2, jpk
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               pdepw_3d(ji,jj,jk) = pdepw_3d(ji,jj,jk-1) + pe3t_3d(ji,jj,jk-1) 
-               pdept_3d(ji,jj,jk) = pdept_3d(ji,jj,jk-1) + pe3w_3d(ji,jj,jk  ) 
-            END DO
-         END DO
-      END DO
-!$OMP END PARALLEL
       !
    END SUBROUTINE e3_to_depth_3d
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 7698)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 7753)
@@ -133,6 +133,4 @@
       CALL dom_msk( ik_top, ik_bot )   ! Masks
       !
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji,ik)
       DO jj = 1, jpj                   ! depth of the iceshelves
          DO ji = 1, jpi
@@ -142,24 +140,12 @@
       END DO
       !
-!$OMP END DO NOWAIT
-!$OMP DO schedule(static) private(jj,ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ht_0(ji,jj) = 0._wp  ! Reference ocean thickness
-            hu_0(ji,jj) = 0._wp
-            hv_0(ji,jj) = 0._wp
-         END DO
+      ht_0(:,:) = 0._wp  ! Reference ocean thickness
+      hu_0(:,:) = 0._wp
+      hv_0(:,:) = 0._wp
+      DO jk = 1, jpk
+         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
+         hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk)
+         hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
       END DO
-      DO jk = 1, jpk
-!$OMP DO schedule(static) private(jj,ji,ik)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ht_0(ji,jj) = ht_0(ji,jj) + e3t_0(ji,jj,jk) * tmask(ji,jj,jk)
-               hu_0(ji,jj) = hu_0(ji,jj) + e3u_0(ji,jj,jk) * umask(ji,jj,jk)
-               hv_0(ji,jj) = hv_0(ji,jj) + e3v_0(ji,jj,jk) * vmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-!$OMP END PARALLEL
       !
       !           !==  time varying part of coordinate system  ==!
@@ -180,11 +166,6 @@
              e3vw_b =  e3vw_0  ;    e3vw_n =  e3vw_0   !        ---          !
          !
-!$OMP PARALLEL DO schedule(static) private(jj,ji)
-         DO jj =1, jpj
-            DO ji=1, jpi
-               z1_hu_0(ji,jj) = ssumask(ji,jj) / ( hu_0(ji,jj) + 1._wp - ssumask(ji,jj) )     ! _i mask due to ISF
-               z1_hv_0(ji,jj) = ssvmask(ji,jj) / ( hv_0(ji,jj) + 1._wp - ssvmask(ji,jj) )
-            END DO
-         END DO
+         z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )     ! _i mask due to ISF
+         z1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) )
          !
          !        before       !          now          !       after         !
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 7698)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 7753)
@@ -40,5 +40,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.7 , NEMO Consortium (2016)
-   !! $Id$
+   !! $Id$ 
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -117,11 +117,6 @@
       IF( iff == 0 ) THEN                 ! Coriolis parameter has not been defined 
          IF(lwp) WRITE(numout,*) '          Coriolis parameter calculated on the sphere from gphif & gphit'
-!$OMP PARALLEL DO schedule(static) private(jj, ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ff_f(ji,jj) = 2. * omega * SIN( rad * gphif(ji,jj) )     ! compute it on the sphere at f-point
-               ff_t(ji,jj) = 2. * omega * SIN( rad * gphit(ji,jj) )     !    -        -       -    at t-point
-            END DO
-         END DO
+         ff_f(:,:) = 2. * omega * SIN( rad * gphif(:,:) )     ! compute it on the sphere at f-point
+         ff_t(:,:) = 2. * omega * SIN( rad * gphit(:,:) )     !    -        -       -    at t-point
       ELSE
          IF( ln_read_cfg ) THEN
@@ -135,39 +130,24 @@
       !                             !==  associated horizontal metrics  ==!
       !
-!$OMP PARALLEL DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            r1_e1t(ji,jj) = 1._wp / e1t(ji,jj)   ;   r1_e2t (ji,jj) = 1._wp / e2t(ji,jj)
-            r1_e1u(ji,jj) = 1._wp / e1u(ji,jj)   ;   r1_e2u (ji,jj) = 1._wp / e2u(ji,jj)
-            r1_e1v(ji,jj) = 1._wp / e1v(ji,jj)   ;   r1_e2v (ji,jj) = 1._wp / e2v(ji,jj)
-            r1_e1f(ji,jj) = 1._wp / e1f(ji,jj)   ;   r1_e2f (ji,jj) = 1._wp / e2f(ji,jj)
-            !
-            e1e2t (ji,jj) = e1t(ji,jj) * e2t(ji,jj)   ;   r1_e1e2t(ji,jj) = 1._wp / e1e2t(ji,jj)
-            e1e2f (ji,jj) = e1f(ji,jj) * e2f(ji,jj)   ;   r1_e1e2f(ji,jj) = 1._wp / e1e2f(ji,jj)
-         END DO
-      END DO
+      r1_e1t(:,:) = 1._wp / e1t(:,:)   ;   r1_e2t (:,:) = 1._wp / e2t(:,:)
+      r1_e1u(:,:) = 1._wp / e1u(:,:)   ;   r1_e2u (:,:) = 1._wp / e2u(:,:)
+      r1_e1v(:,:) = 1._wp / e1v(:,:)   ;   r1_e2v (:,:) = 1._wp / e2v(:,:)
+      r1_e1f(:,:) = 1._wp / e1f(:,:)   ;   r1_e2f (:,:) = 1._wp / e2f(:,:)
+      !
+      e1e2t (:,:) = e1t(:,:) * e2t(:,:)   ;   r1_e1e2t(:,:) = 1._wp / e1e2t(:,:)
+      e1e2f (:,:) = e1f(:,:) * e2f(:,:)   ;   r1_e1e2f(:,:) = 1._wp / e1e2f(:,:)
       IF( ie1e2u_v == 0 ) THEN               ! u- & v-surfaces have not been defined
          IF(lwp) WRITE(numout,*) '          u- & v-surfaces calculated as e1 e2 product'
-!$OMP PARALLEL DO schedule(static) private(jj, ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               e1e2u (ji,jj) = e1u(ji,jj) * e2u(ji,jj)         ! compute them
-               e1e2v (ji,jj) = e1v(ji,jj) * e2v(ji,jj) 
-            END DO
-         END DO
+         e1e2u (:,:) = e1u(:,:) * e2u(:,:)         ! compute them
+         e1e2v (:,:) = e1v(:,:) * e2v(:,:) 
       ELSE
          IF(lwp) WRITE(numout,*) '          u- & v-surfaces have been read in "mesh_mask" file:'
          IF(lwp) WRITE(numout,*) '                     grid size reduction in strait(s) is used'
       ENDIF
-!$OMP PARALLEL DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            r1_e1e2u(ji,jj) = 1._wp / e1e2u(ji,jj)     ! compute their invert in any cases
-            r1_e1e2v(ji,jj) = 1._wp / e1e2v(ji,jj)
-            !   
-            e2_e1u(ji,jj) = e2u(ji,jj) / e1u(ji,jj)
-            e1_e2v(ji,jj) = e1v(ji,jj) / e2v(ji,jj)
-         END DO
-      END DO
+      r1_e1e2u(:,:) = 1._wp / e1e2u(:,:)     ! compute their invert in any cases
+      r1_e1e2v(:,:) = 1._wp / e1e2v(:,:)
+      !   
+      e2_e1u(:,:) = e2u(:,:) / e1u(:,:)
+      e1_e2v(:,:) = e1v(:,:) / e2v(:,:)
       !
       !
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 7698)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 7753)
@@ -47,5 +47,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.2 , LODYC-IPSL  (2009)
-   !! $Id$
+   !! $Id$ 
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -137,14 +137,5 @@
       ! ----------------------------
       !
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jk, jj, ji)
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               tmask(ji,jj,jk) = 0._wp
-            END DO
-         END DO
-      END DO
-!$OMP DO schedule(static) private(jj, ji, iktop, ikbot)
+      tmask(:,:,:) = 0._wp
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -156,5 +147,4 @@
          END DO  
       END DO  
-!$OMP END PARALLEL
 !SF  add here lbc_lnk: bug not still understood : cause now domain configuration is read !
 !!gm I don't understand why...  
@@ -171,5 +161,4 @@
       ! ------------------------
       IF ( ln_bdy .AND. ln_mask_file ) THEN
-!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -184,5 +173,4 @@
       ! ----------------------------------------
       ! NB: at this point, fmask is designed for free slip lateral boundary condition
-!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
       DO jk = 1, jpk
          DO jj = 1, jpjm1
@@ -204,24 +192,12 @@
       ! Ocean/land mask at wu-, wv- and w points    (computed from tmask)
       !-----------------------------------------
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            wmask (ji,jj,1) = tmask(ji,jj,1)     ! surface
-            wumask(ji,jj,1) = umask(ji,jj,1)
-            wvmask(ji,jj,1) = vmask(ji,jj,1)
-         END DO
+      wmask (:,:,1) = tmask(:,:,1)     ! surface
+      wumask(:,:,1) = umask(:,:,1)
+      wvmask(:,:,1) = vmask(:,:,1)
+      DO jk = 2, jpk                   ! interior values
+         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
+         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)   
+         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1)
       END DO
-!$OMP DO schedule(static) private(jk,jj,ji)
-      DO jk = 2, jpk                   ! interior values
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               wmask (ji,jj,jk) = tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
-               wumask(ji,jj,jk) = umask(ji,jj,jk) * umask(ji,jj,jk-1)   
-               wvmask(ji,jj,jk) = vmask(ji,jj,jk) * vmask(ji,jj,jk-1)
-            END DO
-         END DO
-      END DO
-!$OMP END PARALLEL
 
 
@@ -240,10 +216,5 @@
       !
       !                          ! halo mask : 0 on the halo and 1 elsewhere
-!$OMP PARALLEL DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            tmask_h(ji,jj) = 1._wp                  
-         END DO
-      END DO
+      tmask_h(:,:) = 1._wp                  
       tmask_h( 1 :iif,   :   ) = 0._wp      ! first columns
       tmask_h(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns)
@@ -270,10 +241,5 @@
       !
       !                          ! interior mask : 2D ocean mask x halo mask 
-!$OMP PARALLEL DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            tmask_i(ji,jj) = ssmask(ji,jj) * tmask_h(ji,jj)
-         END DO
-      END DO
+      tmask_i(:,:) = ssmask(:,:) * tmask_h(:,:)
 
 
@@ -284,13 +250,6 @@
          CALL wrk_alloc( jpi,jpj,   zwf )
          !
-!$OMP PARALLEL
          DO jk = 1, jpk
-!$OMP DO schedule(static) private(jj, ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zwf(ji,jj) = fmask(ji,jj,jk)         
-               END DO
-            END DO
-!$OMP DO schedule(static) private(jj, ji)
+            zwf(:,:) = fmask(:,:,jk)         
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -301,5 +260,4 @@
                END DO
             END DO
-!$OMP DO schedule(static) private(jj)
             DO jj = 2, jpjm1
                IF( fmask(1,jj,jk) == 0._wp ) THEN
@@ -310,5 +268,4 @@
                ENDIF
             END DO         
-!$OMP DO schedule(static) private(ji)
             DO ji = 2, jpim1
                IF( fmask(ji,1,jk) == 0._wp ) THEN
@@ -320,5 +277,4 @@
             END DO
          END DO
-!$OMP END PARALLEL
          !
          CALL wrk_dealloc( jpi,jpj,   zwf )
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 7698)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 7753)
@@ -135,10 +135,5 @@
       !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
       CALL dom_vvl_rst( nit000, 'READ' )
-!$OMP PARALLEL DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            e3t_a(ji,jj,jpk) = e3t_0(ji,jj,jpk)  ! last level always inside the sea floor set one for all
-         END DO
-      END DO
+      e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
       !
       !                    !== Set of all other vertical scale factors  ==!  (now and before)
@@ -158,17 +153,10 @@
       !
       !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            gdept_n(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1)       ! reference to the ocean surface (used for MLD and light penetration)
-            gdepw_n(ji,jj,1) = 0.0_wp
-            gde3w_n(ji,jj,1) = gdept_n(ji,jj,1) - sshn(ji,jj)  ! reference to a common level z=0 for hpg
-            gdept_b(ji,jj,1) = 0.5_wp * e3w_b(ji,jj,1)
-            gdepw_b(ji,jj,1) = 0.0_wp
-         END DO
-      END DO
+      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)       ! reference to the ocean surface (used for MLD and light penetration)
+      gdepw_n(:,:,1) = 0.0_wp
+      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)  ! reference to a common level z=0 for hpg
+      gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1)
+      gdepw_b(:,:,1) = 0.0_wp
       DO jk = 2, jpk                               ! vertical sum
-!$OMP DO schedule(static) private(jj,ji,zcoef)
          DO jj = 1,jpj
             DO ji = 1,jpi
@@ -190,38 +178,23 @@
       !
       !                    !==  thickness of the water column  !!   (ocean portion only)
-!$OMP DO schedule(static) private(jj,ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ht_n(ji,jj) = e3t_n(ji,jj,1) * tmask(ji,jj,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
-            hu_b(ji,jj) = e3u_b(ji,jj,1) * umask(ji,jj,1)
-            hu_n(ji,jj) = e3u_n(ji,jj,1) * umask(ji,jj,1)
-            hv_b(ji,jj) = e3v_b(ji,jj,1) * vmask(ji,jj,1)
-            hv_n(ji,jj) = e3v_n(ji,jj,1) * vmask(ji,jj,1)
-         END DO
+      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
+      hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
+      hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1)
+      hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
+      hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1)
+      DO jk = 2, jpkm1
+         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
+         hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
+         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
+         hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
+         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
       END DO
-      DO jk = 2, jpkm1
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ht_n(ji,jj) = ht_n(ji,jj) + e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
-               hu_b(ji,jj) = hu_b(ji,jj) + e3u_b(ji,jj,jk) * umask(ji,jj,jk)
-               hu_n(ji,jj) = hu_n(ji,jj) + e3u_n(ji,jj,jk) * umask(ji,jj,jk)
-               hv_b(ji,jj) = hv_b(ji,jj) + e3v_b(ji,jj,jk) * vmask(ji,jj,jk)
-               hv_n(ji,jj) = hv_n(ji,jj) + e3v_n(ji,jj,jk) * vmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
       !
       !                    !==  inverse of water column thickness   ==!   (u- and v- points)
-!$OMP DO schedule(static) private(jj,ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            r1_hu_b(ji,jj) = ssumask(ji,jj) / ( hu_b(ji,jj) + 1._wp - ssumask(ji,jj) )    ! _i mask due to ISF
-            r1_hu_n(ji,jj) = ssumask(ji,jj) / ( hu_n(ji,jj) + 1._wp - ssumask(ji,jj) )
-            r1_hv_b(ji,jj) = ssvmask(ji,jj) / ( hv_b(ji,jj) + 1._wp - ssvmask(ji,jj) )
-            r1_hv_n(ji,jj) = ssvmask(ji,jj) / ( hv_n(ji,jj) + 1._wp - ssvmask(ji,jj) )
-         END DO
-      END DO
-!$OMP END PARALLEL
+      r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
+      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
+      r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
+      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
+
       !                    !==   z_tilde coordinate case  ==!   (Restoring frequencies)
       IF( ln_vvl_ztilde ) THEN
@@ -229,23 +202,12 @@
          !                                   ! Values in days provided via the namelist
          !                                   ! use rsmall to avoid possible division by zero errors with faulty settings
-!$OMP PARALLEL DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               frq_rst_e3t(ji,jj) = 2._wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp )
-               frq_rst_hdv(ji,jj) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
-            END DO
-         END DO
+         frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp )
+         frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
          !
          IF( ln_vvl_ztilde_as_zstar ) THEN   ! z-star emulation using z-tile
-!$OMP PARALLEL DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  frq_rst_e3t(ji,jj) = 0._wp               !Ignore namelist settings
-                  frq_rst_hdv(ji,jj) = 1._wp / rdt
-               END DO
-            END DO
+            frq_rst_e3t(:,:) = 0._wp               !Ignore namelist settings
+            frq_rst_hdv(:,:) = 1._wp / rdt
          ENDIF
          IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator
-!$OMP PARALLEL DO schedule(static) private(jj,ji)
             DO jj = 1, jpj
                DO ji = 1, jpi
@@ -343,21 +305,9 @@
       !                                                ! --------------------------------------------- !
       !
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            z_scale(ji,jj) = ( ssha(ji,jj) - sshb(ji,jj) ) * ssmask(ji,jj) / ( ht_0(ji,jj) + sshn(ji,jj) + 1. - ssmask(ji,jj) )
-         END DO
+      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
+      DO jk = 1, jpkm1
+         ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0)
+         e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
       END DO
-!$OMP DO schedule(static) private(jk,jj,ji)
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0)
-               e3t_a(ji,jj,jk) = e3t_b(ji,jj,jk) + e3t_n(ji,jj,jk) * z_scale(ji,jj) * tmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-!$OMP END PARALLEL
       !
       IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
@@ -368,28 +318,11 @@
          ! 1 - barotropic divergence
          ! -------------------------
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zhdiv(ji,jj) = 0._wp
-               zht(ji,jj)   = 0._wp
-            END DO
-         END DO
+         zhdiv(:,:) = 0._wp
+         zht(:,:)   = 0._wp
          DO jk = 1, jpkm1
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zhdiv(ji,jj) = zhdiv(ji,jj) + e3t_n(ji,jj,jk) * hdivn(ji,jj,jk)
-                  zht  (ji,jj) = zht  (ji,jj) + e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zhdiv(ji,jj) = zhdiv(ji,jj) / ( zht(ji,jj) + 1. - tmask_i(ji,jj) )
-            END DO
-         END DO
-!$OMP END PARALLEL
+            zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
+            zht  (:,:) = zht  (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
+         END DO
+         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
 
          ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
@@ -397,12 +330,7 @@
          IF( ln_vvl_ztilde ) THEN
             IF( kt > nit000 ) THEN
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
                DO jk = 1, jpkm1
-                  DO jj = 1, jpj
-                     DO ji = 1, jpi
-                        hdiv_lf(ji,jj,jk) = hdiv_lf(ji,jj,jk) - rdt * frq_rst_hdv(ji,jj)   &
-                           &          * ( hdiv_lf(ji,jj,jk) - e3t_n(ji,jj,jk) * ( hdivn(ji,jj,jk) - zhdiv(ji,jj) ) )
-                     END DO
-                  END DO
+                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
+                     &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
                END DO
             ENDIF
@@ -411,32 +339,15 @@
          ! II - after z_tilde increments of vertical scale factors
          ! =======================================================
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  tilde_e3t_a(ji,jj,jk) = 0._wp  ! tilde_e3t_a used to store tendency terms
-               END DO
-            END DO
-         END DO
+         tilde_e3t_a(:,:,:) = 0._wp  ! tilde_e3t_a used to store tendency terms
 
          ! 1 - High frequency divergence term
          ! ----------------------------------
          IF( ln_vvl_ztilde ) THEN     ! z_tilde case
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
             DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) - ( e3t_n(ji,jj,jk) * ( hdivn(ji,jj,jk) - zhdiv(ji,jj) ) - hdiv_lf(ji,jj,jk) )
-                  END DO
-               END DO
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
             END DO
          ELSE                         ! layer case
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
             DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) -   e3t_n(ji,jj,jk) * ( hdivn(ji,jj,jk) - zhdiv(ji,jj) ) * tmask(ji,jj,jk)
-                  END DO
-               END DO
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
             END DO
          ENDIF
@@ -445,11 +356,6 @@
          ! ------------------
          IF( ln_vvl_ztilde ) THEN
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
             DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) - frq_rst_e3t(ji,jj) * tilde_e3t_b(ji,jj,jk)
-                  END DO
-               END DO
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
             END DO
          ENDIF
@@ -457,14 +363,7 @@
          ! 3 - Thickness diffusion term
          ! ----------------------------
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zwu(ji,jj) = 0._wp
-               zwv(ji,jj) = 0._wp
-            END DO
-         END DO
+         zwu(:,:) = 0._wp
+         zwv(:,:) = 0._wp
          DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes
-!$OMP DO schedule(static) private(jj,ji)
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
@@ -478,5 +377,4 @@
             END DO
          END DO
-!$OMP DO schedule(static) private(jj,ji)
          DO jj = 1, jpj          ! b - correction for last oceanic u-v points
             DO ji = 1, jpi
@@ -485,5 +383,4 @@
             END DO
          END DO
-!$OMP DO schedule(static) private(jk,jj,ji)
          DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes
             DO jj = 2, jpjm1
@@ -495,5 +392,4 @@
             END DO
          END DO
-!$OMP END PARALLEL
          !                       ! d - thickness diffusion transport: boundary conditions
          !                             (stored for tracer advction and continuity equation)
@@ -511,31 +407,12 @@
          ENDIF
          CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp )
-!$OMP PARALLEL 
-!$OMP DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_b(ji,jj,jk) + z2dt * tmask(ji,jj,jk) * tilde_e3t_a(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
 
          ! Maximum deformation control
          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ze3t(ji,jj,jpk) = 0._wp
-            END DO
-         END DO
-!$OMP DO schedule(static) private(jk,jj,ji)
+         ze3t(:,:,jpk) = 0._wp
          DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               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)
-               END DO
-            END DO
-         END DO
-!$OMP END PARALLEL
+            ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
+         END DO
          z_tmax = MAXVAL( ze3t(:,:,:) )
          IF( lk_mpp )   CALL mpp_max( z_tmax )                 ! max over the global domain
@@ -565,25 +442,12 @@
          ! - ML - end test
          ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  tilde_e3t_a(ji,jj,jk) = MIN( tilde_e3t_a(ji,jj,jk),   rn_zdef_max * e3t_0(ji,jj,jk) )
-                  tilde_e3t_a(ji,jj,jk) = MAX( tilde_e3t_a(ji,jj,jk), - rn_zdef_max * e3t_0(ji,jj,jk) )
-               END DO
-            END DO
-         END DO
+         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
+         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
 
          !
          ! "tilda" change in the after scale factor
          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!$OMP DO schedule(static) private(jk,jj,ji)
          DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  dtilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) - tilde_e3t_b(ji,jj,jk)
-               END DO
-            END DO
+            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
          END DO
          ! III - Barotropic repartition of the sea surface height over the baroclinic profile
@@ -593,44 +457,19 @@
          !        i.e. locally and not spread over the water column.
          !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zht(ji,jj) = 0.
-            END DO
-         END DO
+         zht(:,:) = 0.
          DO jk = 1, jpkm1
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zht(ji,jj)  = zht(ji,jj) + tilde_e3t_a(ji,jj,jk) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               z_scale(ji,jj) =  - zht(ji,jj) / ( ht_0(ji,jj) + sshn(ji,jj) + 1. - ssmask(ji,jj) )
-            END DO
-         END DO
-!$OMP DO schedule(static) private(jk,jj,ji)
+            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
+         END DO
+         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
          DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  dtilde_e3t_a(ji,jj,jk) = dtilde_e3t_a(ji,jj,jk) + e3t_n(ji,jj,jk) * z_scale(ji,jj) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-!$OMP END PARALLEL
+            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+         END DO
+
       ENDIF
 
       IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
       !                                           ! ---baroclinic part--------- !
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
          DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  e3t_a(ji,jj,jk) = e3t_a(ji,jj,jk) + dtilde_e3t_a(ji,jj,jk) * tmask(ji,jj,jk)
-               END DO
-            END DO
+            e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
       ENDIF
@@ -645,60 +484,24 @@
          END IF
          !
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zht(ji,jj) = 0.0_wp
-            END DO
-         END DO
+         zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zht(ji,jj) = zht(ji,jj) + e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-!$OMP END PARALLEL
+            zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
+         END DO
          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax
          !
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zht(ji,jj) = 0.0_wp
-            END DO
-         END DO
+         zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zht(ji,jj) = zht(ji,jj) + e3t_a(ji,jj,jk) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-!$OMP END PARALLEL
+            zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk)
+         END DO
          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax
          !
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zht(ji,jj) = 0.0_wp
-            END DO
-         END DO
+         zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-!$OMP DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zht(ji,jj) = zht(ji,jj) + e3t_b(ji,jj,jk) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-!$OMP END PARALLEL
+            zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk)
+         END DO
          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
@@ -729,31 +532,14 @@
       ! *********************************** !
 
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            hu_a(ji,jj) = e3u_a(ji,jj,1) * umask(ji,jj,1)
-            hv_a(ji,jj) = e3v_a(ji,jj,1) * vmask(ji,jj,1)
-         END DO
-      END DO
+      hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1)
+      hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)
       DO jk = 2, jpkm1
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               hu_a(ji,jj) = hu_a(ji,jj) + e3u_a(ji,jj,jk) * umask(ji,jj,jk)
-               hv_a(ji,jj) = hv_a(ji,jj) + e3v_a(ji,jj,jk) * vmask(ji,jj,jk)
-            END DO
-         END DO
+         hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)
+         hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)
       END DO
       !                                        ! Inverse of the local depth
 !!gm BUG ?  don't understand the use of umask_i here .....
-!$OMP DO schedule(static) private(jj,ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            r1_hu_a(ji,jj) = ssumask(ji,jj) / ( hu_a(ji,jj) + 1._wp - ssumask(ji,jj) )
-            r1_hv_a(ji,jj) = ssvmask(ji,jj) / ( hv_a(ji,jj) + 1._wp - ssvmask(ji,jj) )
-         END DO
-      END DO
-!$OMP END PARALLEL
+      r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )
+      r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )
       !
       CALL wrk_dealloc( jpi,jpj,       zht, z_scale, zwu, zwv, zhdiv )
@@ -810,45 +596,17 @@
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
          IF( neuler == 0 .AND. kt == nit000 ) THEN
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     tilde_e3t_b(ji,jj,jk) = tilde_e3t_n(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
+            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
          ELSE
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     tilde_e3t_b(ji,jj,jk) = tilde_e3t_n(ji,jj,jk) & 
-                     &         + atfp * ( tilde_e3t_b(ji,jj,jk) - 2.0_wp * tilde_e3t_n(ji,jj,jk) + tilde_e3t_a(ji,jj,jk) )
-                  END DO
-               END DO
-            END DO
+            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
+            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
          ENDIF
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  tilde_e3t_n(ji,jj,jk) = tilde_e3t_a(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-      ENDIF
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               gdept_b(ji,jj,jk) = gdept_n(ji,jj,jk)
-               gdepw_b(ji,jj,jk) = gdepw_n(ji,jj,jk)
-        
-               e3t_n(ji,jj,jk) = e3t_a(ji,jj,jk)
-               e3u_n(ji,jj,jk) = e3u_a(ji,jj,jk)
-               e3v_n(ji,jj,jk) = e3v_a(ji,jj,jk)
-            END DO
-         END DO
-      END DO
+         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
+      ENDIF
+      gdept_b(:,:,:) = gdept_n(:,:,:)
+      gdepw_b(:,:,:) = gdepw_n(:,:,:)
+
+      e3t_n(:,:,:) = e3t_a(:,:,:)
+      e3u_n(:,:,:) = e3u_a(:,:,:)
+      e3v_n(:,:,:) = e3v_a(:,:,:)
 
       ! Compute all missing vertical scale factor and depths
@@ -870,15 +628,8 @@
 
       ! t- and w- points depth (set the isf depth as it is in the initial step)
-! !$OMP PARALLEL
-! !$OMP DO schedule(static) private(jj,ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            gdept_n(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1)
-            gdepw_n(ji,jj,1) = 0.0_wp
-            gde3w_n(ji,jj,1) = gdept_n(ji,jj,1) - sshn(ji,jj)
-         END DO
-      END DO
+      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
+      gdepw_n(:,:,1) = 0.0_wp
+      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
       DO jk = 2, jpk
-! !$OMP DO schedule(static) private(jj,ji,zcoef)
          DO jj = 1,jpj
             DO ji = 1,jpi
@@ -896,23 +647,12 @@
       ! Local depth and Inverse of the local depth of the water
       ! -------------------------------------------------------
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            hu_n(ji,jj) = hu_a(ji,jj)   ;   r1_hu_n(ji,jj) = r1_hu_a(ji,jj)
-            hv_n(ji,jj) = hv_a(ji,jj)   ;   r1_hv_n(ji,jj) = r1_hv_a(ji,jj)
-            !
-            ht_n(ji,jj) = e3t_n(ji,jj,1) * tmask(ji,jj,1)
-         END DO
+      hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:)
+      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
+      !
+      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)
+      DO jk = 2, jpkm1
+         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
       END DO
-      DO jk = 2, jpkm1
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ht_n(ji,jj) = ht_n(ji,jj) + e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-!$OMP END PARALLEL
+
       ! write restart file
       ! ==================
@@ -954,5 +694,4 @@
          !
       CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
          DO jk = 1, jpk
             DO jj = 1, jpjm1
@@ -965,15 +704,7 @@
          END DO
          CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  pe3_out(ji,jj,jk) = pe3_out(ji,jj,jk) + e3u_0(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
          !
       CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
          DO jk = 1, jpk
             DO jj = 1, jpjm1
@@ -986,15 +717,7 @@
          END DO
          CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  pe3_out(ji,jj,jk) = pe3_out(ji,jj,jk) + e3v_0(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
          !
       CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
          DO jk = 1, jpk
             DO jj = 1, jpjm1
@@ -1008,84 +731,41 @@
          END DO
          CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  pe3_out(ji,jj,jk) = pe3_out(ji,jj,jk) + e3f_0(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
          !
       CASE( 'W' )                   !* from T- to W-point : vertical simple mean
          !
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               pe3_out(ji,jj,1) = e3w_0(ji,jj,1) + pe3_in(ji,jj,1) - e3t_0(ji,jj,1)
-            END DO
-         END DO
+         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
          ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
 !!gm BUG? use here wmask in case of ISF ?  to be checked
-!$OMP DO schedule(static) private(jk,jj,ji)
          DO jk = 2, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  pe3_out(ji,jj,jk) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * ( tmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) )   &
-                     &                            * ( pe3_in(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )                               &
-                     &                            +            0.5_wp * ( tmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )     &
-                     &                            * ( pe3_in(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
-               END DO
-            END DO
-         END DO
-!$OMP END PARALLEL
+            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   &
+               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               &
+               &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     &
+               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
+         END DO
          !
       CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean
          !
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               pe3_out(ji,jj,1) = e3uw_0(ji,jj,1) + pe3_in(ji,jj,1) - e3u_0(ji,jj,1)
-            END DO
-         END DO
+         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
 !!gm BUG? use here wumask in case of ISF ?  to be checked
-!$OMP DO schedule(static) private(jk,jj,ji)
          DO jk = 2, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  pe3_out(ji,jj,jk) = e3uw_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
-                     &                             * ( pe3_in(ji,jj,jk-1) - e3u_0(ji,jj,jk-1) )                              &
-                     &                             +            0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )    &
-                     &                             * ( pe3_in(ji,jj,jk  ) - e3u_0(ji,jj,jk  ) )
-               END DO
-            END DO
-         END DO
-!$OMP END PARALLEL
+            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
+               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              &
+               &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
+               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
+         END DO
          !
       CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean
          !
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               pe3_out(ji,jj,1) = e3vw_0(ji,jj,1) + pe3_in(ji,jj,1) - e3v_0(ji,jj,1)
-            END DO
-         END DO
+         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
 !!gm BUG? use here wvmask in case of ISF ?  to be checked
-!$OMP DO schedule(static) private(jk,jj,ji)
          DO jk = 2, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  pe3_out(ji,jj,jk) = e3vw_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
-                     &                             * ( pe3_in(ji,jj,jk-1) - e3v_0(ji,jj,jk-1) )                              &
-                     &                             +            0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )    &
-                     &                             * ( pe3_in(ji,jj,jk  ) - e3v_0(ji,jj,jk  ) )
-               END DO
-            END DO
-         END DO
-!$OMP END PARALLEL
+            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
+               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              &
+               &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
+               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
+         END DO
       END SELECT
       !
@@ -1225,5 +905,5 @@
                      sshb(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)           !!gm I don't understand that !
                      sshn(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)
-                     ssha(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)                     
+                     ssha(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)
                   ENDIF
                 ENDDO
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 7698)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 7753)
@@ -72,5 +72,5 @@
       INTEGER, DIMENSION(:,:), INTENT(out) ::   k_top, k_bot   ! ocean first and last level indices
       !
-      INTEGER  ::   ji, jj, jk                  ! dummy loop index
+      INTEGER  ::   jk                  ! dummy loop index
       INTEGER  ::   ioptio, ibat, ios   ! local integer
       REAL(wp) ::   zrefdep             ! depth of the reference level (~10m)
@@ -114,20 +114,8 @@
 !!gm to be remove when removing the OLD definition of e3 scale factors so that gde3w disappears
       ! Compute gde3w_0 (vertical sum of e3w)
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)
-         END DO
+      gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
+      DO jk = 2, jpk
+         gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
       END DO
-      DO jk = 2, jpk
-!$OMP DO schedule(static) private(jj, ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-!$OMP END PARALLEL
       !
       IF(lwp) THEN                     ! Control print
@@ -202,5 +190,5 @@
       INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top , k_bot               ! first & last ocean level
       !
-      INTEGER  ::   jk, jj, ji   ! dummy loop index
+      INTEGER  ::   jk     ! dummy loop index
       INTEGER  ::   inum   ! local logical unit
       REAL(WP) ::   z_zco, z_zps, z_sco, z_cav
@@ -266,17 +254,7 @@
       !                          !* ocean top and bottom level
       CALL iom_get( inum, jpdom_data, 'top_level'    , z2d  , lrowattr=ln_use_jattr )   ! 1st wet T-points (ISF)
-!$OMP PARALLEL DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            k_top(ji,jj) = INT( z2d(ji,jj) )
-         END DO
-      END DO
+      k_top(:,:) = INT( z2d(:,:) )
       CALL iom_get( inum, jpdom_data, 'bottom_level' , z2d  , lrowattr=ln_use_jattr )   ! last wet T-points
-!$OMP PARALLEL DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            k_bot(ji,jj) = INT( z2d(ji,jj) )
-         END DO
-      END DO
+      k_bot(:,:) = INT( z2d(:,:) )
       !
       ! bathymetry with orography (wetting and drying only)
@@ -317,16 +295,10 @@
       IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
       !
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            mikt(ji,jj) = MAX( k_top(ji,jj) , 1 )    ! top    ocean k-index of T-level (=1 over land)
-            !
-            mbkt(ji,jj) = MAX( k_bot(ji,jj) , 1 )    ! bottom ocean k-index of T-level (=1 over land)
-         END DO
-      END DO
+      mikt(:,:) = MAX( k_top(:,:) , 1 )    ! top    ocean k-index of T-level (=1 over land)
+      !
+      mbkt(:,:) = MAX( k_bot(:,:) , 1 )    ! bottom ocean k-index of T-level (=1 over land)
+ 
       !                                    ! N.B.  top     k-index of W-level = mikt
       !                                    !       bottom  k-index of W-level = mbkt+1
-!$OMP DO schedule(static) private(jj, ji)
       DO jj = 1, jpjm1
          DO ji = 1, jpim1
@@ -340,79 +312,10 @@
       END DO
       ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 
-!$OMP DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zk(ji,jj) = REAL( miku(ji,jj), wp )
-         END DO
-      END DO
-!$OMP END PARALLEL
-      CALL lbc_lnk( zk, 'U', 1. )
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            miku(ji,jj) = MAX( INT( zk(ji,jj) ), 1 )
-         END DO
-      END DO
-!$OMP DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zk(ji,jj) = REAL( mikv(ji,jj), wp )
-         END DO
-      END DO
-!$OMP END PARALLEL
-      CALL lbc_lnk( zk, 'V', 1. )
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            mikv(ji,jj) = MAX( INT( zk(ji,jj) ), 1 )
-         END DO
-      END DO
-!$OMP DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zk(ji,jj) = REAL( mikf(ji,jj), wp )
-         END DO
-      END DO
-!$OMP END PARALLEL
-      CALL lbc_lnk( zk, 'F', 1. )
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            mikf(ji,jj) = MAX( INT( zk(ji,jj) ), 1 )
-         END DO
-      END DO
-      !
-!$OMP DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zk(ji,jj) = REAL( mbku(ji,jj), wp )
-         END DO
-      END DO
-!$OMP END PARALLEL
-      CALL lbc_lnk( zk, 'U', 1. )
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            mbku(ji,jj) = MAX( INT( zk(ji,jj) ), 1 )
-         END DO
-      END DO
-!$OMP DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zk(ji,jj) = REAL( mbkv(ji,jj), wp )
-         END DO
-      END DO
-!$OMP END PARALLEL
-      CALL lbc_lnk( zk, 'V', 1. )
-!$OMP PARALLEL DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            mbkv(ji,jj) = MAX( INT( zk(ji,jj) ), 1 )
-         END DO
-      END DO
+      zk(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk( zk, 'U', 1. )   ;   miku(:,:) = MAX( INT( zk(:,:) ), 1 )
+      zk(:,:) = REAL( mikv(:,:), wp )   ;   CALL lbc_lnk( zk, 'V', 1. )   ;   mikv(:,:) = MAX( INT( zk(:,:) ), 1 )
+      zk(:,:) = REAL( mikf(:,:), wp )   ;   CALL lbc_lnk( zk, 'F', 1. )   ;   mikf(:,:) = MAX( INT( zk(:,:) ), 1 )
+      !
+      zk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk( zk, 'U', 1. )   ;   mbku(:,:) = MAX( INT( zk(:,:) ), 1 )
+      zk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk( zk, 'V', 1. )   ;   mbkv(:,:) = MAX( INT( zk(:,:) ), 1 )
       !
       CALL wrk_dealloc( jpi,jpj,   zk )
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90	(revision 7698)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90	(revision 7753)
@@ -161,5 +161,4 @@
          ij0 = 101   ;   ij1 = 109                       ! Reduced T & S in the Alboran Sea
          ii0 = 141   ;   ii1 = 155
-!$OMP PARALLEL DO schedule(static) private(jj, ji)
          DO jj = mj0(ij0), mj1(ij1)
             DO ji = mi0(ii0), mi1(ii1)
@@ -182,13 +181,6 @@
 !!gm end
       !
-!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jk)    ! NO mask
-               ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jk)
-            END DO
-         END DO
-      END DO
+      ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:)    ! NO mask
+      ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:) 
       !
       IF( ln_sco ) THEN                   !==   s- or mixed s-zps-coordinate   ==!
@@ -201,5 +193,4 @@
          ENDIF
          !
-!$OMP PARALLEL DO schedule(static) private(jj, ji, jk, zl, jkk, zi)
          DO jj = 1, jpj                         ! vertical interpolation of T & S
             DO ji = 1, jpi
@@ -235,16 +226,8 @@
       ELSE                                !==   z- or zps- coordinate   ==!
          !                             
-!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ptsd(ji,jj,jk,jp_tem) = ptsd(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)    ! Mask
-                  ptsd(ji,jj,jk,jp_sal) = ptsd(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+         ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) * tmask(:,:,:)    ! Mask
+         ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) * tmask(:,:,:)
          !
          IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level
-!$OMP PARALLEL DO schedule(static) private(jj, ji, ik, zl)
             DO jj = 1, jpj
                DO ji = 1, jpi
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 7698)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 7753)
@@ -59,5 +59,5 @@
       !! ** Purpose :   Initialization of the dynamics and tracer fields.
       !!----------------------------------------------------------------------
-      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
       REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace
       !!----------------------------------------------------------------------
@@ -75,26 +75,9 @@
 !      IF( lk_c1d )   CALL dta_uvd_init        ! Initialization of U & V input data
 !!gm
-!$OMP PARALLEL
-      DO jn = 1, jpts
-!$OMP DO schedule(static) private(jk, jj, ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  tsa  (ji,jj,jk,jn) = 0._wp                                       ! set one for all to 0 at level jpk
-                  rab_b(ji,jj,jk,jn) = 0._wp   ;   rab_n(ji,jj,jk,jn) = 0._wp      ! set one for all to 0 at level jpk
-               END DO
-            END DO
-         END DO
-      END DO
-!$OMP DO schedule(static) private(jk, jj, ji)
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               rhd  (ji,jj,jk  ) = 0._wp   ;   rhop (ji,jj,jk  ) = 0._wp      ! set one for all to 0 at level jpk
-               rn2b (ji,jj,jk  ) = 0._wp   ;   rn2  (ji,jj,jk  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
-            END DO
-         END DO
-      END DO
-!$OMP END PARALLEL
+
+      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
+      rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
+      tsa  (:,:,:,:) = 0._wp                                   ! set one for all to 0 at level jpk
+      rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk
 
       IF( ln_rstart ) THEN                    ! Restart from a file
@@ -114,54 +97,16 @@
             CALL dta_tsd( nit000, tsb )       ! read 3D T and S data at nit000
             !
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj, ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  sshb (ji,jj)   = 0._wp      ! set the ocean at rest
-               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
-                     ub  (ji,jj,jk) = 0._wp
-                     vb  (ji,jj,jk) = 0._wp  
-                  END DO
-               END DO
-            END DO
-!$OMP END PARALLEL
+            sshb(:,:)   = 0._wp               ! set the ocean at rest
+            ub  (:,:,:) = 0._wp
+            vb  (:,:,:) = 0._wp  
             !
          ELSE                                 ! user defined initial T and S
             CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )         
          ENDIF
-!$OMP PARALLEL
-         DO jn = 1, jpts
-!$OMP DO schedule(static) private(jk, jj, ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     tsn  (ji,jj,jk,jn) = tsb (ji,jj,jk,jn)       ! set now values from to before ones
-                  END DO
-               END DO
-            END DO
-         END DO
-!$OMP DO schedule(static) private(jk, jj, ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  un   (ji,jj,jk)   = ub  (ji,jj,jk)
-                  vn   (ji,jj,jk)   = vb  (ji,jj,jk)
-               END DO
-            END DO
-         END DO
-!$OMP DO schedule(static) private(jj, ji)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               sshn (ji,jj)     = sshb(ji,jj)   
-               hdivn(ji,jj,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level
-            END DO
-         END DO
-!$OMP END PARALLEL
+         tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
+         sshn (:,:)     = sshb(:,:)   
+         un   (:,:,:)   = ub  (:,:,:)
+         vn   (:,:,:)   = vb  (:,:,:)
+         hdivn(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level
          CALL div_hor( 0 )                    ! compute interior hdivn value  
 !!gm                                    hdivn(:,:,:) = 0._wp
@@ -197,16 +142,9 @@
       ! Do it whatever the free surface method, these arrays being eventually used
       !
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            un_b(ji,jj) = 0._wp   ;   vn_b(ji,jj) = 0._wp
-            ub_b(ji,jj) = 0._wp   ;   vb_b(ji,jj) = 0._wp
-         END DO
-      END DO
+      un_b(:,:) = 0._wp   ;   vn_b(:,:) = 0._wp
+      ub_b(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
       !
 !!gm  the use of umsak & vmask is not necessary below as un, vn, ub, vb are always masked
       DO jk = 1, jpkm1
-!$OMP DO schedule(static) private(jj, ji)
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -220,15 +158,9 @@
       END DO
       !
-!$OMP DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            un_b(ji,jj) = un_b(ji,jj) * r1_hu_n(ji,jj)
-            vn_b(ji,jj) = vn_b(ji,jj) * r1_hv_n(ji,jj)
-            !
-            ub_b(ji,jj) = ub_b(ji,jj) * r1_hu_b(ji,jj)
-            vb_b(ji,jj) = vb_b(ji,jj) * r1_hv_b(ji,jj)
-         END DO
-      END DO
-!$OMP END PARALLEL
+      un_b(:,:) = un_b(:,:) * r1_hu_n(:,:)
+      vn_b(:,:) = vn_b(:,:) * r1_hv_n(:,:)
+      !
+      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
+      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
       !
       IF( nn_timing == 1 )   CALL timing_stop('istate_init')
