Index: /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
===================================================================
--- /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 4426)
+++ /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 4427)
@@ -203,4 +203,24 @@
 #endif
 !!gm end
+
+!! smp FINISS 
+#if defined key_z_first
+      ! Reset mbkmax to be the first level for which tmask is zero
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            mbkmax(ji,jj) = 1 
+            kloop: DO jk = 1, jpkorig
+               IF (tmask(ji,jj,jk) == 0.0_wp) THEN
+                  mbkmax(ji,jj) = jk
+                  EXIT kloop
+               END IF 
+            END DO kloop
+            IF (mbkmax(ji,jj) > jpk) THEN
+               WRITE (*,*) 'FINISS error: mbkmax(',ji,',',jj,') > jpk (',jpk,') on subdomain ',narea
+               mbkmax(ji,jj) = jpk
+            END IF
+         END DO
+      END DO
+#endif
 
       ! Interior domain mask (used for global sum)
Index: /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 4426)
+++ /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 4427)
@@ -131,5 +131,5 @@
       ! -----------------------------------
       IF( lzoom       )   CALL zgr_bat_zoom     ! correct mbathy in case of zoom subdomain
-      IF( .NOT.lk_c1d )   CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isoated ocean points
+      IF( .NOT.lk_c1d )   CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points
                           CALL zgr_bot_level    ! deepest ocean level for t-, u- and v-points
       !
@@ -819,6 +819,12 @@
       jpkorig = jpk
       IF( domtrim_z ) THEN
-         mbkmax(:,:) = MAX(mbkt(:,:)+1, mbku(:,:), mbkv(:,:))
-!         jpkf = MIN(jpk, 1 + MAXVAL( mbkmax(:,:) ) )
+!        mbkmax(:,:) = MAX(mbkt(:,:)+1, mbku(:,:), mbkv(:,:))
+!        write(*,*) narea, ': SMPDBG: ji, jj, mbkt(ji,jj), mbku(ji,jj), mbkv(ji,jj), mbkmax(ji,jj)'
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               mbkmax(ji,jj) = MIN(jpk, MAX(mbkt(ji,jj)+1, mbku(ji,jj), mbkv(ji,jj)))
+               ! write(*,*) narea, ': SMPDBG: ', ji, jj, mbkt(ji,jj), mbku(ji,jj), mbkv(ji,jj), mbkmax(ji,jj)
+            END DO
+         END DO
          jpkf = MAXVAL( mbkmax(:,:) )
          WRITE(*,*) narea,': ARPDBG: shallowest pt and jpkf = ', &
@@ -830,4 +836,5 @@
       ELSE
          WRITE(*,*) narea,': ARPDBG: NOT trimming domain in z'
+         mbkmax(:,:) = jpk
          jpkf = jpk
       END IF
Index: /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
===================================================================
--- /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 4426)
+++ /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 4427)
@@ -87,8 +87,8 @@
       DO jj = 2, jpj            ! Horizontal kinetic energy at T-point
          DO ji = 2, jpi
-            DO jk = 1, jpkm1
+            DO jk = 1, mbkmax(ji,jj)-1
                zhke(ji,jj,jk) = 0.25 * (   un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
                   &                      + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)   &
-                                         + vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
+                  &                      + vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
                   &                      + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)   )
             END DO  
@@ -97,5 +97,5 @@
       DO jj = 2, jpjm1          ! add the gradient of kinetic energy to the general momentum trends
          DO ji = 2, jpim1
-            DO jk = 1, jpkm1
+            DO jk = 1, mbkmax(ji,jj)-1
                ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
                va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
Index: /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap.F90
===================================================================
--- /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap.F90	(revision 4426)
+++ /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap.F90	(revision 4427)
@@ -83,5 +83,5 @@
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
-            DO jk = 1, jpkm1
+            DO jk = 1, mbkmax(ji,jj)-1
 #else
       !                                                ! ===============
Index: /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
===================================================================
--- /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 4426)
+++ /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 4427)
@@ -99,5 +99,5 @@
             zwuw(ji,jj, 1 ) = 0.e0     ! Surface values set to zero
             zwvw(ji,jj, 1 ) = 0.e0
