Index: /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DYN/dynzdf.F90
===================================================================
--- /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DYN/dynzdf.F90	(revision 10007)
+++ /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DYN/dynzdf.F90	(revision 10008)
@@ -71,5 +71,5 @@
       REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars
       REAL(wp) ::   zzws, ze3va        !   -      -
-      REAL(wp) ::   z1_e3un, z1_e3vn   !   -      -
+      REAL(wp) ::   z1_e3ua, z1_e3va   !   -      -
       REAL(wp) ::   zWu , zWv          !   -      -
       REAL(wp) ::   zWui, zWvi         !   -      -
@@ -166,14 +166,15 @@
                   DO jj = 2, jpjm1 
                      DO ji = fs_2, fs_jpim1   ! vector opt.
-                        z1_e3un = 1._wp / e3u_n(ji,jj,jk)
-                        zzwi = ( ( avm   (ji+1,jj,jk  ) + avm   (ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
-                           &     / e3uw_a(ji  ,jj,jk  ) ) * z1_e3un * wumask(ji,jj,jk  )
-                        zzws = ( ( avm   (ji+1,jj,jk+1) + avm   (ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
-                           &     / e3uw_a(ji  ,jj,jk+1) ) * z1_e3un * wumask(ji,jj,jk+1)
+                        ze3ua   = e3u_a(ji,jj,jk)
+                        z1_e3ua = 1._wp / ze3ua
+                        zzwi = - zdt * ( avm   (ji+1,jj,jk  ) + avm   (ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
+                           &     / ( ze3ua * e3uw_n(ji  ,jj,jk  ) ) * wumask(ji,jj,jk  )
+                        zzws = - zdt * ( avm   (ji+1,jj,jk+1) + avm   (ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
+                           &     / ( ze3ua * e3uw_n(ji  ,jj,jk+1) ) * wumask(ji,jj,jk+1)
                         zWu  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji+1,jj,jk  )   &
                            &               + wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1)   )
-                        zwi(ji,jj,jk) = - 1._wp * ( zzwi + MIN( zWu, 0._wp ) * z1_e3un )
-                        zws(ji,jj,jk) = - 1._wp * ( zzws - MAX( zWu, 0._wp ) * z1_e3un )
-                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWu, 0._wp ) + MIN( zWu, 0._wp ) ) * z1_e3un )
+                        zwi(ji,jj,jk) = zzwi - zdt * MIN( zWu, 0._wp ) * z1_e3ua
+                        zws(ji,jj,jk) = zzws + zdt * MAX( zWu, 0._wp ) * z1_e3ua
+                        zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWu, 0._wp ) + MIN( zWu, 0._wp ) ) * z1_e3ua
                      END DO
                   END DO
@@ -183,14 +184,15 @@
                   DO jj = 2, jpjm1 
                      DO ji = fs_2, fs_jpim1   ! vector opt.
-                        z1_e3un = 1._wp / e3u_n(ji,jj,jk)
-                        zzwi = ( ( avm   (ji+1,jj,jk  ) + avm(ji,jj,jk  ) )   &
-                           &     / e3uw_a(ji  ,jj,jk  ) ) * z1_e3un * wumask(ji,jj,jk  )
-                        zzws = ( ( avm   (ji+1,jj,jk+1) + avm(ji,jj,jk+1) )   &
-                           &     / e3uw_a(ji  ,jj,jk+1) ) * z1_e3un * wumask(ji,jj,jk+1)
+                        ze3ua   = e3u_a(ji,jj,jk)
+                        z1_e3ua = 1._wp / ze3ua
+                        zzwi = - zdt * ( avm   (ji+1,jj,jk  ) + avm(ji,jj,jk  ) )   &
+                           &     / ( ze3ua * e3uw_n(ji  ,jj,jk  ) ) * wumask(ji,jj,jk  )
+                        zzws = - zdt * ( avm   (ji+1,jj,jk+1) + avm(ji,jj,jk+1) )   &
+                           &     / ( ze3ua * e3uw_n(ji  ,jj,jk+1) ) * wumask(ji,jj,jk+1)
                         zWu  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji+1,jj,jk  )   &
                            &               + wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1)   )
