Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/LIM_SRC_2/limhdf_2.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/LIM_SRC_2/limhdf_2.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/LIM_SRC_2/limhdf_2.F90	(revision 7525)
@@ -125,5 +125,5 @@
          IF( iter == 1 ) THEN
 !$OMP PARALLEL DO schedule(static) private(jj, ji)
-             DO jj= 2, jpjm1                                       ! diffusive trend : divergence of the fluxes
+             DO jj= 2, jpjm1
                 DO ji = fs_2 , fs_jpim1   ! vector opt. 
                    zdiv0(ji,jj) = zdiv(ji,jj)              ! save the 1st evaluation of the diffusive trend in zdiv0
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 7525)
@@ -132,8 +132,9 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-      zqnsoce(ji,jj) = qns(ji,jj)
+            zqnsoce(ji,jj) = qns(ji,jj)
          END DO
       END DO
 !$OMP END DO NOWAIT
+
 !$OMP DO schedule(static) private(jj,ji,zinda,ifvt,i1mfr,idfr,iflt,ial,iadv,ifral,ifrdv,zqsr,zqns,zqhc,zemp,zemp_snw,zfmm,zfsalt,zcd) 
       DO jj = 1, jpj
@@ -437,5 +438,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  utau_oce(ji,jj) = utau(ji,jj)                 !* save the air-ocean stresses at ice time-step
+                  utau_oce(ji,jj) = utau(ji,jj)       !* save the air-ocean stresses at ice time-step
                   vtau_oce(ji,jj) = vtau(ji,jj)
                END DO
@@ -497,5 +498,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            soce_0(ji,jj) = soce                     ! constant SSS and ice salinity used in levitating sea-ice case
+            soce_0(ji,jj) = soce             ! constant SSS and ice salinity used in levitating sea-ice case
             sice_0(ji,jj) = sice
          END DO
@@ -520,21 +521,21 @@
       ELSE
 !$OMP PARALLEL DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-         snwice_mass  (ji,jj) = 0.e0           ! no mass exchanges
-         snwice_mass_b(ji,jj) = 0.e0           ! no mass exchanges
-         snwice_fmass (ji,jj) = 0.e0           ! no mass exchanges
-         END DO
-      END DO
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               snwice_mass  (ji,jj) = 0.e0           ! no mass exchanges
+               snwice_mass_b(ji,jj) = 0.e0           ! no mass exchanges
+               snwice_fmass (ji,jj) = 0.e0           ! no mass exchanges
+            END DO
+         END DO
       ENDIF
       IF( nn_ice_embd == 2 .AND.          &  ! full embedment (case 2) & no restart : 
          &   .NOT.ln_rstart ) THEN           ! deplete the initial ssh below sea-ice area
 !$OMP PARALLEL DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-         sshn(ji,jj) = sshn(ji,jj) - snwice_mass(ji,jj) * r1_rau0
-         sshb(ji,jj) = sshb(ji,jj) - snwice_mass(ji,jj) * r1_rau0
-         END DO
-      END DO
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               sshn(ji,jj) = sshn(ji,jj) - snwice_mass(ji,jj) * r1_rau0
+               sshb(ji,jj) = sshb(ji,jj) - snwice_mass(ji,jj) * r1_rau0
+            END DO
+         END DO
 !!gm I really don't like this staff here...  Find a way to put that elsewhere or differently
 !!gm
@@ -542,9 +543,13 @@
 
 !$OMP PARALLEL
-!$OMP DO schedule(static) private(jk)
-            do jk = 1,jpkm1                     ! adjust initial vertical scale factors
-               e3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
-               e3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
-            end do
+!$OMP DO schedule(static) private(jk, jj, ji)
+            DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk)*( 1._wp + sshn(ji,jj)*tmask(ji,jj,1)/(ht_0(ji,jj) + 1.0 - tmask(ji,jj,1)) )
+                     e3t_b(ji,jj,jk) = e3t_0(ji,jj,jk)*( 1._wp + sshb(ji,jj)*tmask(ji,jj,1)/(ht_0(ji,jj) + 1.0 - tmask(ji,jj,1)) )
+                  END DO
+               END DO
+            END DO
 !$OMP DO schedule(static) private(jk, jj, ji)
             DO jk = 1, jpk
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90	(revision 7525)
@@ -446,5 +446,4 @@
             fr_i  (ji,jj) = 1.0 - frld(ji,jj)  
             hicifp(ji,jj) = hicif(ji,jj) * fr_i(ji,jj) - hicifp(ji,jj)
-            ztmp(ji,jj) = 1. - pfrld(ji,jj)                              ! fraction of ice after the dynamic, before the thermodynamic
          END DO
       END DO
@@ -452,4 +451,10 @@
       ! Outputs
       !--------------------------------------------------------------------------------
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ztmp(ji,jj) = 1. - pfrld(ji,jj)                              ! fraction of ice after the dynamic, before the thermodynamic
+         END DO
+      END DO
       IF( iom_use('ist_cea'    ) )   CALL iom_put( 'ist_cea', (sist(:,:) - rt0) * ztmp(:,:) )   ! Ice surface temperature   [Celius]
       IF( iom_use('qsr_ai_cea' ) )   CALL iom_put( 'qsr_ai_cea', qsr_ice(:,:,1) * ztmp(:,:) )   ! Solar flux over the ice     [W/m2]
@@ -469,5 +474,10 @@
       IF( iom_use('bicemel_cea') )   CALL iom_put( 'bicemel_cea' , zdvomif(:,:) * zztmp     )   ! Melt at Sea Ice bottom   [kg/m2/s]
       IF( iom_use('licepro_cea') ) THEN
-         zlicegr(:,:) = MAX( 0.e0, rdm_ice(:,:)-zlicegr(:,:) )
+!$OMP PARALLEL DO schedule(static) private(jj,ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zlicegr(ji,jj) = MAX( 0.e0, rdm_ice(ji,jj)-zlicegr(ji,jj) )
+            END DO
+         END DO
                                      CALL iom_put( 'licepro_cea' , zlicegr(:,:) * zztmp     )   ! Lateral sea ice growth   [kg/m2/s]
       ENDIF
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90	(revision 7525)
@@ -55,6 +55,7 @@
       REAL(wp) :: zarea, zvol, zwei
       REAL(wp) :: ztemi(4), ztemo(4), zsali(4), zsalo(4), zflxi(4), zflxo(4)
-      REAL(wp) :: zt, zs, zu  
+      REAL(wp) :: zt, zs, zu
       REAL(wp) :: zsm0, zfwfnew
+      REAL(wp), DIMENSION(:,:)     ::   ztmp      ! 2D workspace
       IF( cp_cfg == "orca" .AND. jp_cfg == 1 .OR. jp_cfg == 2 .OR. jp_cfg == 4 ) THEN
       !!----------------------------------------------------------------------
@@ -63,5 +64,4 @@
       ! Mean global salinity
       zsm0 = 34.72654
-
       ! To compute fwf mean value mean fwf
 
@@ -72,5 +72,11 @@
          a_salb   = 0.e0 ! valeur de sal au debut de la simulation
          ! sshb used because diafwb called after tranxt (i.e. after the swap)
-         a_sshb = SUM( e1e2t(:,:) * sshb(:,:) * tmask_i(:,:) )
+!$OMP PARALLEL DO schedule(static) private(jj,ji) REDUCTION(+:a_sshb)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ztmp(ji,jj) = e1e2t(ji,jj) * sshb(ji,jj) * tmask_i(ji,jj)
+               a_sshb = a_sshb + ztmp(ji,jj)
+            END DO
+         END DO
          IF( lk_mpp )   CALL mpp_sum( a_sshb )      ! sum over the global domain
 
@@ -86,6 +92,13 @@
          IF( lk_mpp )   CALL mpp_sum( a_salb )      ! sum over the global domain
       ENDIF
-      
-      a_fwf    = SUM( e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * tmask_i(:,:) ) 
+        
+!$OMP PARALLEL DO schedule(static) private(jj,ji) REDUCTION(+:a_fwf)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ztmp(ji,jj) = e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj)
+            a_fwf = a_fwf + ztmp(ji,jj)
+         END DO
+      END DO
+
       IF( lk_mpp )   CALL mpp_sum( a_fwf    )       ! sum over the global domain
 
@@ -97,7 +110,19 @@
          zfwfnew = 0.e0
          ! Mean sea level at nitend
-         a_sshn = SUM( e1e2t(:,:) * sshn(:,:) * tmask_i(:,:) )
+!$OMP PARALLEL DO schedule(static) private(jj,ji) REDUCTION(+:a_sshn)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ztmp(ji,jj) = e1e2t(ji,jj) * sshn(ji,jj) * tmask_i(ji,jj)
+               a_sshn = a_sshn + ztmp(ji,jj)
+            END DO
+         END DO
          IF( lk_mpp )   CALL mpp_sum( a_sshn )      ! sum over the global domain
-         zarea  = SUM( e1e2t(:,:) *             tmask_i(:,:) )
+!$OMP PARALLEL DO schedule(static) private(jj,ji) REDUCTION(+:zarea)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ztmp(ji,jj) = e1e2t(ji,jj) *            tmask_i(ji,jj)
+               zarea = zarea + ztmp(ji,jj)
+            END DO
+         END DO
          IF( lk_mpp )   CALL mpp_sum( zarea  )      ! sum over the global domain
          
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 7525)
@@ -241,7 +241,11 @@
             END DO
          END DO
-!$OMP DO schedule(static) private(jk)
+!$OMP DO schedule(static) private(jk,jj,ji)
          DO jk = 1, jpk
-            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:) 
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  z3d(ji,jj,jk) = wn(ji,jj,jk) * z2d(ji,jj) 
+               END DO
+            END DO
          END DO
 !$OMP END DO NOWAIT
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90	(revision 7525)
@@ -63,39 +63,39 @@
       CASE( 'U' )
 !$OMP PARALLEL DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zglam(ji,jj) = glamu(ji,jj) ; zgphi(ji,jj) = gphiu(ji,jj)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zglam(ji,jj) = glamu(ji,jj) ; zgphi(ji,jj) = gphiu(ji,jj)
+            END DO
          END DO
-      END DO
-      zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,ik)
+         zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,ik)
       CASE( 'V' )  
 !$OMP PARALLEL DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zglam(ji,jj) = glamv(ji,jj) ; zgphi(ji,jj) = gphiv(ji,jj)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zglam(ji,jj) = glamv(ji,jj) ; zgphi(ji,jj) = gphiv(ji,jj)
+            END DO
          END DO
-      END DO
-      zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,ik)
+         zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,ik)
       CASE( 'F' )  
 !$OMP PARALLEL DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zglam(ji,jj) = glamf(ji,jj) ; zgphi(ji,jj) = gphif(ji,jj)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zglam(ji,jj) = glamf(ji,jj) ; zgphi(ji,jj) = gphif(ji,jj)
+            END DO
          END DO
-      END DO
-      zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,ik)
+         zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,ik)
       CASE DEFAULT 
 !$OMP PARALLEL DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zglam(ji,jj) = glamt(ji,jj) ; zgphi(ji,jj) = gphit(ji,jj)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zglam(ji,jj) = glamt(ji,jj) ; zgphi(ji,jj) = gphit(ji,jj)