-            DO jk = 2, jpkm1
+            DO jk = 2, mbkmax(ji,jj)-1
                zwuw(ji,jj,jk) =   ( zww(ji+1,jj  )*wn(ji+1,jj  ,jk) + zww(ji,jj)*wn(ji,jj,jk) )   &
                   &             * ( un(ji,jj,jk-1)-un(ji,jj,jk) ) 
@@ -105,6 +105,6 @@
                   &             * ( vn(ji,jj,jk-1)-vn(ji,jj,jk) )
             END DO  
-            zwuw(ji,jj,jpk) = 0.e0     ! Bottom values set to zero
-            zwvw(ji,jj,jpk) = 0.e0
+            zwuw(ji,jj,mbkmax(ji,jj)) = 0.e0     ! Bottom values set to zero
+            zwvw(ji,jj,mbkmax(ji,jj)) = 0.e0
          END DO   
       END DO
@@ -136,5 +136,5 @@
       DO jj = 2, jpjm1              ! Vertical momentum advection at u- and v-points
          DO ji = 2, jpim1
-            DO jk = 1, jpkm1
+            DO jk = 1, mbkmax(ji,jj)-1
 #else
       DO jk = 1, jpkm1              ! Vertical momentum advection at u- and v-points
Index: /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 4426)
+++ /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 4427)
@@ -94,4 +94,7 @@
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
+#if defined key_z_first
+      INTEGER  ::   klim         ! upper bound on k loop
+#endif
       REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars
       !!----------------------------------------------------------------------
@@ -107,5 +110,15 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
          !
+#if defined key_z_first
+         DO jj=1,jpj
+            DO ji=1,jpi
+               DO jk=mbkmax(ji,jj), jpk
+                  wn(ji,jj,jk) = 0._wp
+               END DO
+            END DO
+         END DO
+#else
          wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
+#endif
          !
          IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only)
@@ -161,76 +174,19 @@
             ! gdept_1(1:jpkm1,:,:) = (gdept(1:jpkm1,:,:)*(1.+sshn(:,:)*mut(1:jpkm1,:,:))) 
             ! which contains non-conforming array expressions.
