Index: branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 4750)
+++ branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 5965)
@@ -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,11 +99,6 @@
 
       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)
@@ -114,4 +107,14 @@
             END DO
          END DO
+         IF ( ln_isfcav ) THEN
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  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
+         END IF
       ENDIF
 
@@ -138,6 +141,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.
@@ -150,4 +153,16 @@
             END DO
          END DO
+         IF ( ln_isfcav ) THEN
+            DO jj = 2, jpjm1        
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  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
+         END IF
       ENDIF
 #endif
@@ -164,9 +179,9 @@
                ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl   * fse3u_a(ji,jj,jk)   ! after scale factor at T-point
                zcoef = - p2dt / ze3ua      
-               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)
-               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
-               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
+               zzwi          = zcoef * avmu  (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
+               zwi(ji,jj,jk) = zzwi  * wumask(ji,jj,jk  )
+               zzws          = zcoef * avmu  (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
+               zws(ji,jj,jk) = zzws  * wumask(ji,jj,jk+1)
+               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
             END DO
          END DO
@@ -194,27 +209,29 @@
       !-----------------------------------------------------------------------
       !
-      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.
-               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
-            END DO
-         END DO
-      END DO
-      !
-      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) 
-#if defined key_dynspg_ts
-            ua(ji,jj,1) = ua(ji,jj,1) + 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
+      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
       DO jk = 2, jpkm1
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
+               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      !
+      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+#if defined key_dynspg_ts
+            ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1) 
+            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
+               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
+#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     ) * umask(ji,jj,1) ) 
+#endif
+         END DO
+      END DO
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
 #if defined key_dynspg_ts
                zrhs = ua(ji,jj,jk)   ! zrhs=right hand side
@@ -233,6 +250,6 @@
       END DO
       DO jk = jpk-2, 1, -1
-         DO jj = 2, jpjm1   
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
                ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
             END DO
@@ -263,8 +280,8 @@
                zcoef = - p2dt / ze3va
                zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
-               zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk)
+               zwi(ji,jj,jk) =  zzwi * wvmask(ji,jj,jk)
                zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
-               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
-               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
+               zws(ji,jj,jk) =  zzws * wvmask(ji,jj,jk+1)
+               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
             END DO
          END DO
@@ -292,5 +309,6 @@
       !-----------------------------------------------------------------------
       !
-      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
+      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
+      DO jk = 2, jpkm1        
          DO jj = 2, jpjm1   
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -302,10 +320,11 @@
       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
          DO ji = fs_2, fs_jpim1   ! vector opt.
+#if defined key_dynspg_ts            
             ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1) 
-#if defined key_dynspg_ts            
             va(ji,jj,1) = va(ji,jj,1) + 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) )   &
+            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
@@ -331,6 +350,6 @@
       END DO
       DO jk = jpk-2, 1, -1
-         DO jj = 2, jpjm1   
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
                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 +371,22 @@
       !! 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
+           END DO
+        END DO
+        IF (ln_isfcav) THEN
+           DO jj = 2, jpjm1
+              DO ji = 2, jpim1
+                 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
+        END IF
       ENDIF
       !