+            END DO
          END DO
-      END DO
-      zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,ik)
+         zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,ik)
       END SELECT
 
       IF (jphgr_msh /= 2 .AND. jphgr_msh /= 3) THEN
          zlon       = MOD( plon       + 720., 360. )                                     ! plon between    0 and 360
-!$OMP PARALLEL DO schedule(static) private(jj, ji, zglam, zlon)
+!$OMP PARALLEL DO schedule(static) private(jj, ji, zlon)
          DO jj = 1, jpj
             DO ji = 1, jpi
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90	(revision 7525)
@@ -66,12 +66,12 @@
            CALL wrk_alloc( jpi,jpj,jpk,   ztrdu, ztrdv )
 !$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
-        DO jk = 1, jpk
-           DO jj = 1, jpj
-              DO ji = 1, jpi
-                 ztrdu(ji,jj,jk) = ua(ji,jj,jk)
-                 ztrdv(ji,jj,jk) = va(ji,jj,jk)
+           DO jk = 1, jpk
+              DO jj = 1, jpj
+                 DO ji = 1, jpi
+                    ztrdu(ji,jj,jk) = ua(ji,jj,jk)
+                    ztrdv(ji,jj,jk) = va(ji,jj,jk)
+                 END DO
               END DO
            END DO
-        END DO
         ENDIF
 
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 7525)
@@ -93,12 +93,12 @@
          CALL wrk_alloc( jpi,jpj,jpk,   ztrdu, ztrdv )
 !$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
-        DO jk = 1, jpk
-           DO jj = 1, jpj
-              DO ji = 1, jpi
-                 ztrdu(ji,jj,jk) = ua(ji,jj,jk)
-                 ztrdv(ji,jj,jk) = va(ji,jj,jk)
-              END DO
-           END DO
-        END DO
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdu(ji,jj,jk) = ua(ji,jj,jk)
+                  ztrdv(ji,jj,jk) = va(ji,jj,jk)
+               END DO
+            END DO
+         END DO
       ENDIF
 !$OMP PARALLEL DO schedule(static) private(jj, ji)
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 7525)
@@ -132,8 +132,12 @@
             END DO
          END DO
-!$OMP DO schedule(static) private(jk)
+!$OMP DO schedule(static) private(jk,jj,ji)
          DO jk = 1, jpkm1
-            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
-            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ua(ji,jj,jk) = ( ua(ji,jj,jk) - zue(ji,jj) * r1_hu_a(ji,jj) + ua_b(ji,jj) ) * umask(ji,jj,jk)
+                  va(ji,jj,jk) = ( va(ji,jj,jk) - zve(ji,jj) * r1_hv_a(ji,jj) + va_b(ji,jj) ) * vmask(ji,jj,jk)
+               END DO
+            END DO
          END DO
 !$OMP END DO NOWAIT
@@ -145,8 +149,12 @@
             ! In the forward case, this is done below after asselin filtering   
             ! so that asselin contribution is removed at the same time 
-!$OMP PARALLEL DO schedule(static) private(jk)
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
             DO jk = 1, jpkm1
-               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
-               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     un(ji,jj,jk) = ( un(ji,jj,jk) - un_adv(ji,jj) + un_b(ji,jj) )*umask(ji,jj,jk)
+                     vn(ji,jj,jk) = ( vn(ji,jj,jk) - vn_adv(ji,jj) + vn_b(ji,jj) )*vmask(ji,jj,jk)
+                  END DO
+               END DO
             END DO  
          ENDIF
@@ -198,8 +206,8 @@
                   zua(ji,jj,jk) = un(ji,jj,jk)             ! save the now velocity before the asselin filter
                   zva(ji,jj,jk) = vn(ji,jj,jk)             ! (caution: there will be a shift by 1 timestep in the
-               END DO
-            END DO
-         END DO
-         !                                  !  computation of the asselin filter trends)
+                                                           !  computation of the asselin filter trends)
+               END DO
+            END DO
+         END DO
       ENDIF
 
@@ -208,16 +216,24 @@
       IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
 !$OMP PARALLEL
-!$OMP DO schedule(static) private(jk)
+!$OMP DO schedule(static) private(jk,jj,ji)
          DO jk = 1, jpkm1
-            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
-            vn(:,:,jk) = va(:,:,jk)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  un(ji,jj,jk) = ua(ji,jj,jk)                          ! un <-- ua
+                  vn(ji,jj,jk) = va(ji,jj,jk)
+               END DO
+            END DO
          END DO
 !$OMP END DO NOWAIT
          IF(.NOT.ln_linssh ) THEN
-!$OMP DO schedule(static) private(jk)
+!$OMP DO schedule(static) private(jk,jj,ji)
             DO jk = 1, jpkm1
-               e3t_b(:,:,jk) = e3t_n(:,:,jk)
-               e3u_b(:,:,jk) = e3u_n(:,:,jk)
-               e3v_b(:,:,jk) = e3v_n(:,:,jk)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     e3t_b(ji,jj,jk) = e3t_n(ji,jj,jk)
+                     e3u_b(ji,jj,jk) = e3u_n(ji,jj,jk)
+                     e3v_b(ji,jj,jk) = e3v_n(ji,jj,jk)
+                  END DO
+               END DO
             END DO
 !$OMP END DO NOWAIT
@@ -256,7 +272,11 @@
                END DO
             ELSE
-!$OMP PARALLEL DO schedule(static) private(jk)
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
                DO jk = 1, jpkm1
-                  e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
+                  DO jj = 1, jpj
+                     DO ji = 1, jpi    
+                        e3t_b(ji,jj,jk) = e3t_n(ji,jj,jk) + atfp * ( e3t_b(ji,jj,jk) - 2._wp * e3t_n(ji,jj,jk) + e3t_a(ji,jj,jk) )
+                     END DO
+                  END DO
                END DO
                ! Add volume filter correction: compatibility with tracer advection scheme
@@ -365,8 +385,12 @@
                END DO
             END DO
-!$OMP DO schedule(static) private(jk)
+!$OMP DO schedule(static) private(jk,jj,ji)
             DO jk = 1, jpkm1
-               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
-               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ub(ji,jj,jk) = ub(ji,jj,jk) - (zue(ji,jj) * r1_hu_n(ji,jj) - un_b(ji,jj)) * umask(ji,jj,jk)
+                     vb(ji,jj,jk) = vb(ji,jj,jk) - (zve(ji,jj) * r1_hv_n(ji,jj) - vn_b(ji,jj)) * vmask(ji,jj,jk)
+                  END DO
+               END DO
             END DO
 !$OMP END DO NOWAIT
@@ -446,5 +470,5 @@
       ENDIF
       IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
-!$OMP DO schedule(static) private(jk, jj, ji)
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
          DO jk = 1, jpkm1
             DO jj = 1, jpj
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 7525)
@@ -1314,8 +1314,12 @@
       ! Update barotropic trend:
       IF( ln_dynadv_vec .OR. ln_linssh ) THEN
-!$OMP PARALLEL DO schedule(static) private(jk)
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
          DO jk=1,jpkm1
-            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
-            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ua(ji,jj,jk) = ua(ji,jj,jk) + ( ua_b(ji,jj) - ub_b(ji,jj) ) * z1_2dt_b
+                  va(ji,jj,jk) = va(ji,jj,jk) + ( va_b(ji,jj) - vb_b(ji,jj) ) * z1_2dt_b
+               END DO
+            END DO
          END DO
       ELSE
@@ -1335,8 +1339,12 @@
          !
 !$OMP PARALLEL
-!$OMP DO schedule(static) private(jk)
+!$OMP DO schedule(static) private(jk,jj,ji)
          DO jk=1,jpkm1
-            ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b
-            va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ua(ji,jj,jk) = ua(ji,jj,jk) + r1_hu_n(ji,jj) * ( ua_b(ji,jj) - ub_b(ji,jj) * hu_b(ji,jj) ) * z1_2dt_b
+                  va(ji,jj,jk) = va(ji,jj,jk) + r1_hv_n(ji,jj) * ( va_b(ji,jj) - vb_b(ji,jj) * hv_b(ji,jj) ) * z1_2dt_b
+               END DO
+            END DO
          END DO
 !$OMP END DO NOWAIT
@@ -1352,10 +1360,14 @@
       ENDIF
       !
-!$OMP PARALLEL DO schedule(static) private(jk)
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
       DO jk = 1, jpkm1
-         ! Correct velocities:
-         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk)
-         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk)
-         !
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ! Correct velocities:
+               un(ji,jj,jk) = ( un(ji,jj,jk) + un_adv(ji,jj) - un_b(ji,jj) ) * umask(ji,jj,jk)
+               vn(ji,jj,jk) = ( vn(ji,jj,jk) + vn_adv(ji,jj) - vn_b(ji,jj) ) * vmask(ji,jj,jk)
+               !
+            END DO
+         END DO
       END DO
       !
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 7525)
@@ -108,42 +108,42 @@
          IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztrdu(ji,jj,jk) = ua(ji,jj,jk)