-         DO jj=1,jpj,1
-            DO ji=1,jpi,1
-               DO jk=1,jpk,1
-                  fsdept(ji,jj,jk) = fsdept_n(ji,jj,jk)   ! now local depths stored in fsdep. arrays
-               END DO
-            END DO
-         END DO
-         DO jj=1,jpj,1
-            DO ji=1,jpi,1
-               DO jk=1,jpk,1
-                  fsdepw(ji,jj,jk) = fsdepw_n(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         DO jj=1,jpj,1
-            DO ji=1,jpi,1
-               DO jk=1,jpk,1
-                  fsde3w(ji,jj,jk) = fsde3w_n(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-            !
-         DO jj=1,jpj,1
-            DO ji=1,jpi,1
-               DO jk=1,jpk,1
-                  fse3t (ji,jj,jk) = fse3t_n (ji,jj,jk)   ! vertical scale factors stored in fse3. arrays
-               END DO
-            END DO
-         END DO
-         DO jj=1,jpj,1
-            DO ji=1,jpi,1
-               DO jk=1,jpk,1
-                  fse3u (ji,jj,jk) = fse3u_n (ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         DO jj=1,jpj,1
-            DO ji=1,jpi,1
-               DO jk=1,jpk,1
-                  fse3v (ji,jj,jk) = fse3v_n (ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         DO jj=1,jpj,1
-            DO ji=1,jpi,1
-               DO jk=1,jpk,1
-                  fse3f (ji,jj,jk) = fse3f_n (ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         DO jj=1,jpj,1
-            DO ji=1,jpi,1
-               DO jk=1,jpk,1
-                  fse3w (ji,jj,jk) = fse3w_n (ji,jj,jk)
-               END DO
-            END DO
-         END DO
-
-
-         DO jj=1,jpj,1
-            DO ji=1,jpi,1
-               DO jk=1,jpk,1
-                  fse3uw(ji,jj,jk) = fse3uw_n(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-
-         DO jj=1,jpj,1
-            DO ji=1,jpi,1
-               DO jk=1,jpk,1
-                  fse3vw(ji,jj,jk) = fse3vw_n(ji,jj,jk)
-               END DO
+         DO jj=1,jpj
+            DO ji=1,jpi
+               klim=mbkmax(ji,jj)
+               ! now local depths stored in fsdep. arrays
+               fsdept(ji,jj,1:klim) = fsdept_n(ji,jj,1:klim)
+               fsdepw(ji,jj,1:klim) = fsdepw_n(ji,jj,1:klim)
+               fsde3w(ji,jj,1:klim) = fsde3w_n(ji,jj,1:klim)
+               ! vertical scale factors stored in fse3. arrays
+               fse3t (ji,jj,1:klim) = fse3t_n (ji,jj,1:klim)
+               fse3u (ji,jj,1:klim) = fse3u_n (ji,jj,1:klim)
+               fse3v (ji,jj,1:klim) = fse3v_n (ji,jj,1:klim)
+               fse3f (ji,jj,1:klim) = fse3f_n (ji,jj,1:klim)
+               fse3w (ji,jj,1:klim) = fse3w_n (ji,jj,1:klim)
+               fse3uw(ji,jj,1:klim) = fse3uw_n(ji,jj,1:klim)
+               fse3vw(ji,jj,1:klim) = fse3vw_n(ji,jj,1:klim)
             END DO
          END DO
@@ -279,5 +235,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            DO jk = 1, jpkm1                           ! Horizontal divergence of barotropic transports
+            DO jk = 1, mbkmax(ji,jj)-1                 ! Horizontal divergence of barotropic transports
                zhdiv(ji,jj) = zhdiv(ji,jj) + fse3t(ji,jj,jk) * hdivn(ji,jj,jk)
             END DO
@@ -355,5 +311,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            DO jk = jpkm1, 1, -1                      ! integrate from the bottom the hor. divergence
+            DO jk = mbkmax(ji,jj)-1, 1, -1                      ! integrate from the bottom the hor. divergence
                 wn(ji,jj,jk) = wn(ji,jj,jk+1)                               &
                    &         -   fse3t_n(ji,jj,jk) * hdivn(ji,jj,jk)        &
@@ -390,5 +346,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               DO jk = 1, jpk
+               DO jk = 1, mbkmax(ji,jj)
                   z3d(ji,jj,jk) = wn(ji,jj,jk) * z2d(ji,jj)
                END DO
Index: /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 4426)
+++ /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 4427)
@@ -815,6 +815,11 @@
       !
       !                          !==   surface mixed layer mask   !
+#if defined key_z_first
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            DO jk = 1, jpkf          ! =1 inside the mixed layer, =0 otherwise
+#else
       DO jk = 1, jpkf                ! =1 inside the mixed layer, =0 otherwise
-# if ( defined key_vectopt_loop ) && ! ( defined key_z_first )
+# if defined key_vectopt_loop 
          DO jj = 1, 1
             DO ji = 1, jpij   ! vector opt. (forced unrolling)
@@ -823,4 +828,5 @@
             DO ji = 1, jpi
 # endif
+#endif
                ik = nmln(ji,jj) - 1
                IF( jk <= ik ) THEN   ;   omlmask(ji,jj,jk) = 1._wp
@@ -967,5 +973,5 @@
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
-                  DO jk = 1, jpkf
+                  DO jk = 1, mbkmax(ji,jj)
 #else
             DO jk = 1, jpkf
Index: /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
===================================================================
--- /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 4426)
+++ /branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 4427)
@@ -146,4 +146,5 @@
 
 #if defined key_z_first
+!! SMP Should not need to reserve or release 9 and 10 any more. 
       IF( wrk_in_use(3, 6,7,8,9,10) .OR. wrk_in_use(2, 3) ) THEN
 #else
@@ -162,4 +163,10 @@
       !                                                          ! ===========
 !DIR$ SHORTLOOP
+
+#if defined key_z_first
+      zdit(:,:,:) = 0.0_wp
+      zdjt(:,:,:) = 0.0_wp
+#endif
+
       DO jn = 1, kjpt                                            ! tracer loop
          !                                                       ! ===========
@@ -169,30 +176,23 @@
          !!----------------------------------------------------------------------
          !CALL timing_start('traldf_iso_I')
+
+         ! Horizontal tracer gradient 
+#if defined key_z_first
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
+                  zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
+                  zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+#else
          !!bug ajout.... why?   ( 1,jpj,:) and (jpi,1,:) should be sufficient....
-#if defined key_z_first
-         DO jj=1,jpj,1
-            DO jk=1,jpkf,1
-               zdit(1  ,jj,jk) = 0.0_wp
-               zdit(jpi,jj,jk) = 0.0_wp
-               zdjt(1  ,jj,jk) = 0.0_wp
-               zdjt(jpi,jj,jk) = 0.0_wp
-            END DO
-         END DO
-#else
          zdit (1,:,1:jpkf) = 0.e0     ;     zdit (jpi,:,1:jpkf) = 0.e0
          zdjt (1,:,1:jpkf) = 0.e0     ;     zdjt (jpi,:,1:jpkf) = 0.e0
-#endif
          !!end
-
-         ! Horizontal tracer gradient 
-#if defined key_z_first
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               DO jk = 1, jpkfm1 ! jpkm1
-#else
          DO jk = 1, jpkfm1 ! jpkm1
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-#endif
                   zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
                   zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
@@ -200,4 +200,5 @@
             END DO
          END DO
+#endif
          IF( ln_zps ) THEN      ! partial steps correction at the last ocean level 
             DO jj = 1, jpjm1
@@ -257,5 +258,5 @@
          DO jj = 2 , jpjm1
             DO ji = 2, jpim1
-               DO jk = 1, jpkfm1 ! jpkm1
+               DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
 
                   ! 1. Vertical tracer gradient at level jk and jk+1
@@ -392,5 +393,5 @@
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
-                  DO jk = 1, jpkfm1 ! jpkm1
+                  DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
 #else
             DO jk = 1, jpkfm1 ! jpkm1
@@ -409,5 +410,5 @@
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
-                  DO jk = 1, jpkfm1 ! jpkm1
+                  DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
 #else
             DO jk = 1, jpkfm1 ! jpkm1
@@ -433,14 +434,4 @@
          ! Local constant initialization
          ! -----------------------------
-#if defined key_z_first
-         DO jj=1,jpj,1
-            DO jk=1,jpkf,1
-               ztfw(1  ,jj,jk) = 0.0_wp
-               ztfw(jpi,jj,jk) = 0.0_wp
-            END DO
-         END DO
-#else
-         ztfw(1,:,1:jpkf) = 0.e0     ;     ztfw(jpi,:,1:jpkf) = 0.e0
-#endif
          ! Vertical fluxes
          ! ---------------
@@ -448,12 +439,7 @@
          ! Surface and bottom vertical fluxes set to zero
 #if defined key_z_first
-         DO ji=1,jpi,1
-            DO jj=1,jpj,1
-               ztfw(ji,jj,1  ) = 0.0_wp
-               ztfw(ji,jj,jpkf) = 0.0_wp ! ARPDBG - should this be jpk anyway
-                                         ! since may be below ocean floor?
-            END DO
-         END DO
-#else
+         ztfw(:,:,:) = 0.0_wp
+#else
+         ztfw(1,:,1:jpkf) = 0.e0     ;     ztfw(jpi,:,1:jpkf) = 0.e0
          ztfw(:,:, 1 ) = 0.e0      ;      ztfw(:,:,jpkf) = 0.e0
 #endif
@@ -463,5 +449,5 @@
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               DO jk = 2, jpkfm1
+               DO jk = 2, mbkmax(ji,jj)-1
 #else
          DO jk = 2, jpkfm1
@@ -493,5 +479,5 @@
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               DO jk = 1, jpkfm1
+               DO jk = 1, mbkmax(ji,jj)-1
 #else
          DO jk = 1, jpkfm1