-                        zwi(ji,jj,jk) = - 1._wp * ( zzwi + MIN( zWu, 0._wp ) * z1_e3un )
-                        zws(ji,jj,jk) = - 1._wp * ( zzws - MAX( zWu, 0._wp ) * z1_e3un )
-                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWu, 0._wp ) + MIN( zWu, 0._wp ) ) * z1_e3un )
+                        zwi(ji,jj,jk) = zzwi - zdt * MIN( zWu, 0._wp ) * z1_e3ua 
+                        zws(ji,jj,jk) = zzws + zdt * MAX( zWu, 0._wp ) * z1_e3ua 
+                        zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWu, 0._wp ) + MIN( zWu, 0._wp ) ) * z1_e3ua
                      END DO
                   END DO
@@ -204,13 +206,13 @@
                      DO ji = fs_2, fs_jpim1   ! vector opt.
                         ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
-                        zzwi = ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
-                           &         / ( ze3ua * e3uw_a(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
-                        zzws = ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
-                           &         / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+                        zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
+                           &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
+                        zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
+                           &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
                         zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) )
                         zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
-                        zwi(ji,jj,jk) = - zdt * ( zzwi +  MIN( zWui, 0._wp ) ) 
-                        zws(ji,jj,jk) = - zdt * ( zzws -  MAX( zWus, 0._wp ) )
-                        zwd(ji,jj,jk) = 1._wp + zzwi + zzws - zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
+                        zwi(ji,jj,jk) = zzwi - zdt * MIN( zWui, 0._wp ) 
+                        zws(ji,jj,jk) = zzws + zdt * MAX( zWus, 0._wp )
+                        zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
                      END DO
                   END DO
@@ -221,11 +223,11 @@
                      DO ji = fs_2, fs_jpim1   ! vector opt.
                         ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
-                        zzwi = ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_a(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
-                        zzws = ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+                        zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
+                        zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
                         zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) )
                         zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
-                        zwi(ji,jj,jk) = - zdt * ( zzwi +  MIN( zWui, 0._wp ) ) 
-                        zws(ji,jj,jk) = - zdt * ( zzws -  MAX( zWus, 0._wp ) )
-                        zwd(ji,jj,jk) = 1._wp + zzwi + zzws - zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
+                        zwi(ji,jj,jk) = zzwi - zdt * MIN( zWui, 0._wp )
+                        zws(ji,jj,jk) = zzws + zdt * MAX( zWus, 0._wp )
+                        zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
                      END DO
                   END DO
@@ -240,11 +242,11 @@
                   DO ji = fs_2, fs_jpim1   ! vector opt.
                      ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