-                  ztrdv(ji,jj,jk) = va(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdu(ji,jj,jk) = ua(ji,jj,jk)
+                     ztrdv(ji,jj,jk) = va(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
             CALL vor_ene( kt, nrvm, ua, va )                      ! relative vorticity or metric trend
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
-                  ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
+                     ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztrdu(ji,jj,jk) = ua(ji,jj,jk)
-                  ztrdv(ji,jj,jk) = va(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdu(ji,jj,jk) = ua(ji,jj,jk)
+                     ztrdv(ji,jj,jk) = va(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
             CALL vor_ene( kt, ncor, ua, va )                      ! planetary vorticity trend
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
-                  ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
+                     ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
@@ -154,42 +154,42 @@
          IF( l_trddyn ) THEN                                ! trend diagnostics: splitthe trend in two    
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztrdu(ji,jj,jk) = ua(ji,jj,jk)
-                  ztrdv(ji,jj,jk) = va(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdu(ji,jj,jk) = ua(ji,jj,jk)
+                     ztrdv(ji,jj,jk) = va(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
             CALL vor_ens( kt, nrvm, ua, va )                      ! relative vorticity or metric trend
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
-                  ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
+                     ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztrdu(ji,jj,jk) = ua(ji,jj,jk)
-                  ztrdv(ji,jj,jk) = va(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdu(ji,jj,jk) = ua(ji,jj,jk)
+                     ztrdv(ji,jj,jk) = va(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
             CALL vor_ens( kt, ncor, ua, va )                      ! planetary vorticity trend
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
-                  ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
+                     ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
@@ -200,42 +200,42 @@
          IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztrdu(ji,jj,jk) = ua(ji,jj,jk)
-                  ztrdv(ji,jj,jk) = va(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdu(ji,jj,jk) = ua(ji,jj,jk)
+                     ztrdv(ji,jj,jk) = va(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
             CALL vor_ens( kt, nrvm, ua, va )                      ! relative vorticity or metric trend (ens)
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
-                  ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
+                     ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztrdu(ji,jj,jk) = ua(ji,jj,jk)
-                  ztrdv(ji,jj,jk) = va(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdu(ji,jj,jk) = ua(ji,jj,jk)
+                     ztrdv(ji,jj,jk) = va(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
             CALL vor_ene( kt, ncor, ua, va )                      ! planetary vorticity trend (ene)
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
-                  ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
+                     ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
@@ -247,42 +247,42 @@
          IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztrdu(ji,jj,jk) = ua(ji,jj,jk)
-                  ztrdv(ji,jj,jk) = va(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdu(ji,jj,jk) = ua(ji,jj,jk)
+                     ztrdv(ji,jj,jk) = va(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
             CALL vor_een( kt, nrvm, ua, va )                      ! relative vorticity or metric trend
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
-                  ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
+                     ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztrdu(ji,jj,jk) = ua(ji,jj,jk)
-                  ztrdv(ji,jj,jk) = va(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdu(ji,jj,jk) = ua(ji,jj,jk)
+                     ztrdv(ji,jj,jk) = va(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
             CALL vor_een( kt, ncor, ua, va )                      ! planetary vorticity trend
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
-                  ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
+                     ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 7525)
@@ -78,12 +78,12 @@
          CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
 !$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
-        DO jk = 1, jpk
-           DO jj = 1, jpj
-              DO ji = 1, jpi
-                 ztrdu(ji,jj,jk) = ua(ji,jj,jk)
-                 ztrdv(ji,jj,jk) = va(ji,jj,jk)
-              END DO
-           END DO
-        END DO
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdu(ji,jj,jk) = ua(ji,jj,jk)
+                  ztrdv(ji,jj,jk) = va(ji,jj,jk)
+               END DO
+            END DO
+         END DO
       ENDIF
       
@@ -145,12 +145,12 @@
       IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
 !$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
-           DO jk = 1, jpk
-              DO jj = 1, jpj
-                 DO ji = 1, jpi
-                    ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
-                    ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
-                 END DO
-              END DO
-           END DO
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
+                  ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
+               END DO
+            END DO
+         END DO
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
          CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 7525)
@@ -68,12 +68,12 @@
          CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ztrdu(ji,jj,jk) = ua(ji,jj,jk)
-               ztrdv(ji,jj,jk) = va(ji,jj,jk)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztrdu(ji,jj,jk) = ua(ji,jj,jk)
+                  ztrdv(ji,jj,jk) = va(ji,jj,jk)
+               END DO
             END DO
          END DO
-      END DO
       ENDIF
 
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 7525)
@@ -105,9 +105,9 @@
       DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
 !$OMP DO schedule(static) private(jj, ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi   ! vector opt.
-                  zhdiv(ji,jj) = zhdiv(ji,jj) + e3t_n(ji,jj,jk) * hdivn(ji,jj,jk)
-               END DO
-            END DO           
+         DO jj = 1, jpj
+            DO ji = 1, jpi   ! vector opt.
+               zhdiv(ji,jj) = zhdiv(ji,jj) + e3t_n(ji,jj,jk) * hdivn(ji,jj,jk)
+            END DO
+         END DO           
       END DO
 !$OMP END PARALLEL
@@ -120,9 +120,9 @@
       IF(ln_wd) CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt)
 !$OMP PARALLEL DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ssha(ji,jj) = (  sshb(ji,jj) - z2dt * ( zcoef * ( emp_b(ji,jj) + emp(ji,jj) ) + zhdiv(ji,jj) )  ) * ssmask(ji,jj)
-         END DO
-      END DO           
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ssha(ji,jj) = (  sshb(ji,jj) - z2dt * ( zcoef * ( emp_b(ji,jj) + emp(ji,jj) ) + zhdiv(ji,jj) )  ) * ssmask(ji,jj)
+            END DO
+         END DO           
       IF ( .NOT.ln_dynspg_ts ) THEN
          ! These lines are not necessary with time splitting since
@@ -143,9 +143,9 @@
          CALL ssh_asm_inc( kt )
 !$OMP PARALLEL DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ssha(ji,jj) = ssha(ji,jj) + z2dt * ssh_iau(ji,jj)
-         END DO
-      END DO           
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ssha(ji,jj) = ssha(ji,jj) + z2dt * ssh_iau(ji,jj)
+            END DO
+         END DO           
       ENDIF
 #endif
@@ -193,5 +193,10 @@
          IF(lwp) WRITE(numout,*) '~~~~~ '
          !
-         wn(ji,jj,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               wn(ji,jj,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
+            END DO
+         END DO 
       ENDIF
       !                                           !------------------------------!
@@ -221,6 +226,6 @@
             DO jj = 1, jpj
                DO ji = 1, jpi   ! vector opt.
-            wn(ji,jj,jk) = wn(ji,jj,jk+1) - ( e3t_n(ji,jj,jk) * hdivn(ji,jj,jk) + zhdiv(ji,jj,jk)    &
-               &                         + z1_2dt * ( e3t_a(ji,jj,jk) - e3t_b(ji,jj,jk) )     ) * tmask(ji,jj,jk)
+                  wn(ji,jj,jk) = wn(ji,jj,jk+1) - ( e3t_n(ji,jj,jk) * hdivn(ji,jj,jk) + zhdiv(ji,jj,jk)    &
+                  &                         + z1_2dt * ( e3t_a(ji,jj,jk) - e3t_b(ji,jj,jk) )     ) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -233,7 +238,7 @@
             DO jj = 1, jpj
                DO ji = 1, jpi   ! vector opt.
-            ! computation of w
-            wn(ji,jj,jk) = wn(ji,jj,jk+1) - (  e3t_n(ji,jj,jk) * hdivn(ji,jj,jk)                 &
-               &                         + z1_2dt * ( e3t_a(ji,jj,jk) - e3t_b(ji,jj,jk) )  ) * tmask(ji,jj,jk)
+                  ! computation of w
+                  wn(ji,jj,jk) = wn(ji,jj,jk+1) - (  e3t_n(ji,jj,jk) * hdivn(ji,jj,jk)                 &
+                  &                         + z1_2dt * ( e3t_a(ji,jj,jk) - e3t_b(ji,jj,jk) )  ) * tmask(ji,jj,jk)
                 END DO
             END DO
@@ -291,10 +296,10 @@
          & ( ln_bt_fw    .AND. ln_dynspg_ts )      ) THEN 
 !$OMP PARALLEL DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            sshb(ji,jj) = sshn(ji,jj)                              ! before <-- now
-            sshn(ji,jj) = ssha(ji,jj)                              ! now    <-- after  (before already = now)
-         END DO
-      END DO           
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               sshb(ji,jj) = sshn(ji,jj)                              ! before <-- now
+               sshn(ji,jj) = ssha(ji,jj)                              ! now    <-- after  (before already = now)
+            END DO
+         END DO           
          !
       ELSE           !==  Leap-Frog time-stepping: Asselin filter + swap  ==!
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90	(revision 7525)
@@ -160,13 +160,13 @@
             IF(lwp) WRITE(numout,*) '          momentum mixing coef. = constant '
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ahmt(ji,jj,jk) = zah0 * tmask(ji,jj,jk)
-               ahmf(ji,jj,jk) = zah0 * fmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-            !
+               DO jk = 1, jpk
+                  DO jj = 1, jpj
+                     DO ji = 1, jpi
+                        ahmt(ji,jj,jk) = zah0 * tmask(ji,jj,jk)
+                        ahmf(ji,jj,jk) = zah0 * fmask(ji,jj,jk)
+                     END DO
+                  END DO
+               END DO
+               !
          CASE(  10  )      !==  fixed profile  ==!
             IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( depth )'
@@ -189,8 +189,12 @@
 !!              do we introduce a scaling by the max value of the array, and then multiply by zah0 ????
 !!              better:  check that the max is <=1  i.e. it is a shape from 0 to 1, not a coef that has physical dimension
-!$OMP PARALLEL DO schedule(static) private(jk)
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
             DO jk = 2, jpkm1
-               ahmt(:,:,jk) = ahmt(:,:,1) * tmask(:,:,jk)
-               ahmf(:,:,jk) = ahmf(:,:,1) * fmask(:,:,jk)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ahmt(ji,jj,jk) = ahmt(ji,jj,1) * tmask(ji,jj,jk)
+                     ahmf(ji,jj,jk) = ahmf(ji,jj,1) * fmask(ji,jj,jk)
+                  END DO
+               END DO
             END DO
             !
@@ -208,8 +212,12 @@
 !!gm Question : info for LAP or BLP case  to take into account the SQRT in the bilaplacian case ????
 !!              do we introduce a scaling by the max value of the array, and then multiply by zah0 ????
-!$OMP PARALLEL DO schedule(static) private(jk)
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
             DO jk = 1, jpkm1
-               ahmt(:,:,jk) = ahmt(:,:,jk) * tmask(:,:,jk)
-               ahmf(:,:,jk) = ahmf(:,:,jk) * fmask(:,:,jk)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ahmt(ji,jj,jk) = ahmt(ji,jj,jk) * tmask(ji,jj,jk)
+                     ahmf(ji,jj,jk) = ahmf(ji,jj,jk) * fmask(ji,jj,jk)
+                  END DO
+               END DO
             END DO
             !
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90	(revision 7525)
@@ -158,5 +158,6 @@
       !
 !$OMP PARALLEL
-!$OMP DO schedule(static) private(jj,ji,zlam,zphi,zxnpt,zynpt,znnpt,zxnpu,zynpu,znnpu,zxnpv,zynpv,znnpv,zxnpf,zynpf,znnpf,zlan,zphh,zxvvt,zyvvt,znvvt,zxffu,zyffu,znffu,zxffv,zyffv,znffv,zxuuf,zyuuf,znuuf)
+!$OMP DO schedule(static) private(jj,ji,zlam,zphi,zxnpt,zynpt,znnpt,zxnpu,zynpu,znnpu,zxnpv,zynpv,znnpv,zxnpf) &
+!$OMP& private(zynpf,znnpf,zlan,zphh,zxvvt,zyvvt,znvvt,zxffu,zyffu,znffu,zxffv,zyffv,znffv,zxuuf,zyuuf,znuuf)
       DO jj = 2, jpjm1
          DO ji = fs_2, jpi   ! vector opt.
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90	(revision 7525)
@@ -272,5 +272,5 @@
       ztau_sais = 0.015
       ztaun = ztau - ztau_sais * COS( (ztime - ztimemax) / (ztimemin - ztimemax) * rpi )
-      ! module of wind stress and wind speed at T-point
+
       zcoef = 1. / ( zrhoa * zcdrag ) 
 !$OMP PARALLEL
@@ -285,4 +285,5 @@
       END DO
 
+      ! module of wind stress and wind speed at T-point
 !$OMP DO schedule(static) private(jj, ji, ztx, zty, zmod)
       DO jj = 2, jpjm1
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 7525)
@@ -282,13 +282,14 @@
                zst(ji,jj) = pst(ji,jj) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
 
+      ! ... components ( U10m - U_oce ) at T-point (unmasked)
+              zwnd_i(ji,jj) = 0.e0  
+              zwnd_j(ji,jj) = 0.e0
+            END DO
+         END DO
+
       ! ----------------------------------------------------------------------------- !
       !      0   Wind components and module at T-point relative to the moving ocean   !
       ! ----------------------------------------------------------------------------- !
 
-      ! ... components ( U10m - U_oce ) at T-point (unmasked)
-              zwnd_i(ji,jj) = 0.e0  
-              zwnd_j(ji,jj) = 0.e0
-            END DO
-         END DO
 #if defined key_cyclone
       CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012)
@@ -325,23 +326,29 @@
       ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave
       zztmp = 1. - albo
-      IF( ln_dm2dc ) THEN
-         qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
-      ELSE
-         qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
+      IF( ln_dm2dc ) THEN    ;    qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
+      ELSE                   ;    qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
       ENDIF
 
-!$OMP PARALLEL DO schedule(static) private(jj, ji)
+!$OMP PARALLEL
+!$OMP DO schedule(static) private(jj, ji)
       DO jj = 1, jpj
          DO ji = 1, jpi
             zqlw(ji,jj) = (  sf(jp_qlw)%fnow(ji,jj,1) - Stef * zst(ji,jj)*zst(ji,jj)*zst(ji,jj)*zst(ji,jj)  ) * tmask(ji,jj,1)   ! Long  Wave
+         END DO
+      END DO
+!OMP END DO NOWAIT
             ! ----------------------------------------------------------------------------- !
             !     II    Turbulent FLUXES                                                    !
             ! ----------------------------------------------------------------------------- !
 
+!$OMP DO schedule(static) private(jj, ji)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
             ! ... specific humidity at SST and IST
             zqsatw(ji,jj) = zcoef_qsatw * EXP( -5107.4 / zst(ji,jj) )
-
-         END DO
-      END DO
+         END DO
+      END DO
+!$OMP END PARALLEL
+
       ! ... NCAR Bulk formulae, computation of Cd, Ch, Ce at T-point :
       CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm,   &
@@ -388,17 +395,27 @@
       !  Turbulent fluxes over ocean
       ! -----------------------------
+      IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !! q_air and t_air are (or "are almost") given at 10m (wind reference height)
+               zevap(ji,jj) = rn_efac*MAX( 0._wp,     rhoa*Ce(ji,jj)*( zqsatw(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) )*wndm(ji,jj) ) ! Evaporation
+               zqsb (ji,jj) =                     cpa*rhoa*Ch(ji,jj)*( zst   (ji,jj) - sf(jp_tair)%fnow(ji,jj,1) )*wndm(ji,jj)   ! Sensible Heat
+            END DO
+         END DO
+      ELSE
+!$OMP PARALLEL DO schedule(static) private(jj, ji)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !! q_air and t_air are not given at 10m (wind reference height)
+               ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!!
+               zevap(ji,jj) = rn_efac*MAX( 0._wp,     rhoa*Ce(ji,jj)*( zqsatw(ji,jj) - zq_zu(ji,jj) )*wndm(ji,jj) )   ! Evaporation
+               zqsb (ji,jj) =                     cpa*rhoa*Ch(ji,jj)*( zst   (ji,jj) - zt_zu(ji,jj) )*wndm(ji,jj)     ! Sensible Heat
+            END DO
+         END DO
+      ENDIF
 !$OMP PARALLEL DO schedule(static) private(jj, ji)
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
-            !! q_air and t_air are (or "are almost") given at 10m (wind reference height)
-               zevap(ji,jj) = rn_efac*MAX( 0._wp,     rhoa*Ce(ji,jj)*( zqsatw(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) )*wndm(ji,jj) ) ! Evaporation
-               zqsb (ji,jj) =                     cpa*rhoa*Ch(ji,jj)*( zst   (ji,jj) - sf(jp_tair)%fnow(ji,jj,1) )*wndm(ji,jj)   ! Sensible Heat
-            ELSE
-            !! q_air and t_air are not given at 10m (wind reference height)
-            ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!!
-               zevap(ji,jj) = rn_efac*MAX( 0._wp,     rhoa*Ce(ji,jj)*( zqsatw(ji,jj) - zq_zu(ji,jj) )*wndm(ji,jj) )   ! Evaporation
-               zqsb (ji,jj) =                     cpa*rhoa*Ch(ji,jj)*( zst   (ji,jj) - zt_zu(ji,jj) )*wndm(ji,jj)     ! Sensible Heat
-            ENDIF
             zqla (ji,jj) = Lv * zevap(ji,jj)                                                              ! Latent Heat
          END DO
@@ -422,14 +439,14 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-      emp (ji,jj) = (  zevap(ji,jj)                                          &   ! mass flux (evap. - precip.)
-         &         - sf(jp_prec)%fnow(ji,jj,1) * rn_pfac  ) * tmask(ji,jj,1)
-      !
-      qns(ji,jj) = zqlw(ji,jj) - zqsb(ji,jj) - zqla(ji,jj)                                &   ! Downward Non Solar 
-         &     - sf(jp_snow)%fnow(ji,jj,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip
-         &     - zevap(ji,jj) * pst(ji,jj) * rcp                                      &   ! remove evap heat content at SST
-         &     + ( sf(jp_prec)%fnow(ji,jj,1) - sf(jp_snow)%fnow(ji,jj,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair
-         &     * ( sf(jp_tair)%fnow(ji,jj,1) - rt0 ) * rcp                          &
-         &     + sf(jp_snow)%fnow(ji,jj,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
-         &     * ( MIN( sf(jp_tair)%fnow(ji,jj,1), rt0_snow ) - rt0 ) * cpic * tmask(ji,jj,1)
+            emp (ji,jj) = (  zevap(ji,jj)                                          &   ! mass flux (evap. - precip.)
+               &         - sf(jp_prec)%fnow(ji,jj,1) * rn_pfac  ) * tmask(ji,jj,1)
+            !
+            qns(ji,jj) = zqlw(ji,jj) - zqsb(ji,jj) - zqla(ji,jj)                                &   ! Downward Non Solar 
+               &     - sf(jp_snow)%fnow(ji,jj,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip
+               &     - zevap(ji,jj) * pst(ji,jj) * rcp                                      &   ! remove evap heat content at SST
+               &     + ( sf(jp_prec)%fnow(ji,jj,1) - sf(jp_snow)%fnow(ji,jj,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair
+               &     * ( sf(jp_tair)%fnow(ji,jj,1) - rt0 ) * rcp                          &
+               &     + sf(jp_snow)%fnow(ji,jj,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
+               &     * ( MIN( sf(jp_tair)%fnow(ji,jj,1), rt0_snow ) - rt0 ) * cpic * tmask(ji,jj,1)
          END DO
       END DO
@@ -454,10 +471,10 @@
          CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean
 !$OMP PARALLEL DO schedule(static) private(jj, ji)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            tprecip(ji,jj) = sf(jp_prec)%fnow(ji,jj,1) * rn_pfac   ! output total precipitation [kg/m2/s]
-            sprecip(ji,jj) = sf(jp_snow)%fnow(ji,jj,1) * rn_pfac   ! output solid precipitation [kg/m2/s]
-         END DO
-      END DO
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               tprecip(ji,jj) = sf(jp_prec)%fnow(ji,jj,1) * rn_pfac   ! output total precipitation [kg/m2/s]
+               sprecip(ji,jj) = sf(jp_snow)%fnow(ji,jj,1) * rn_pfac   ! output solid precipitation [kg/m2/s]
+           END DO
+         END DO
          CALL iom_put( 'snowpre', sprecip * 86400. )        ! Snow
          CALL iom_put( 'precip' , tprecip * 86400. )        ! Total precipitation
@@ -598,5 +615,5 @@
       REAL(wp) ::   zst2, zst3
       REAL(wp) ::   zcoef_dqlw, zcoef_dqla, zcoef_dqsb
-      REAL(wp) ::   zztmp, z1_lsub                               ! temporary variable
+      REAL(wp) ::   zztmp, z1_lsub, ztmp1, ztmp2                 ! temporary variable
       !!
       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw             ! long wave heat flux over ice
@@ -706,5 +723,4 @@
 !$OMP END PARALLEL
       CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing 
-   
       emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
       emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
@@ -712,10 +728,10 @@
 
       ! --- heat flux associated with emp --- !
-      qemp_oce(:,:) = - pfrld(:,:) * zevap(:,:) * sst_m(:,:) * rcp                            & ! evap at sst
-      &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair
-      &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                    & ! solid precip at min(Tair,Tsnow)
-      &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
-      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                   & ! solid precip (only)
-      &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
+      qemp_oce(:,:) = - pfrld(:,:) * zevap(:,:) * sst_m(:,:) * rcp                               & ! evap at sst
+         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair
+         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow)
+         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
+      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
+         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
 
       ! --- total solar and non solar fluxes --- !
@@ -723,6 +739,5 @@
       qsr_tot(:,:) = pfrld(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
 
-      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo)
-      ! --- !
+      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
       qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
 
@@ -741,9 +756,11 @@
       ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 )
       !
+      ztmp1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
+      ztmp2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
 !$OMP PARALLEL DO schedule(static) private(jj, ji)
       DO jj = 1, jpj
          DO ji = 1, jpi
-            fr1_i0(ji,jj) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
-            fr2_i0(ji,jj) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
+            fr1_i0(ji,jj) = ztmp1
+            fr2_i0(ji,jj) = ztmp2
          END DO
       END DO
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7525)
@@ -195,6 +195,6 @@
             END DO
 !$OMP END DO NOWAIT
-!$OMP DO schedule(static) private(jp,jj,ji)
             DO jp = 1, jpts
+!$OMP DO schedule(static) private(jj,ji)
                DO jj = 1, jpj
                   DO ji = 1, jpi
@@ -336,13 +336,13 @@
          !                                         ! ---------------------------------------- !
 !$OMP PARALLEL DO schedule(static) private(jj,ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  utau_b(ji,jj) = utau(ji,jj)                         ! Swap the ocean forcing fields
-                  vtau_b(ji,jj) = vtau(ji,jj)                         ! (except at nit000 where before fields
-                  qns_b (ji,jj) = qns (ji,jj)                         !  are set at the end of the routine)
-                  emp_b (ji,jj) = emp (ji,jj)
-                  sfx_b (ji,jj) = sfx (ji,jj)
-               END DO
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               utau_b(ji,jj) = utau(ji,jj)                         ! Swap the ocean forcing fields
+               vtau_b(ji,jj) = vtau(ji,jj)                         ! (except at nit000 where before fields
+               qns_b (ji,jj) = qns (ji,jj)                         !  are set at the end of the routine)
+               emp_b (ji,jj) = emp (ji,jj)
+               sfx_b (ji,jj) = sfx (ji,jj)
             END DO
+         END DO
          IF ( ln_rnf ) THEN
 !$OMP PARALLEL
@@ -354,6 +354,6 @@
             END DO
 !$OMP END DO NOWAIT
-!$OMP DO schedule(static) private(jp,jj,ji)
             DO jp = 1, jpts
+!$OMP DO schedule(static) private(jj,ji)
                DO jj = 1, jpj
                   DO ji = 1, jpi
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 7525)
@@ -100,14 +100,18 @@
       !
       !                                         !==  effective transport  ==!
-!$OMP PARALLEL DO schedule(static) private(jk)
+!$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
       !
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN                                ! add z-tilde and/or vvl corrections
 !$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
-         DO jk = 1, jpk
+         DO jk = 1, jpkm1
             DO jj = 1, jpj
                DO ji = 1, jpi
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90	(revision 7525)
@@ -339,9 +339,9 @@
          IF( ln_linssh ) THEN    ! top ocean value: high order = upstream  ==>>  zwz=0
 !$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
+            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
          !
@@ -368,5 +368,5 @@
          !
          IF( l_trd ) THEN     ! trend diagnostics (contribution of upstream fluxes)
-!$OMP DO schedule(static) private(jk, jj, ji)
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
             DO jk = 1, jpk
                DO jj = 1, jpj
@@ -375,7 +375,7 @@
                      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
+                  END DO
+               END DO
+            END DO
             !
             CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90	(revision 7525)
@@ -114,15 +114,19 @@
 !$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
+            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)
+!$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
+               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
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 7525)
@@ -328,5 +328,5 @@
                 ztfw(1,jj,jk) = 0._wp     ;     ztfw(jpi,jj,jk) = 0._wp
              END DO
-          END DO
+         END DO
          !
          ! Vertical fluxes
@@ -338,5 +338,5 @@
                 ztfw(ji,jj, 1 ) = 0._wp      ;      ztfw(ji,jj,jpk) = 0._wp
              END DO
-          END DO
+         END DO
          
 !$OMP DO schedule(static) private(jk, jj, ji, zmsku, zmskv, zahu_w, zahv_w, zcoef3, zcoef4)
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 7525)
@@ -128,7 +128,12 @@
       IF( l_trdtra ) THEN      ! trends diagnostic: save the input temperature trend
          CALL wrk_alloc( jpi,jpj,jpk,   ztrdt ) 
-!$OMP PARALLEL WORKSHARE
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
-!$OMP END PARALLEL WORKSHARE
+!$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
       !
@@ -144,13 +149,23 @@
          ELSE                                           ! No restart or restart not found: Euler forward time stepping
             z1_2 = 1._wp
-!$OMP PARALLEL WORKSHARE
-            qsr_hc_b(:,:,:) = 0._wp
-!$OMP END PARALLEL WORKSHARE
+!$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
-!$OMP PARALLEL WORKSHARE
-         qsr_hc_b(:,:,:) = qsr_hc(:,:,:)
-!$OMP END PARALLEL WORKSHARE
+!$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
       !
@@ -161,7 +176,11 @@
       CASE( np_BIO )                   !==  bio-model fluxes  ==!
          !
-!$OMP PARALLEL DO schedule(static) private(jk)
+!$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
          !
@@ -198,7 +217,11 @@
             END DO
          ELSE                                !* constant chrlorophyll
-!$OMP PARALLEL DO schedule(static) private(jk)
+!$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
@@ -305,7 +328,10 @@
          !
 !$OMP PARALLEL
-!$OMP WORKSHARE
-         zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero
-!$OMP END WORKSHARE
+!$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
 !$OMP DO schedule(static) private(jj,ji)
@@ -329,7 +355,12 @@
       !
       IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
-!$OMP PARALLEL WORKSHARE
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
-!$OMP END PARALLEL WORKSHARE
+!$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 ) 
@@ -458,7 +489,12 @@
       END SELECT
       !
-!$OMP PARALLEL WORKSHARE
-      qsr_hc(:,:,:) = 0._wp     ! now qsr heat content set to zero where it will not be computed
-!$OMP END PARALLEL WORKSHARE
+!$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)
@@ -466,7 +502,10 @@
          CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev'  , fraqsr_1lev  )
       ELSE
-!$OMP PARALLEL WORKSHARE
-         fraqsr_1lev(:,:) = 1._wp   ! default : no penetration
-!$OMP END PARALLEL WORKSHARE
+!$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: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 7525)
@@ -900,10 +900,14 @@
       !                               !* set vertical eddy coef. to the background value
 !$OMP PARALLEL
-!$OMP DO schedule(static) private(jk)
+!$OMP DO schedule(static) private(jk,jj,ji)
       DO jk = 1, jpk
-         avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
-         avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
-         avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
-         avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               avt (ji,jj,jk) = avtb(jk) * wmask (ji,jj,jk)
+               avm (ji,jj,jk) = avmb(jk) * wmask (ji,jj,jk)
+               avmu(ji,jj,jk) = avmb(jk) * wumask(ji,jj,jk)
+               avmv(ji,jj,jk) = avmb(jk) * wvmask(ji,jj,jk)
+            END DO
+         END DO
       END DO
 !$OMP END DO NOWAIT
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90	(revision 7525)
@@ -529,11 +529,11 @@
 !$OMP PARALLEL
 !$OMP DO schedule(static) private(jk, jj, ji) 
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zav_tide(ji,jj,jk) = 0.e0
-            END DO
-         END DO
-      END DO
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zav_tide(ji,jj,jk) = 0.e0
+               END DO
+            END DO
+         END DO
 !$OMP DO schedule(static) private(jk)
          DO jk = 2, jpkm1
@@ -1024,5 +1024,6 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            znu_t(ji,jj,jk) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(ji,jj,jk,jp_tem) + 0.00694_wp * tsn(ji,jj,jk,jp_tem) * tsn(ji,jj,jk,jp_tem)  &
+            znu_t(ji,jj,jk) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(ji,jj,jk,jp_tem) &
+         &                                  + 0.00694_wp * tsn(ji,jj,jk,jp_tem) * tsn(ji,jj,jk,jp_tem)  &
          &                                  + 0.02305_wp * tsn(ji,jj,jk,jp_sal)  ) * tmask(ji,jj,jk) * r1_rau0
          END DO
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7525)
@@ -73,5 +73,5 @@
       !!              -8- Outputs and diagnostics
       !!----------------------------------------------------------------------
-      INTEGER ::   jk, jj, ji, jt       ! dummy loop indice
+      INTEGER ::   jn, jk, jj, ji       ! dummy loop indice
       INTEGER ::   indic    ! error indicator if < 0
       INTEGER ::   kcall    ! optional integer argument (dom_vvl_sf_nxt)
@@ -202,26 +202,13 @@
                &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level
 !!jc: fs simplification
-!$OMP PARALLEL
-!$OMP DO schedule(static) private(jk, jj, ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ua(ji,jj,jk) = 0._wp            ! set dynamics trends to zero
-                  va(ji,jj,jk) = 0._wp
-               END DO
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ua(ji,jj,jk) = 0._wp            ! set dynamics trends to zero
+               va(ji,jj,jk) = 0._wp
             END DO
          END DO
-!$OMP END DO NOWAIT
-!$OMP DO schedule(static) private(jt, jk, jj, ji)
-         DO jt = 1, jpts
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     tsa(ji,jj,jk,jt) = 0._wp         ! set tracer trends to zero
-                  END DO
-               END DO
-            END DO
-         END DO
-!$OMP END PARALLEL
+      END DO
       IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   &
                &         CALL dyn_asm_inc   ( kstp )  ! apply dynamics assimilation increment
@@ -276,4 +263,14 @@
       ! Active tracers                              
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      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
+                  tsa(ji,jj,jk,jn) = 0._wp         ! set tracer trends to zero
+               END DO
+            END DO
+         END DO
+      END DO
       IF(  lk_asminc .AND. ln_asmiau .AND. &
          & ln_trainc )   CALL tra_asm_inc   ( kstp )  ! apply tracer assimilation increment
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90	(revision 7525)
@@ -199,5 +199,6 @@
       ! CHEMICAL CONSTANTS - DEEP OCEAN
       ! -------------------------------
-!$OMP DO schedule(static) private(jk,jj,ji,zpres,ztkel,zsal,zsqrt,zsal15,zlogt,ztr,zis,zis2,zisqrt,ztc,zcl,zst,zft,zcks,zckf,zckb,zck1,zck2,zckw,zaksp0,zak1,zak2,zakb,zakw,zaksp1,zcpexp,zcpexp2,zbuf1,zbuf2)
+!$OMP DO schedule(static) private(jk,jj,ji,zpres,ztkel,zsal,zsqrt,zsal15,zlogt,ztr,zis,zis2,zisqrt,ztc,zcl,zst) &
+!$OMP& private(zft,zcks,zckf,zckb,zck1,zck2,zckw,zaksp0,zak1,zak2,zakb,zakw,zaksp1,zcpexp,zcpexp2,zbuf1,zbuf2)
       DO jk = 1, jpk
          DO jj = 1, jpj
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 7525)
@@ -136,5 +136,6 @@
             END DO
          END DO
-!$OMP DO schedule(static) private(jk,jj,ji,ztligand,zionic,zph,zoxy,zkox,zkph2,zkph1,ztfe,za,zb,zc,zkappa1,zkappa2,za2,za1,za0,zp,zq,zp3,zq2,zd,zr,zphi,zxs,zfff,zfunc)
+!$OMP DO schedule(static) private(jk,jj,ji,ztligand,zionic,zph,zoxy,zkox,zkph2,zkph1,ztfe,za) &
+!$OMP& private(zb,zc,zkappa1,zkappa2,za2,za1,za0,zp,zq,zp3,zq2,zd,zr,zphi,zxs,zfff,zfunc)
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -240,5 +241,6 @@
       zdust = 0.         ! if no dust available
       !
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,zkeq,zfesatur,ztfe,zstep,zfeequi,zfecoll,ztrc,zlam1b,zscave,zdenom1,zdenom2,zlamfac,zdep,zcoag,zlam1a,zaggdfea,zaggdfeb)
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,zkeq,zfesatur,ztfe,zstep,zfeequi,zfecoll) &
+!$OMP& private(ztrc,zlam1b,zscave,zdenom1,zdenom2,zlamfac,zdep,zcoag,zlam1a,zaggdfea,zaggdfeb)
       DO jk = 1, jpkm1
          DO jj = 1, jpj
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 7525)
@@ -84,5 +84,6 @@
       !
 !$OMP PARALLEL
-!$OMP DO schedule(static) private(jk,jj,ji,zno3,zferlim,zconcd,zconcd2,zconcn,zconcn2,z1_trbphy,z1_trbdia,zconc1d,zconc1dnh4,zconc0n,zconc0nnh4,zdenom,zlim1,zlim2,zlim3,zlim4,zratio,zironmin)
+!$OMP DO schedule(static) private(jk,jj,ji,zno3,zferlim,zconcd,zconcd2,zconcn,zconcn2,z1_trbphy,z1_trbdia) &
+!$OMP& private(zconc1d,zconc1dnh4,zconc0n,zconc0nnh4,zdenom,zlim1,zlim2,zlim3,zlim4,zratio,zironmin)
       DO jk = 1, jpkm1
          DO jj = 1, jpj
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90	(revision 7525)
@@ -99,5 +99,8 @@
       ENDIF
 
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,zcompam,zstep,zfact,zrespz2,ztortz2,zcompadi,zcompaz,zcompaph,zcompapoc,zfood,zfoodlim,zdenom,zdenom2,zgraze2,zgrazd,zgrazz,zgrazn,zgrazpoc,zgraznf,zgrazf,zgrazpof,zgrazffeg,zgrazfffg,zgrazffep,zgrazfffp,zgraztot,zproport,zratio,zratio2,zfrac,zfracfe,zgraztotf,zgrasrat,zgrasratn,zepshert,zepsherv,zgrarem2,zgrafer2,zgrapoc2,zgrarsig,zmortz2,zmortzgoc,zprcaca)
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,zcompam,zstep,zfact,zrespz2,ztortz2,zcompadi,zcompaz,zcompaph) &
+!$OMP& private(zcompapoc,zfood,zfoodlim,zdenom,zdenom2,zgraze2,zgrazd,zgrazz,zgrazn,zgrazpoc,zgraznf,zgrazf,zgrazpof) &
+!$OMP& private(zgrazffeg,zgrazfffg,zgrazffep,zgrazfffp,zgraztot,zproport,zratio,zratio2,zfrac,zfracfe,zgraztotf,zgrasrat) &
+!$OMP& private(zgrasratn,zepshert,zepsherv,zgrarem2,zgrafer2,zgrapoc2,zgrarsig,zmortz2,zmortzgoc,zprcaca)
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -269,11 +272,11 @@
          IF( iom_use( "PCAL" ) ) THEN
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zw3d(ji,jj,jk) = prodcal(ji,jj,jk) * 1.e+3 * rfact2r * tmask(ji,jj,jk)   !  Calcite production
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zw3d(ji,jj,jk) = prodcal(ji,jj,jk) * 1.e+3 * rfact2r * tmask(ji,jj,jk)   !  Calcite production
+                  END DO
+               END DO
             END DO
-         END DO
-      END DO
             CALL iom_put( "PCAL", zw3d )  
          ENDIF
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90	(revision 7525)
@@ -85,5 +85,7 @@
       IF( lk_iomput )  CALL wrk_alloc( jpi, jpj, jpk, zgrazing )
       !
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,zcompaz,zstep,zfact,zrespz,ztortz,zcompadi,zcompaph,zcompapoc,zfood,zfoodlim,zdenom,zdenom2,zgraze,zgrazp,zgrazm,zgrazsd,zgrazpf,zgrazmf,zgrazsf,zgraztot,zgraztotf,zgraztotn,zgrasrat,zgrasratn,zepshert,zepsherv,zgrafer,zgrarem,zgrapoc,zgrarsig,zmortz,zprcaca)
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,zcompaz,zstep,zfact,zrespz,ztortz,zcompadi,zcompaph,zcompapoc) &
+!$OMP& private(zfood,zfoodlim,zdenom,zdenom2,zgraze,zgrazp,zgrazm,zgrazsd,zgrazpf,zgrazmf,zgrazsf,zgraztot,zgraztotf) &
+!$OMP& private(zgraztotn,zgrasrat,zgrasratn,zepshert,zepsherv,zgrafer,zgrarem,zgrapoc,zgrarsig,zmortz,zprcaca)
       DO jk = 1, jpkm1
          DO jj = 1, jpj
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 7525)
@@ -90,13 +90,13 @@
 !$OMP PARALLEL
 !$OMP DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     ze1(ji,jj,jk) = 0._wp
-                     ze2(ji,jj,jk) = 0._wp
-                     ze3(ji,jj,jk) = 0._wp
-                  END DO
-               END DO
-            END DO
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ze1(ji,jj,jk) = 0._wp
+               ze2(ji,jj,jk) = 0._wp
+               ze3(ji,jj,jk) = 0._wp
+            END DO
+         END DO
+      END DO
 !$OMP END DO NOWAIT
       !                                        !* attenuation coef. function of Chlorophyll and wavelength (Red-Green-Blue)
@@ -249,9 +249,9 @@
       !
 !$OMP DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-      emoy(ji,jj,jk) = etot(ji,jj,jk)       ! remineralisation
-      zpar(ji,jj,jk) = etot_ndcy(ji,jj,jk)  ! diagnostic : PAR with no diurnal cycle 
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               emoy(ji,jj,jk) = etot(ji,jj,jk)       ! remineralisation
+               zpar(ji,jj,jk) = etot_ndcy(ji,jj,jk)  ! diagnostic : PAR with no diurnal cycle 
             END DO
          END DO
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 7525)
@@ -129,12 +129,12 @@
 
       IF( lk_degrad ) THEN
