anchor.diff:                    + (eta(i,j)   + eta(i,j+2) - 2*eta(i,j+1)) /(dy*dy)
               temp2 = ((  eta(i+1,j) - eta(i-1,j)) / (2*dx))**2               &
                     + ((4*eta(i,j+1) - eta(i,j+2) - 3*eta(i,j)) / (2*dy))**2
  +
               b(i,j) = (Bu+Ro*eta(i,j))*temp1 + Ro/2*temp2
               q(i,j) = temp1
   
               j = pasy
               temp1 = (eta(i-1,j) + eta(i+1,j) - 2*eta(i,j))   /(dx*dx)       &
                     + (eta(i,j)   + eta(i,j-2) - 2*eta(i,j-1)) /(dy*dy)
  -            temp2 = ((eta(i+1,j) - eta(i-1,j))/(2*dx))**2                     &
  +            temp2 = ((eta(i+1,j) - eta(i-1,j))/(2*dx))**2                   &
                     + ((-4*eta(i,j-1) + eta(i,j-2) + 3*eta(i,j)) / (2*dy))**2
  +
               b(i,j) = (Bu+Ro*eta(i,j))*temp1 + Ro/2*temp2
               q(i,j) = temp1
           end do
  
positive.code:                                - Ro*(temp2*temp4 - temp1*temp5)
  
                      ! v (j=pasy)     
  
                      temp3 = (-4*leta(i,j-1) + leta(i,j-2) + 3*leta(i,j))/2/dy 
                      temp5 = (eta(i,j-2) + eta(i,j) - 2*eta(i,j-1))/(dy*dy)
                      v2(i,j) = temp2 - Et*(temp1 - temp3)/dt         &
  
negative.code:                               -Ro*(temp2*temp4 - temp1*temp5)
  
                      ! v (j=1)      
                      temp3 = (4*leta(i,j+1) - leta(i,j+2) - 3*leta(i,j))/2/dy 
                      temp5 = (eta(i,j+2) + eta(i,j) - 2*eta(i,j+1))/(dy*dy)
                      v2(i,j) = temp2-Et*(temp1 - temp3)/dt             &
                               - Ro*(temp2*temp5 - temp1*temp4)
                      
                      ! u (j=pasy)            
                      j = pasy
                      temp1 = (-4*eta(i,j-1) + eta(i,j-2) + 3*eta(i,j))/2/dy
                      temp2 = (eta(i+1,j) - eta(i-1,j))/2/dx
                      temp3 = (leta(i+1,j) - leta(i-1,j))/2/dx
  