-                     zzwi = ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
-                        &         / ( ze3ua * e3uw_a(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
-                     zzws = ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
-                        &         / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
-                     zwi(ji,jj,jk) = - zdt * zzwi
-                     zws(ji,jj,jk) = - zdt * zzws
-                     zwd(ji,jj,jk) = 1._wp + zzwi + zzws
+                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
+                        &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
+                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
+                        &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+                     zwi(ji,jj,jk) = zzwi
+                     zws(ji,jj,jk) = zzws
+                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
                   END DO
                END DO
@@ -255,9 +257,9 @@
                   DO ji = fs_2, fs_jpim1   ! vector opt.
                      ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
-                     zzwi = ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_a(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
-                     zzws = ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
-                     zwi(ji,jj,jk) = - zdt * zzwi
-                     zws(ji,jj,jk) = - zdt * zzws
-                     zwd(ji,jj,jk) = 1._wp + zzwi + zzws
+                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
+                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+                     zwi(ji,jj,jk) = zzwi
+                     zws(ji,jj,jk) = zzws
+                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
                   END DO
                END DO
@@ -361,14 +363,15 @@
                   DO jj = 2, jpjm1 
                      DO ji = fs_2, fs_jpim1   ! vector opt.
-                        z1_e3vn = 1._wp / e3v_n(ji,jj,jk)
-                        zzwi = ( ( avm   (ji,jj+1,jk  ) + avm   (ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
-                           &     / e3vw_a(ji,jj  ,jk  ) ) * z1_e3vn * wvmask(ji,jj,jk  )
-                        zzws = ( ( avm   (ji,jj+1,jk+1) + avm   (ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
-                           &     / e3vw_a(ji,jj  ,jk+1) ) * z1_e3vn * wvmask(ji,jj,jk+1)
+                        ze3va   = e3v_a(ji,jj,jk)
+                        z1_e3va = 1._wp / ze3va
+                        zzwi = - zdt * ( avm   (ji,jj+1,jk  ) + avm   (ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
+                           &             / ( ze3va * e3vw_n(ji,jj  ,jk  ) ) * wvmask(ji,jj,jk  )
+                        zzws = - zdt * ( avm   (ji,jj+1,jk+1) + avm   (ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
+                           &             / ( ze3va * e3vw_n(ji,jj  ,jk+1) ) * wvmask(ji,jj,jk+1)
                         zWv  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji,jj+1,jk  )   &
                            &               + wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1)   )
-                        zwi(ji,jj,jk) = - 1._wp * ( zzwi + MIN( zWv, 0._wp ) * z1_e3vn )
-                        zws(ji,jj,jk) = - 1._wp * ( zzws - MAX( zWv, 0._wp ) * z1_e3vn )
-                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWv, 0._wp ) + MIN( zWv, 0._wp ) ) * z1_e3vn )
+                        zwi(ji,jj,jk) = zzwi - zdt * MIN( zWv, 0._wp ) * z1_e3va
+                        zws(ji,jj,jk) = zzws + zdt * MAX( zWv, 0._wp ) * z1_e3va
+                        zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWv, 0._wp ) + MIN( zWv, 0._wp ) ) * z1_e3va
                      END DO
                   END DO
@@ -378,14 +381,15 @@
                   DO jj = 2, jpjm1 
                      DO ji = fs_2, fs_jpim1   ! vector opt.
-                        z1_e3vn = 1._wp / e3v_n(ji,jj,jk)
-                        zzwi = ( ( avm   (ji,jj+1,jk  ) + avm(ji,jj,jk  ) )   &
-                           &     / e3vw_a(ji,jj  ,jk  ) ) * z1_e3vn * wvmask(ji,jj,jk  )
-                        zzws = ( ( avm   (ji,jj+1,jk+1) + avm(ji,jj,jk+1) )   &
-                           &     / e3vw_a(ji  ,jj,jk+1) ) * z1_e3vn * wvmask(ji,jj,jk+1)
+                        ze3va   = e3v_a(ji,jj,jk)
+                        z1_e3va = 1._wp / ze3va
+                        zzwi = - zdt * ( avm   (ji,jj+1,jk  ) + avm(ji,jj,jk  ) )   &
+                           &             / ( ze3va * e3vw_n(ji,jj  ,jk  ) ) * wvmask(ji,jj,jk  )
+                        zzws = - zdt * ( avm   (ji,jj+1,jk+1) + avm(ji,jj,jk+1) )   &
+                           &             / ( ze3va * e3vw_n(ji  ,jj,jk+1) ) * wvmask(ji,jj,jk+1)
                         zWv  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji,jj+1,jk  )   &
                            &               + wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1)   )
-                        zwi(ji,jj,jk) = - 1._wp * ( zzwi + MIN( zWv, 0._wp ) * z1_e3vn )
-                        zws(ji,jj,jk) = - 1._wp * ( zzws - MAX( zWv, 0._wp ) * z1_e3vn )
-                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWv, 0._wp ) + MIN( zWv, 0._wp ) ) * z1_e3vn )
+                        zwi(ji,jj,jk) = zzwi - zdt * MIN( zWv, 0._wp ) * z1_e3va
+                        zws(ji,jj,jk) = zzws + zdt * MAX( zWv, 0._wp ) * z1_e3va
+                        zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWv, 0._wp ) + MIN( zWv, 0._wp ) ) * z1_e3va
                      END DO
                   END DO