-!$OMP DO schedule(static) private(jk,jj,ji)
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               prmax(ji,jj,jk) = prmax(ji,jj,jk) * facvol(ji,jj,jk) 
-            END DO
-         END DO
-      END DO 
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  prmax(ji,jj,jk) = prmax(ji,jj,jk) * facvol(ji,jj,jk) 
+               END DO
+            END DO
+         END DO 
       END IF
 
@@ -453,161 +453,161 @@
           IF( iom_use( "PPPHY" ) .OR. iom_use( "PPPHY2" ) )  THEN
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = zprorca (ji,jj,jk) * zfact * tmask(ji,jj,jk)  ! primary production by nanophyto
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "PPPHY"  , zw3d )
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = zprorcad (ji,jj,jk) * zfact * tmask(ji,jj,jk)  ! primary production by nanophyto
-                  END DO
-               END DO
-            END DO
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = zprorca (ji,jj,jk) * zfact * tmask(ji,jj,jk)  ! primary production by nanophyto
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "PPPHY"  , zw3d )
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = zprorcad (ji,jj,jk) * zfact * tmask(ji,jj,jk)  ! primary production by nanophyto
+                   END DO
+                END DO
+             END DO
+             !
+             CALL iom_put( "PPPHY2"  , zw3d )
+          ENDIF
+          IF( iom_use( "PPNEWN" ) .OR. iom_use( "PPNEWD" ) )  THEN
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = zpronew (ji,jj,jk) * zfact * tmask(ji,jj,jk)  ! new primary production by nanophyto
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "PPNEWN"  , zw3d )
+             !
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = zpronewd (ji,jj,jk) * zfact * tmask(ji,jj,jk)  ! new primary production by nanophyto
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "PPNEWD"  , zw3d )
+          ENDIF
+          IF( iom_use( "PBSi" ) )  THEN
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = zprorcad(ji,jj,jk) * zfact * tmask(ji,jj,jk) * zysopt(ji,jj,jk) ! biogenic silica production
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "PBSi"  , zw3d )
+          ENDIF
+          IF( iom_use( "PFeN" ) .OR. iom_use( "PFeD" ) )  THEN
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = zprofen(ji,jj,jk) * zfact * tmask(ji,jj,jk)  ! biogenic iron production by nanophyto
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "PFeN"  , zw3d )
+             !
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = zprofed(ji,jj,jk) * zfact * tmask(ji,jj,jk)  ! biogenic iron production by nanophyto
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "PFeD"  , zw3d )
+          ENDIF
+          IF( iom_use( "Mumax" ) )  THEN
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = prmax(ji,jj,jk) * tmask(ji,jj,jk)   ! Maximum growth rate
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "Mumax"  , zw3d )
+          ENDIF
+          IF( iom_use( "MuN" ) .OR. iom_use( "MuD" ) )  THEN
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = zprbio(ji,jj,jk) * xlimphy(ji,jj,jk) * tmask(ji,jj,jk)  ! Realized growth rate for nanophyto
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "MuN"  , zw3d )
               !
