Index: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90
===================================================================
--- branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90	(revision 4431)
+++ branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90	(revision 4432)
@@ -137,5 +137,5 @@
       DO jj = 2, jpjm1               ! add the diffusive trend to the general momentum trends
          DO ji = 2, jpim1
-            DO jk = 1, jpkm1
+            DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
 #else
       DO jj = 2, jpjm1               ! add the diffusive trend to the general momentum trends
@@ -238,238 +238,7 @@
 
 #if defined key_z_first
-      !                               ! ********** !   ! ===============
-!      DO jk = 1, jpkm1                ! First step !   ! Horizontal slab
-         !                            ! ********** !   ! ===============
-
-         ! I.1 Vertical gradient of pu and pv at level jk and jk+1
-         ! -------------------------------------------------------
-         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
-         !                             zdkv(jk=1)=zdkv(jk=2)
-#if 0
-         DO jjf = 1, jpj
-            DO jif = 1, jpi
-
-!!$               jj= jjf
-!!$               ji = jif
-!!$
-!!$               zdk1u = ( pu(ji,jj,jk) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1)
-!!$               zdk1v = ( pv(ji,jj,jk) - pv(ji,jj,jk+1) ) * vmask(ji,jj,jk+1)
-!!$
-!!$               IF( jk == 1 ) THEN
-!!$                  zdku = zdk1u
-!!$                  zdkv = zdk1v
-!!$               ELSE
-!!$                  zdku = ( pu(ji,jj,jk-1) - pu(ji,jj,jk) ) * umask(ji,jj,jk)
-!!$                  zdkv = ( pv(ji,jj,jk-1) - pv(ji,jj,jk) ) * vmask(ji,jj,jk)
-!!$               ENDIF
-!            END DO
-!         END DO
-         !                                -----f-----
-         ! I.2 Horizontal fluxes on U          |
-         ! ------------------------===     t   u   t
-         !                                     |
-         ! i-flux at t-point              -----f-----
-!         DO jj = 1, jpjm1
-!            DO ji = 2, jpi
-!         DO jjf = 1, jpj
-!            DO jif = 1, jpi
-
-               jj = MIN(jjf, jpjm1)
-               ji = MAX(2, jif)
-
-!               zabe1 = e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
-
-!               zmkt  = 1._wp/MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
-!                                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1._wp )
-
-!               zcof1 = -e2t(ji,jj) / MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
-!                                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1._wp )   &
-!                     * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
-
-               IF( jk == 1 )THEN
-                  ziut(ji,jj) = tmask(ji,jj,jk) *   &
-                           (  (e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) )   &
-!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
-                            + (-e2t(ji,jj) / MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
-                                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1._wp )   &
-                     * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ) * &
-( ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji,jj,jk+1) + &
-  ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + &
-!                                       +zdk1u(ji,jj) + zdku (ji-1,jj) )  )
-  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji,jj,jk+1) + &
-  ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) &
-)  )
-               ELSE
-                  ziut(ji,jj) = tmask(ji,jj,jk) *   &
-                           (  (e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) )   &
-!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
-                            + (-e2t(ji,jj) / MAX(umask(ji-1,jj,jk)+umask(ji,jj,jk+1)   &
-                                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk), 1._wp )   &
-                     * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ) * &
-( ( pu(ji  ,jj,jk-1) - pu(ji  ,jj,jk  ) ) * umask(ji,jj,jk) + &
-  ( pu(ji-1,jj,jk  ) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) +     &
-!                                       +zdk1u(ji,jj) + zdku (ji-1,jj) )  )
-  ( pu(ji  ,jj,jk  ) - pu(ji  ,jj,jk+1) ) * umask(ji,jj,jk+1) + &
-  ( pu(ji-1,jj,jk-1) - pu(ji-1,jj,jk  ) ) * umask(ji-1,jj,jk) &
-)  )
-               END IF
-!            END DO
-!         END DO
-
-         ! j-flux at f-point
-!         DO jj = 1, jpjm1
-!            DO ji = 1, jpim1
-!         DO jjf = 1, jpj
-!            DO jif = 1, jpi
-
-               jj = MIN(jjf, jpjm1)
-               ji = MIN(jif, jpim1)
-
-               !zabe2 = e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
-
-               !zmkf  = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
-               !               + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )
-
-               !zcof2 = -e1f(ji,jj) /MAX(umask(ji,jj+1,jk)+umask(ji,jj,jk+1)   &
-               !                       + umask(ji,jj+1,jk+1)+umask(ji,jj,jk), 1. )   &
-               !      * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
-
-               IF(jk == 1)THEN
-                  ! zdku = zdk1u
-                  zjuf(ji,jj) = fmask(ji,jj,jk) *   &
-                           (  (e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) ) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) )   &
-!                            + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
-                             + (-e1f(ji,jj) /MAX(umask(ji,jj+1,jk)+umask(ji,jj,jk+1)   &
-                                      + umask(ji,jj+1,jk+1)+umask(ji,jj,jk), 1. )   &
-                     * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) ) * &
-( &
-             ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1)  &
-           + ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj  ,jk+1)  &
-!                                       +zdk1u(ji,jj+1) + zdku (ji,jj) )  )
-           + ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) &
-           + ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1)   &
-)  )
-               
-            ELSE
-               zjuf(ji,jj) = fmask(ji,jj,jk) *   &
-                           (  (e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) ) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) )   &
-!                            + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
-                             + (-e1f(ji,jj) /MAX(umask(ji,jj+1,jk)+umask(ji,jj,jk+1)   &
-                                      + umask(ji,jj+1,jk+1)+umask(ji,jj,jk), 1. )   &
-                     * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) ) * &
-(  &
-  (pu(ji,jj+1,jk-1) - pu(ji,jj+1,jk  ) ) * umask(ji,jj+1,jk) +   &
-  (pu(ji,jj  ,jk  ) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1) +   &
-!                                       +zdk1u(ji,jj+1) + zdku (ji,jj) )  )
-  (pu(ji,jj+1,jk  ) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) + &
-  (pu(ji,jj  ,jk-1) - pu(ji,jj  ,jk  ) ) * umask(ji,jj,jk) &
- )  )
-            ENDIF
-!            END DO
-!         END DO
-
-         !                                 |   t   |
-         ! I.3 Horizontal fluxes on V      |       |
-         ! ------------------------===     f---v---f
-         !                                 |       |
-         ! i-flux at f-point               |   t   |
-!         DO jj = 1, jpjm1
-!            DO ji = 1, jpim1
-!         DO jjf = 1, jpj
-!            DO jif = 1, jpi
-
-               jj = MIN(jjf, jpjm1)
-               ji = MIN(jif, jpim1)
-
-!               zabe1 = e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
-
-!               zmkf  = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
-!                              + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. )
-
-!               zcof1 = (-e2f(ji,jj) / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
-!                              + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. ))   &
-!                     * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
-
-               IF (jk == 1)THEN
-                  ! zdku == zdk1u
-                  zivf(ji,jj) = fmask(ji,jj,jk) *   &
-                           (  (e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) )   &
-!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj)     &
-                             + ((-e2f(ji,jj) / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
-                              + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. ))   &
-                     * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * ( &
-  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) + &
-  ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) +    &
-!                                       +zdk1u(ji,jj) + zdku (ji+1,jj) )  ) 
-  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) +  &
-  ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) ) )
-               ELSE
-                  zivf(ji,jj) = fmask(ji,jj,jk) *   &
-                           (  (e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) )   &
-!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj)     &
-+ ((-e2f(ji,jj) / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
-                              + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. ))   &
-                     * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * ( &
-  ( pu(ji  ,jj,jk-1) - pu(ji  ,jj,jk  ) ) * umask(ji  ,jj,jk  ) + &
-  ( pu(ji+1,jj,jk  ) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) +    &
-!                                       +zdk1u(ji,jj) + zdku (ji+1,jj) )  )
-  ( pu(ji  ,jj,jk  ) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) + &
-  ( pu(ji+1,jj,jk-1) - pu(ji+1,jj,jk  ) ) * umask(ji+1,jj,jk  ) )  )
-               END IF
-!            END DO
-!         END DO
-
-         ! j-flux at t-point
-!         DO jj = 2, jpj
-!            DO ji = 1, jpim1
-!         DO jjf = 1, jpj
-!            DO jif = 1, jpi
-
-               jj = MAX(2,jjf)
-               ji = MIN(jif, jpim1)
-
-               !zabe2 = e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
-
-               !zmkt  = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
-               !               + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
-
-               !zcof2 = (-e1t(ji,jj)/MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
-               !               + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) )   &
-               !      * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
-
-               IF( jk == 1 )THEN
-                  zjvt(ji,jj) = tmask(ji,jj,jk) *   &
-                           (  (e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) )   &
-!                            + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj)     &
-                            + ((-e1t(ji,jj)/MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
-                              + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) )   &
-                     * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * ( &
-  ( pu(ji,jj-1,jk) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + &
-  ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1)   +  &
-!                                       +zdk1u(ji,jj-1) + zdku (ji,jj) )  )
-  ( pu(ji,jj-1,jk) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + &
-  ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1)  )  )
-               ELSE
-                  zjvt(ji,jj) = tmask(ji,jj,jk) *   &
-                           (  (e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) )   &
-!                            + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj)     &
-                            + ((-e1t(ji,jj)/MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
-                              + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) )   &
-                     * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * &
-( &
-  ( pu(ji,jj-1,jk-1) - pu(ji,jj-1,jk  ) ) * umask(ji,jj-1,jk) + &
-  ( pu(ji,jj  ,jk  ) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1)  +    &
-!                                       +zdk1u(ji,jj-1) + zdku (ji,jj) )  )
-  ( pu(ji,jj-1,jk  ) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + &
-  ( pu(ji,jj  ,jk-1) - pu(ji,jj  ,jk  ) ) * umask(ji,jj,jk) )  )
-               END IF
-            END DO
-         END DO
-#endif
-
-         ! I.4 Second derivative (divergence) (not divided by the volume)
-         ! ---------------------
-
-
+      !                               ! ********** ! 
+                                      ! First step ! 
+         !                            ! ********** ! 
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
@@ -590,5 +359,5 @@
 
 
-               DO jk = 2, jpkm1
+               DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
 
 !               plu(ji,jj,jk) = ziut (ji+1,jj) - &