@@ -399,13 +403,13 @@
                      DO ji = fs_2, fs_jpim1   ! vector opt.
                         ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at U-point
-                        zzwi = ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
-                           &         / ( ze3va * e3vw_a(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
-                        zzws = ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
-                           &         / ( ze3va * e3vw_a(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
-                        zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) )
-                        zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) )
-                        zwi(ji,jj,jk) = - 1._wp * ( zzwi +  MIN( zWvi, 0._wp ) ) 
-                        zws(ji,jj,jk) = - 1._wp * ( zzws -  MAX( zWvs, 0._wp ) )
-                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
+                        zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
+                           &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
+                        zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
+                           &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                        zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  )
+                        zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1)
+                        zwi(ji,jj,jk) = zzwi - zdt * MIN( zWvi, 0._wp )
+                        zws(ji,jj,jk) = zzws + zdt * MAX( zWvs, 0._wp )
+                        zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
                      END DO
                   END DO
@@ -416,11 +420,11 @@
                      DO ji = fs_2, fs_jpim1   ! vector opt.
                         ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at U-point
-                        zzwi = ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_a(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
-                        zzws = ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_a(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
-                        zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) )
-                        zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) )
-                        zwi(ji,jj,jk) = - 1._wp * ( zzwi +  MIN( zWvi, 0._wp ) ) 
-                        zws(ji,jj,jk) = - 1._wp * ( zzws -  MAX( zWvs, 0._wp ) )
-                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
+                        zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
+                        zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                        zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  )
+                        zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1)
+                        zwi(ji,jj,jk) = zzwi  - zdt * MIN( zWvi, 0._wp )
+                        zws(ji,jj,jk) = zzws  + zdt * MAX( zWvs, 0._wp )
+                        zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
                      END DO
                   END DO
@@ -439,6 +443,6 @@
                      zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
                         &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
-                     zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
-                     zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
+                     zwi(ji,jj,jk) = zzwi
+                     zws(ji,jj,jk) = zzws
                      zwd(ji,jj,jk) = 1._wp - zzwi - zzws
                   END DO
@@ -452,6 +456,6 @@
                      zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
                      zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
-                     zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
-                     zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
+                     zwi(ji,jj,jk) = zzwi
+                     zws(ji,jj,jk) = zzws
                      zwd(ji,jj,jk) = 1._wp - zzwi - zzws
                   END DO
Index: /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/TRA/trazdf.F90
===================================================================
--- /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/TRA/trazdf.F90	(revision 10007)
+++ /NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/TRA/trazdf.F90	(revision 10008)
@@ -136,5 +136,5 @@
       !
       INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp) ::  zrhs             ! local scalars
+      REAL(wp) ::  zrhs, zzwi, zzws ! local scalars
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwi, zwt, zwd, zws
       !!---------------------------------------------------------------------
@@ -180,9 +180,11 @@
                DO jk = 1, jpkm1
                   DO jj = 2, jpjm1
-                     DO ji = fs_2, fs_jpim1   ! vector opt.
-                        zwi(ji,jj,jk) = - p2dt * ( zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  )  + MIN( wi(ji,jj,jk  ) , 0._wp ) )
-                        zws(ji,jj,jk) = - p2dt * ( zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1)  - MAX( wi(ji,jj,jk+1) , 0._wp ) )
-                        zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk)   &
-                           &                   - p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )
+                     DO ji = fs_2, fs_jpim1   ! vector opt. (ensure same order of calculation as below if wi=0.)
+                        zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  )
+                        zzws = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1)
+                        zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zzwi - zzws   &
+                           &                 - p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )
+                        zwi(ji,jj,jk) = zzwi - p2dt *   MIN( wi(ji,jj,jk  ) , 0._wp )
+                        zws(ji,jj,jk) = zzws + p2dt *   MAX( wi(ji,jj,jk+1) , 0._wp )
                     END DO
                   END DO