-              CALL iom_put( "PPPHY2"  , zw3d )
-          ENDIF
-          IF( iom_use( "PPNEWN" ) .OR. iom_use( "PPNEWD" ) )  THEN
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = zpronew (ji,jj,jk) * zfact * tmask(ji,jj,jk)  ! new primary production by nanophyto
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "PPNEWN"  , zw3d )
-              !
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = zpronewd (ji,jj,jk) * zfact * tmask(ji,jj,jk)  ! new primary production by nanophyto
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "PPNEWD"  , zw3d )
-          ENDIF
-          IF( iom_use( "PBSi" ) )  THEN
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = zprorcad(ji,jj,jk) * zfact * tmask(ji,jj,jk) * zysopt(ji,jj,jk) ! biogenic silica production
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "PBSi"  , zw3d )
-          ENDIF
-          IF( iom_use( "PFeN" ) .OR. iom_use( "PFeD" ) )  THEN
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = zprofen(ji,jj,jk) * zfact * tmask(ji,jj,jk)  ! biogenic iron production by nanophyto
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "PFeN"  , zw3d )
-              !
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = zprofed(ji,jj,jk) * zfact * tmask(ji,jj,jk)  ! biogenic iron production by nanophyto
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "PFeD"  , zw3d )
-          ENDIF
-          IF( iom_use( "Mumax" ) )  THEN
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = prmax(ji,jj,jk) * tmask(ji,jj,jk)   ! Maximum growth rate
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "Mumax"  , zw3d )
-          ENDIF
-          IF( iom_use( "MuN" ) .OR. iom_use( "MuD" ) )  THEN
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = zprbio(ji,jj,jk) * xlimphy(ji,jj,jk) * tmask(ji,jj,jk)  ! Realized growth rate for nanophyto
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "MuN"  , zw3d )
-              !
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) =  zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * tmask(ji,jj,jk)  ! Realized growth rate for diatoms
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "MuD"  , zw3d )
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) =  zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * tmask(ji,jj,jk)  ! Realized growth rate for diatoms
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "MuD"  , zw3d )
           ENDIF
           IF( iom_use( "LNlight" ) .OR. iom_use( "LDlight" ) )  THEN
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = zprbio (ji,jj,jk) / (prmax(ji,jj,jk) + rtrn) * tmask(ji,jj,jk) ! light limitation term
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "LNlight"  , zw3d )
-              !
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) =  zprdia (ji,jj,jk) / (prmax(ji,jj,jk) + rtrn) * tmask(ji,jj,jk)  ! light limitation term
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "LDlight"  , zw3d )
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = zprbio (ji,jj,jk) / (prmax(ji,jj,jk) + rtrn) * tmask(ji,jj,jk) ! light limitation term
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "LNlight"  , zw3d )
+             !
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) =  zprdia (ji,jj,jk) / (prmax(ji,jj,jk) + rtrn) * tmask(ji,jj,jk)  ! light limitation term
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "LDlight"  , zw3d )
           ENDIF
           IF( iom_use( "TPP" ) )  THEN
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = ( zprorca(ji,jj,jk) + zprorcad(ji,jj,jk) ) * zfact * tmask(ji,jj,jk)  ! total primary production
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "TPP"  , zw3d )
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = ( zprorca(ji,jj,jk) + zprorcad(ji,jj,jk) ) * zfact * tmask(ji,jj,jk)  ! total primary production
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "TPP"  , zw3d )
           ENDIF
           IF( iom_use( "TPNEW" ) )  THEN
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = ( zpronew(ji,jj,jk) + zpronewd(ji,jj,jk) ) * zfact * tmask(ji,jj,jk)  ! total new production
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "TPNEW"  , zw3d )
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = ( zpronew(ji,jj,jk) + zpronewd(ji,jj,jk) ) * zfact * tmask(ji,jj,jk)  ! total new production
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "TPNEW"  , zw3d )
           ENDIF
           IF( iom_use( "TPBFE" ) )  THEN
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = ( zprofen(ji,jj,jk) + zprofed(ji,jj,jk) ) * zfact * tmask(ji,jj,jk)  ! total biogenic iron production
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "TPBFE"  , zw3d )
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = ( zprofen(ji,jj,jk) + zprofed(ji,jj,jk) ) * zfact * tmask(ji,jj,jk)  ! total biogenic iron production
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "TPBFE"  , zw3d )
           ENDIF
           IF( iom_use( "INTPPPHY" ) .OR. iom_use( "INTPPPHY2" ) ) THEN  
