Index: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
===================================================================
--- branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 5006)
+++ branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 5602)
@@ -106,14 +106,15 @@
       ENDIF
       !
-      zwi(:,:,:) = 0.e0 ; zwz(:,:,:) = 0.e0
+      zwi(:,:,:) = 0.e0 ; 
       !
       !                                                          ! ===========
       DO jn = 1, kjpt                                            ! tracer loop
          !                                                       ! ===========
-         ! 1. Bottom value : flux set to zero
+         ! 1. Bottom and k=1 value : flux set to zero
          ! ----------------------------------
          zwx(:,:,jpk) = 0.e0    ;    zwz(:,:,jpk) = 0.e0
          zwy(:,:,jpk) = 0.e0    ;    zwi(:,:,jpk) = 0.e0
-
+          
+         zwz(:,:,1  ) = 0._wp
          ! 2. upstream advection with initial mass fluxes & intermediate update
          ! --------------------------------------------------------------------
@@ -134,28 +135,36 @@
 
          ! upstream tracer flux in the k direction
+         ! Interior value
+         DO jk = 2, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
+                  zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
+                  zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
          ! Surface value
          IF( lk_vvl ) THEN   
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zwz(ji,jj, mikt(ji,jj) ) = 0.e0                         ! volume variable
-               END DO
-            END DO
+            IF ( ln_isfcav ) THEN
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zwz(ji,jj, mikt(ji,jj) ) = 0.e0          ! volume variable
+                  END DO
+               END DO
+            ELSE
+               zwz(:,:,1) = 0.e0          ! volume variable
+            END IF
          ELSE                
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface 
-               END DO
-            END DO   
+            IF ( ln_isfcav ) THEN
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface 
+                  END DO
+               END DO   
+            ELSE
+               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface
+            END IF
          ENDIF
-         ! Interior value
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               DO jk = mikt(ji,jj)+1, jpkm1
-                  zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
-                  zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
-                  zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) )
-               END DO
-            END DO
-         END DO
 
          ! total advective trend
@@ -184,7 +193,7 @@
          END IF
          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN  
-           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
-           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
+         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
+           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) )
+           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) )
          ENDIF
 
@@ -202,17 +211,22 @@
       
          ! antidiffusive flux on k
-         zwz(:,:,1) = 0.e0         ! Surface value
-         !
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ik=mikt(ji,jj)
-               ! surface value
-               zwz(ji,jj,1:ik) = 0.e0
-               ! Interior value
-               DO jk = mikt(ji,jj)+1, jpkm1                    
+         ! Interior value
+         DO jk = 2, jpkm1                    
+            DO jj = 1, jpj
+               DO ji = 1, jpi
                   zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - zwz(ji,jj,jk)
                END DO
             END DO
          END DO
+         ! surface value
+         IF ( ln_isfcav ) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zwz(ji,jj,mikt(ji,jj)) = 0.e0
+               END DO
+            END DO
+         ELSE
+            zwz(:,:,1) = 0.e0
+         END IF
          CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions
          CALL lbc_lnk( zwz, 'W',  1. )
@@ -250,7 +264,7 @@
          END IF
          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN  
-           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) ) + htr_adv(:)
-           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) ) + str_adv(:)
+         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
+           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) + htr_adv(:)
+           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) + str_adv(:)
          ENDIF
          !
@@ -358,8 +372,4 @@
 
          ! upstream tracer flux in the k direction
-         ! Surface value
-         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0._wp                        ! volume variable
-         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface 
-         ENDIF
          ! Interior value
          DO jk = 2, jpkm1
@@ -372,4 +382,26 @@
             END DO
          END DO
+         ! Surface value
+         IF( lk_vvl ) THEN
+            IF ( ln_isfcav ) THEN
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zwz(ji,jj, mikt(ji,jj) ) = 0.e0          ! volume variable +    isf
+                  END DO
+               END DO
+            ELSE
+               zwz(:,:,1) = 0.e0                              ! volume variable + no isf
+            END IF
+         ELSE
+            IF ( ln_isfcav ) THEN
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface +    isf
+                  END DO
+               END DO
+            ELSE
+               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)                                               ! linear free surface + no isf
+            END IF
+         ENDIF
 
          ! total advective trend
@@ -398,7 +430,7 @@
          END IF
          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN  
-           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
-           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
+         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
+           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) )
+           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) )
          ENDIF
 
@@ -524,7 +556,7 @@
          END IF
          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN  
-           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) ) + htr_adv(:)
-           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) ) + str_adv(:)
+         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
+           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) + htr_adv(:)
+           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) + str_adv(:)
          ENDIF
          !
@@ -580,10 +612,10 @@
          &        paft * tmask + zbig * ( 1._wp - tmask )  )
 
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            DO jk = mikt(ji,jj), jpkm1
-               ikm1 = MAX(jk-1,mikt(ji,jj))
-               z2dtt = p2dt(jk)
-               
+      DO jk = 1, jpkm1
+         ikm1 = MAX(jk-1,1)
+         z2dtt = p2dt(jk)
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+
                ! search maximum in neighbourhood
                zup = MAX(  zbup(ji  ,jj  ,jk  ),   &