@@ -715,7 +484,7 @@
          END DO
 
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
+         !                                       
+      END DO
+
 #else
       !                               ! ********** !   ! ===============
@@ -850,5 +619,5 @@
 #if defined key_z_first
          DO ji = 2, jpi
-            DO jk = 1, jpk
+            DO jk = 1, mbkmax(ji,jj) ! jpk
                ! i-gradient of u at jj
                zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( pu(ji,jj  ,jk) - pu(ji-1,jj  ,jk) )
@@ -862,5 +631,5 @@
          END DO
          DO ji = 1, jpim1
-            DO jk = 1, jpk
+            DO jk = 1, mbkmax(ji,jj) ! jpk
                ! i-gradient of v at jj
                zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( pv(ji+1,jj,jk) - pv(ji  ,jj  ,jk) )
@@ -893,14 +662,22 @@
          ! Surface and bottom vertical fluxes set to zero
 
+#if defined key_z_first
+         DO ji=1, jpi
+            zfuw(ji, 1               ) = 0.e0
+            zfvw(ji, 1               ) = 0.e0
+            zfuw(ji,mbkmax(ji,jj):jpk) = 0.e0
+            zfvw(ji,mbkmax(ji,jj):jpk) = 0.e0
+         END DO
+#else
          zfuw(:, 1 ) = 0.e0
          zfvw(:, 1 ) = 0.e0
          zfuw(:,jpk) = 0.e0
          zfvw(:,jpk) = 0.e0
-
+#endif
          ! interior (2=<jk=<jpk-1) on pu field
 
 #if defined key_z_first
          DO ji = 2, jpim1
-            DO jk = 2, jpkm1
+            DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
 #else
          DO jk = 2, jpkm1
@@ -935,5 +712,5 @@
 #if defined key_z_first
          DO ji = 2, jpim1
-            DO jk = 2, jpkm1
+            DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
 #else
          DO jk = 2, jpkm1
@@ -972,5 +749,5 @@
 #if defined key_z_first
             DO ji = 2, jpim1
-               DO jk = 1, jpkm1
+               DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
 #else
             DO jk = 1, jpkm1
@@ -992,5 +769,5 @@
 #if defined key_z_first
             DO ji = 2, jpim1
-               DO jk = 1, jpkm1
+               DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
 #else
             DO jk = 1, jpkm1