@@ -708,5 +708,5 @@
              ENDDO
 !$OMP END PARALLEL
-            CALL iom_put( "INTPBFE" , zw2d )
+             CALL iom_put( "INTPBFE" , zw2d )
           ENDIF
           IF( iom_use( "INTPBSI" ) ) THEN           !   total biogenic silica production  ( vertically integrated )
@@ -739,19 +739,19 @@
            zfact = 1.e+3 * rfact2r
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     trc3d(ji,jj,jk,jp_pcs0_3d + 4)  = zprorca (ji,jj,jk) * zfact * tmask(ji,jj,jk)
-                     trc3d(ji,jj,jk,jp_pcs0_3d + 5)  = zprorcad(ji,jj,jk) * zfact * tmask(ji,jj,jk)
-                     trc3d(ji,jj,jk,jp_pcs0_3d + 6)  = zpronew (ji,jj,jk) * zfact * tmask(ji,jj,jk)
-                     trc3d(ji,jj,jk,jp_pcs0_3d + 7)  = zpronewd(ji,jj,jk) * zfact * tmask(ji,jj,jk)
-                     trc3d(ji,jj,jk,jp_pcs0_3d + 8)  = zprorcad(ji,jj,jk) * zfact * tmask(ji,jj,jk) * zysopt(ji,jj,jk)
-                     trc3d(ji,jj,jk,jp_pcs0_3d + 9)  = zprofed (ji,jj,jk) * zfact * tmask(ji,jj,jk)
+           DO jk = 1, jpk
+              DO jj = 1, jpj
+                 DO ji = 1, jpi
+                    trc3d(ji,jj,jk,jp_pcs0_3d + 4)  = zprorca (ji,jj,jk) * zfact * tmask(ji,jj,jk)
+                    trc3d(ji,jj,jk,jp_pcs0_3d + 5)  = zprorcad(ji,jj,jk) * zfact * tmask(ji,jj,jk)
+                    trc3d(ji,jj,jk,jp_pcs0_3d + 6)  = zpronew (ji,jj,jk) * zfact * tmask(ji,jj,jk)
+                    trc3d(ji,jj,jk,jp_pcs0_3d + 7)  = zpronewd(ji,jj,jk) * zfact * tmask(ji,jj,jk)
+                    trc3d(ji,jj,jk,jp_pcs0_3d + 8)  = zprorcad(ji,jj,jk) * zfact * tmask(ji,jj,jk) * zysopt(ji,jj,jk)
+                    trc3d(ji,jj,jk,jp_pcs0_3d + 9)  = zprofed (ji,jj,jk) * zfact * tmask(ji,jj,jk)
 #  if ! defined key_kriest
