Index: branches/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- branches/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 4771)
+++ branches/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 5075)
@@ -70,6 +70,4 @@
       REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars
       REAL(wp) ::   ze3ua, ze3va
-      !!----------------------------------------------------------------------
-
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws
       !!----------------------------------------------------------------------
@@ -101,15 +99,14 @@
 
       IF( ln_bfrimp ) THEN
-# if defined key_vectopt_loop
-         DO jj = 1, 1
-            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
-# else
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-# endif
                ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points 
                ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
                avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)
                avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
+               ikbu = miku(ji,jj)       ! ocean top level at u- and v-points 
+               ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
+               IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu)
+               IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv)
             END DO
          END DO
@@ -138,6 +135,6 @@
             ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk)
             va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk)
-         ENDDO
-         ! Add bottom stress due to barotropic component only:
+         END DO
+         ! Add bottom/top stress due to barotropic component only:
          DO jj = 2, jpjm1        
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -148,4 +145,10 @@
                ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua
                va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va
+               ikbu = miku(ji,jj)         ! top ocean level at u- and v-points 
+               ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
+               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu)
+               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv)
+               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua
+               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va
             END DO
          END DO
@@ -166,5 +169,5 @@
                zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
                zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
-               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
+               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
                zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
@@ -194,7 +197,8 @@
       !-----------------------------------------------------------------------
       !
-      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
-         DO jj = 2, jpjm1   
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
+      DO jj = 2, jpjm1   
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jk = miku(ji,jj)+1, jpkm1
                zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
             END DO
@@ -204,17 +208,14 @@
       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1) 
+            ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,miku(ji,jj)) + r_vvl   * fse3u_a(ji,jj,miku(ji,jj)) 
 #if defined key_dynspg_ts
-            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
+            ua(ji,jj,miku(ji,jj)) = ua(ji,jj,miku(ji,jj)) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
                &                                      / ( ze3ua * rau0 ) 
 #else
-            ua(ji,jj,1) = ub(ji,jj,1) + p2dt *(ua(ji,jj,1) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
-               &                                                     / ( fse3u(ji,jj,1) * rau0     ) ) 
-#endif
-         END DO
-      END DO
-      DO jk = 2, jpkm1
-         DO jj = 2, jpjm1   
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+            ua(ji,jj,miku(ji,jj)) = ub(ji,jj,miku(ji,jj)) &
+               &                   + p2dt *(ua(ji,jj,miku(ji,jj)) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
+               &                                  / ( fse3u(ji,jj,miku(ji,jj)) * rau0     ) ) 
+#endif
+            DO jk = miku(ji,jj)+1, jpkm1
 #if defined key_dynspg_ts
                zrhs = ua(ji,jj,jk)   ! zrhs=right hand side
@@ -230,9 +231,5 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
-         END DO
-      END DO
-      DO jk = jpk-2, 1, -1
-         DO jj = 2, jpjm1   
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jk = jpk-2, miku(ji,jj), -1
                ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
             END DO
@@ -292,7 +289,8 @@
       !-----------------------------------------------------------------------
       !
-      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
-         DO jj = 2, jpjm1   
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
+      DO jj = 2, jpjm1   
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jk = mikv(ji,jj)+1, jpkm1        
                zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
             END DO
@@ -302,17 +300,14 @@
       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1) 
+            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,mikv(ji,jj)) + r_vvl   * fse3v_a(ji,jj,mikv(ji,jj)) 
 #if defined key_dynspg_ts            
-            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
+            va(ji,jj,mikv(ji,jj)) = va(ji,jj,mikv(ji,jj)) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
                &                                      / ( ze3va * rau0 ) 
 #else
-            va(ji,jj,1) = vb(ji,jj,1) + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
-               &                                                       / ( fse3v(ji,jj,1) * rau0     )  )
-#endif
-         END DO
-      END DO
-      DO jk = 2, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+            va(ji,jj,mikv(ji,jj)) = vb(ji,jj,mikv(ji,jj)) &
+               &                   + p2dt *(va(ji,jj,mikv(ji,jj)) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
+               &                                                       / ( fse3v(ji,jj,mikv(ji,jj)) * rau0     )  )
+#endif
+            DO jk = mikv(ji,jj)+1, jpkm1
 #if defined key_dynspg_ts
                zrhs = va(ji,jj,jk)   ! zrhs=right hand side
@@ -328,9 +323,5 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
-         END DO
-      END DO
-      DO jk = jpk-2, 1, -1
-         DO jj = 2, jpjm1   
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jk = jpk-2, mikv(ji,jj), -1
                va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
             END DO
@@ -352,17 +343,16 @@
       !! restore bottom layer avmu(v) 
       IF( ln_bfrimp ) THEN
-# if defined key_vectopt_loop
-      DO jj = 1, 1
-         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
-# else
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1
-# endif
-            ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
-            ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
-            avmu(ji,jj,ikbu+1) = 0.e0
-            avmv(ji,jj,ikbv+1) = 0.e0
-         END DO
-      END DO
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1
+              ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
+              ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
+              avmu(ji,jj,ikbu+1) = 0.e0
+              avmv(ji,jj,ikbv+1) = 0.e0
+              ikbu = miku(ji,jj)         ! ocean top level at u- and v-points 
+              ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
+              IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0
+              IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0
+           END DO
+        END DO
       ENDIF
       !