-                     trc3d(ji,jj,jk,jp_pcs0_3d + 10) = zprofen (ji,jj,jk) * zfact * tmask(ji,jj,jk)
+                    trc3d(ji,jj,jk,jp_pcs0_3d + 10) = zprofen (ji,jj,jk) * zfact * tmask(ji,jj,jk)
 #  endif
-                  END DO
-               END DO
-            END DO
+                 END DO
+              END DO
+           END DO
         ENDIF
      ENDIF
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 7525)
@@ -115,16 +115,16 @@
             IF( nn_ice_tr == -1 .AND. .NOT. ln_ironice ) THEN
 !$OMP PARALLEL DO schedule(static) private(jj, ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  dust(ji,jj) = sf_dust(1)%fnow(ji,jj,1)
-               END DO
-            END DO
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     dust(ji,jj) = sf_dust(1)%fnow(ji,jj,1)
+                  END DO
+               END DO
             ELSE
 !$OMP PARALLEL DO schedule(static) private(jj, ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  dust(ji,jj) = sf_dust(1)%fnow(ji,jj,1) * ( 1.0 - fr_i(ji,jj) )
-               END DO
-            END DO
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     dust(ji,jj) = sf_dust(1)%fnow(ji,jj,1) * ( 1.0 - fr_i(ji,jj) )
+                  END DO
+               END DO
             ENDIF
          ENDIF
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90	(revision 7525)
@@ -209,5 +209,6 @@
       !  Exchange between organic matter compartments due to coagulation/disaggregation
       !  ---------------------------------------------------
-!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zstep, zfact, zagg, zagg1, zagg2, zagg3, zagg4, zaggfe, zaggdoc, zaggdoc2, zaggdoc3)
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zstep, zfact, zagg, zagg1, zagg2) &
+!$OMP& private(zagg3, zagg4, zaggfe, zaggdoc, zaggdoc2, zaggdoc3)
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -259,5 +260,5 @@
      ! Total carbon export per year
      IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  ) &
- &        t_oce_co2_exp = glob_sum( ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * e1e2t(:,:) * tmask(:,:,1) )
+        &        t_oce_co2_exp = glob_sum( ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * e1e2t(:,:) * tmask(:,:,1) )
      !
      IF( lk_iomput ) THEN
@@ -269,81 +270,81 @@
           IF( iom_use( "EPC100" ) )  THEN
 !$OMP DO schedule(static) private(jj, ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zw2d(ji,jj) = ( sinking(ji,jj,ik100) + sinking2(ji,jj,ik100) ) * zfact * tmask(ji,jj,1) ! Export of carbon at 100m
-               END DO
-            END DO
-              CALL iom_put( "EPC100"  , zw2d )
+             DO jj = 1, jpj
+                DO ji = 1, jpi
+                   zw2d(ji,jj) = ( sinking(ji,jj,ik100) + sinking2(ji,jj,ik100) ) * zfact * tmask(ji,jj,1) ! Export of carbon at 100m
+                END DO
+             END DO
+             CALL iom_put( "EPC100"  , zw2d )
           ENDIF
           IF( iom_use( "EPFE100" ) )  THEN
 !$OMP DO schedule(static) private(jj, ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zw2d(ji,jj) = ( sinkfer(ji,jj,ik100) + sinkfer2(ji,jj,ik100) ) * zfact * tmask(ji,jj,1) ! Export of iron at 100m
-               END DO
-            END DO
-              CALL iom_put( "EPFE100"  , zw2d )
+             DO jj = 1, jpj
+                DO ji = 1, jpi
+                   zw2d(ji,jj) = ( sinkfer(ji,jj,ik100) + sinkfer2(ji,jj,ik100) ) * zfact * tmask(ji,jj,1) ! Export of iron at 100m
+                END DO
+             END DO
+             CALL iom_put( "EPFE100"  , zw2d )
           ENDIF
           IF( iom_use( "EPCAL100" ) )  THEN
 !$OMP DO schedule(static) private(jj, ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zw2d(ji,jj) = sinkcal(ji,jj,ik100) * zfact * tmask(ji,jj,1) ! Export of calcite at 100m
-               END DO
-            END DO
-              CALL iom_put( "EPCAL100"  , zw2d )
+             DO jj = 1, jpj
+                DO ji = 1, jpi
+                   zw2d(ji,jj) = sinkcal(ji,jj,ik100) * zfact * tmask(ji,jj,1) ! Export of calcite at 100m
+                END DO
+             END DO
+             CALL iom_put( "EPCAL100"  , zw2d )
           ENDIF
           IF( iom_use( "EPSI100" ) )  THEN
 !$OMP DO schedule(static) private(jj, ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zw2d(ji,jj) =  sinksil(ji,jj,ik100) * zfact * tmask(ji,jj,1) ! Export of bigenic silica at 100m
-               END DO
-            END DO
-              CALL iom_put( "EPSI100"  , zw2d )
+             DO jj = 1, jpj
+                DO ji = 1, jpi
+                   zw2d(ji,jj) =  sinksil(ji,jj,ik100) * zfact * tmask(ji,jj,1) ! Export of bigenic silica at 100m
+                END DO
+             END DO
+             CALL iom_put( "EPSI100"  , zw2d )
           ENDIF
           IF( iom_use( "EXPC" ) )  THEN
 !$OMP DO schedule(static) private(jk, jj, ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = ( sinking(ji,jj,jk) + sinking2(ji,jj,jk) ) * zfact * tmask(ji,jj,jk) ! Export of carbon in the water column
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "EXPC"  , zw3d )
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = ( sinking(ji,jj,jk) + sinking2(ji,jj,jk) ) * zfact * tmask(ji,jj,jk) ! Export of carbon in the water column
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "EXPC"  , zw3d )
           ENDIF
           IF( iom_use( "EXPFE" ) )  THEN
 !$OMP DO schedule(static) private(jk, jj, ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = ( sinkfer(ji,jj,jk) + sinkfer2(ji,jj,jk) ) * zfact * tmask(ji,jj,jk) ! Export of iron 
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "EXPFE"  , zw3d )
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = ( sinkfer(ji,jj,jk) + sinkfer2(ji,jj,jk) ) * zfact * tmask(ji,jj,jk) ! Export of iron 
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "EXPFE"  , zw3d )
           ENDIF
           IF( iom_use( "EXPCAL" ) )  THEN
 !$OMP DO schedule(static) private(jk, jj, ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = sinkcal(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of calcite 
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "EXPCAL"  , zw3d )
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = sinkcal(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of calcite 
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "EXPCAL"  , zw3d )
           ENDIF
           IF( iom_use( "EXPSI" ) )  THEN
 !$OMP DO schedule(static) private(jk, jj, ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = sinksil(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of bigenic silica
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "EXPSI"  , zw3d )
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = sinksil(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of bigenic silica
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "EXPSI"  , zw3d )
           ENDIF
           IF( iom_use( "tcexp" ) )  CALL iom_put( "tcexp" , t_oce_co2_exp * zfact )   ! molC/s
@@ -526,5 +527,6 @@
       zval4 = 4. + xkr_eta
 
-!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, znum, zeps, zfm, zsm, zdiv, zdiv1, zdiv2, zdiv3, zdiv4, zdiv5, zagg, zagg1, zagg2, zagg3, zagg4, zagg5, zfract, zaggdoc, zaggdoc1, zaggsh, zaggsi, znumdoc)
+!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, znum, zeps, zfm, zsm, zdiv, zdiv1, zdiv2, zdiv3, zdiv4) &
+!$OMP& private(zdiv5, zagg, zagg1, zagg2, zagg3, zagg4, zagg5, zfract, zaggdoc, zaggdoc1, zaggsh, zaggsi, znumdoc)
       DO jk = 1,jpkm1
          DO jj = 1,jpj
@@ -632,114 +634,114 @@
           IF( iom_use( "EPC100" ) )  THEN
 !$OMP DO schedule(static) private(jj, ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zw2d(ji,jj) = sinking(ji,jj,ik100) * zfact * tmask(ji,jj,1) ! Export of carbon at 100m
-               END DO
-            END DO
-              CALL iom_put( "EPC100"  , zw2d )
+             DO jj = 1, jpj
+                DO ji = 1, jpi
+                   zw2d(ji,jj) = sinking(ji,jj,ik100) * zfact * tmask(ji,jj,1) ! Export of carbon at 100m
+                END DO
+             END DO
+             CALL iom_put( "EPC100"  , zw2d )
           ENDIF
           IF( iom_use( "EPN100" ) )  THEN
 !$OMP DO schedule(static) private(jj, ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zw2d(ji,jj) = sinking2(ji,jj,ik100) * zfact * tmask(ji,jj,1) ! Export of number of aggregates ?
-               END DO
-            END DO
-              CALL iom_put( "EPN100"  , zw2d )
+             DO jj = 1, jpj
+                DO ji = 1, jpi
+                   zw2d(ji,jj) = sinking2(ji,jj,ik100) * zfact * tmask(ji,jj,1) ! Export of number of aggregates ?
+                END DO
+             END DO
+             CALL iom_put( "EPN100"  , zw2d )
           ENDIF
           IF( iom_use( "EPCAL100" ) )  THEN
 !$OMP DO schedule(static) private(jj, ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zw2d(ji,jj) = sinkcal(ji,jj,ik100) * zfact * tmask(ji,jj,1) !Export of calcite at 100m
-               END DO
-            END DO
-              CALL iom_put( "EPCAL100"  , zw2d )
+             DO jj = 1, jpj
+                DO ji = 1, jpi
+                   zw2d(ji,jj) = sinkcal(ji,jj,ik100) * zfact * tmask(ji,jj,1) !Export of calcite at 100m
+                END DO
+             END DO
+             CALL iom_put( "EPCAL100"  , zw2d )
           ENDIF
           IF( iom_use( "EPSI100" ) )  THEN
 !$OMP DO schedule(static) private(jj, ji)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zw2d(ji,jj) =  sinksil(ji,jj,ik100) * zfact * tmask(ji,jj,1) ! Export of bigenic silica at 100m
-               END DO
-            END DO
-              CALL iom_put( "EPSI100"  , zw2d )
+             DO jj = 1, jpj
+                DO ji = 1, jpi
+                   zw2d(ji,jj) =  sinksil(ji,jj,ik100) * zfact * tmask(ji,jj,1) ! Export of bigenic silica at 100m
+                END DO
+             END DO
+             CALL iom_put( "EPSI100"  , zw2d )
           ENDIF
           IF( iom_use( "EXPC" ) )  THEN
 !$OMP DO schedule(static) private(jk, jj, ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = sinking(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of carbon in the water column
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "EXPC"  , zw3d )
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = sinking(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of carbon in the water column
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "EXPC"  , zw3d )
           ENDIF
           IF( iom_use( "EXPN" ) )  THEN
 !$OMP DO schedule(static) private(jk, jj, ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = sinking(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of carbon in the water column
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "EXPN"  , zw3d )
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = sinking(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of carbon in the water column
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "EXPN"  , zw3d )
           ENDIF
           IF( iom_use( "EXPCAL" ) )  THEN
 !$OMP DO schedule(static) private(jk, jj, ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = sinkcal(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of calcite 
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "EXPCAL"  , zw3d )
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = sinkcal(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of calcite 
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "EXPCAL"  , zw3d )
           ENDIF
           IF( iom_use( "EXPSI" ) )  THEN
 !$OMP DO schedule(static) private(jk, jj, ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = sinksil(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of bigenic silica
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "EXPSI"  , zw3d )
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = sinksil(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of bigenic silica
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "EXPSI"  , zw3d )
           ENDIF
           IF( iom_use( "XNUM" ) )  THEN
 !$OMP DO schedule(static) private(jk, jj, ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = znum3d(ji,jj,jk) * tmask(ji,jj,jk) !  Number of particles on aggregats
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "XNUM"  , zw3d )
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = znum3d(ji,jj,jk) * tmask(ji,jj,jk) !  Number of particles on aggregats
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "XNUM"  , zw3d )
           ENDIF
           IF( iom_use( "WSC" ) )  THEN
 !$OMP DO schedule(static) private(jk, jj, ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = wsbio3(ji,jj,jk) * tmask(ji,jj,jk) ! Sinking speed of carbon particles
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "WSC"  , zw3d )
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = wsbio3(ji,jj,jk) * tmask(ji,jj,jk) ! Sinking speed of carbon particles
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "WSC"  , zw3d )
           ENDIF
           IF( iom_use( "WSN" ) )  THEN
 !$OMP DO schedule(static) private(jk, jj, ji)
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zw3d(ji,jj,jk) = wsbio4(ji,jj,jk) * tmask(ji,jj,jk) ! Sinking speed of particles number
-                  END DO
-               END DO
-            END DO
-              CALL iom_put( "WSN"  , zw3d )
+             DO jk = 1, jpk
+                DO jj = 1, jpj
+                   DO ji = 1, jpi
+                      zw3d(ji,jj,jk) = wsbio4(ji,jj,jk) * tmask(ji,jj,jk) ! Sinking speed of particles number
+                   END DO
+                END DO
+             END DO
+             CALL iom_put( "WSN"  , zw3d )
           ENDIF
           !
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/TRP/trcbbl.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/TRP/trcbbl.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/TRP/trcbbl.F90	(revision 7525)
@@ -61,6 +61,6 @@
       IF( l_trdtrc )  THEN
          CALL wrk_alloc( jpi, jpj, jpk, jptra, ztrtrd ) ! temporary save of trends
+!$OMP PARALLEL DO schedule(static) private(jn,jk,jj,ji)
          DO jn = 1, jptra
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
             DO jk = 1, jpk
                DO jj = 1, jpj
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90	(revision 7525)
@@ -76,6 +76,6 @@
       IF( l_trdtrc )  THEN
          CALL wrk_alloc( jpi,jpj,jpk,jptra,   ztrtrd )
+!$OMP PARALLEL DO schedule(static) private(jn,jk,jj,ji)
          DO jn = 1, jptra
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
             DO jk = 1, jpk
                DO jj = 1, jpj
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90	(revision 7525)
@@ -108,6 +108,6 @@
       IF( l_trdtrc )  THEN             ! trends: store now fields before the Asselin filter application
          CALL wrk_alloc( jpi, jpj, jpk, jptra, ztrdt )
+!$OMP PARALLEL DO schedule(static) private(jn,jk,jj,ji)
          DO jn = 1, jptra
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
             DO jk = 1, jpk
                DO jj = 1, jpj
@@ -121,8 +121,12 @@
       !                                ! Leap-Frog + Asselin filter time stepping
       IF( neuler == 0 .AND. kt == nittrc000 ) THEN    ! Euler time-stepping at first time-step (only swap)
+!$OMP PARALLEL DO schedule(static) private(jn,jk,jj,ji)
          DO jn = 1, jptra
-!$OMP PARALLEL DO schedule(static) private(jk)
             DO jk = 1, jpkm1
-               trn(:,:,jk,jn) = tra(:,:,jk,jn)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     trn(ji,jj,jk,jn) = tra(ji,jj,jk,jn)
+                  END DO
+               END DO
             END DO
          END DO
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90	(revision 7525)
@@ -109,6 +109,6 @@
          ELSE                                         ! No restart or restart not found: Euler forward time stepping
            zfact = 1._wp
+!$OMP PARALLEL DO schedule(static) private(jn,jj,ji)
            DO jn = 1, jptra
-!$OMP PARALLEL DO schedule(static) private(jj,ji)
               DO jj = 1, jpj
                  DO ji = 1, jpi
@@ -121,6 +121,6 @@
          IF( ln_top_euler ) THEN
             zfact = 1._wp
+!$OMP PARALLEL DO schedule(static) private(jn,jj,ji)
            DO jn = 1, jptra
-!$OMP PARALLEL DO schedule(static) private(jj,ji)
               DO jj = 1, jpj
                  DO ji = 1, jpi
@@ -131,6 +131,6 @@
          ELSE
             zfact = 0.5_wp
+!$OMP PARALLEL DO schedule(static) private(jn,jj,ji)
            DO jn = 1, jptra
-!$OMP PARALLEL DO schedule(static) private(jj,ji)
               DO jj = 1, jpj
                  DO ji = 1, jpi
@@ -220,11 +220,11 @@
          IF( l_trdtrc ) THEN
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-           DO jk = 1, jpk
-              DO jj = 1, jpj
-                 DO ji = 1, jpi
-                    ztrtrd(ji,jj,jk) = tra(ji,jj,jk,jn) - ztrtrd(ji,jj,jk)
-                 END DO
-              END DO
-           END DO                                      ! online coupling free surface or offline with free surface
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrtrd(ji,jj,jk) = tra(ji,jj,jk,jn) - ztrtrd(ji,jj,jk)
+                  END DO
+               END DO
+            END DO                                      ! online coupling free surface or offline with free surface
             CALL trd_tra( kt, 'TRC', jn, jptra_nsr, ztrtrd )
          END IF
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90	(revision 7525)
@@ -71,6 +71,6 @@
       IF( l_trdtrc )  THEN
          CALL wrk_alloc( jpi, jpj, jpk, jptra, ztrtrd )
+!$OMP PARALLEL DO schedule(static) private(jn,jk,jj,ji)
          DO jn = 1, jptra
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
             DO jk = 1, jpk
                DO jj = 1, jpj
@@ -90,7 +90,11 @@
       IF( l_trdtrc )   THEN                      ! save the vertical diffusive trends for further diagnostics
          DO jn = 1, jptra
-!$OMP PARALLEL DO schedule(static) private(jk)
+!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
             DO jk = 1, jpkm1
-               ztrtrd(:,:,jk,jn) = ( ( tra(:,:,jk,jn) - trb(:,:,jk,jn) ) / r2dttrc ) - ztrtrd(:,:,jk,jn)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztrtrd(ji,jj,jk,jn) = ( ( tra(ji,jj,jk,jn) - trb(ji,jj,jk,jn) ) / r2dttrc ) - ztrtrd(ji,jj,jk,jn)
+                  END DO
+               END DO
             END DO
             CALL trd_tra( kt, 'TRC', jn, jptra_zdf, ztrtrd(:,:,:,jn) )
Index: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 7517)
+++ branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 7525)
@@ -119,15 +119,19 @@
 !$OMP PARALLEL DO schedule(static) private(jk)
       DO jk = 1, jpk
-         cvol(:,:,jk) = e1e2t(:,:) * e3t_n(:,:,jk) * tmask(:,:,jk)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               cvol(ji,jj,jk) = e1e2t(ji,jj) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
+            END DO
+         END DO
       END DO
       IF( lk_degrad ) THEN  
 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               cvol(ji,jj,jk) = cvol(ji,jj,jk) * facvol(ji,jj,jk)    ! degrad option: reduction by facvol
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  cvol(ji,jj,jk) = cvol(ji,jj,jk) * facvol(ji,jj,jk)    ! degrad option: reduction by facvol
+               END DO
             END DO
          END DO
-      END DO
       END IF
       !                                                              ! total volume of the ocean 
@@ -255,6 +259,6 @@
         ENDIF
         !
+!$OMP PARALLEL DO schedule(static) private(jn,jk,jj,ji)
       DO jn = 1, jptra
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
          DO jk = 1, jpk
             DO jj = 1, jpj
@@ -268,6 +272,6 @@
       ENDIF
  
+!$OMP PARALLEL DO schedule(static) private(jn,jk,jj,ji)
       DO jn = 1, jptra
-!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
          DO jk = 1, jpk
             DO jj = 1, jpj
