anchor.diff:         call write_netcdf(grid(lev)%cA,vname='ca',netcdf_file_name='cA.nc',rank=myrank,iter=lev)
       endif
   
  -  end subroutine define_matrix
  +  end subroutine define_matrix_orig
   
   end module mg_define_matrix
  
positive.code:           cA(3,k,j,i) = qrt * ( &
                 ( hlf * (zr(k+1,j+1,i  )-zr(k+1,j-1,i  )) / dy(j,i) ) * dx(j,i) + &
                 ( hlf * (zr(k,j,i  )-zr(k,j-2,i  )) / dy(j-1,i) ) * dx(j-1,i) )
                 !! zydx(k+1,j,i) + &
                 !! zydx(k,j-1,i) ) ! couples with k+1 j-1
            cA(4,k,j,i) =                                                                &
                                  ! couples with j-1
                 ( qrt                                                                 * &
  
negative.code:                ( zw(k+1,j,i) - zw(k,j,i) + zw(k+1,j-1,i) - zw(k,j-1,i) ) * &
                 (dx(j,i)+dx(j-1,i)) ) / ( hlf * (dy(j,i)+dy(j-1,i)) ) & ! couples with j-1
                 + qrt * ( &
                 - (( hlf * (zr(k,j  ,i)-zr(k,j-2,i)) / dy(j-1,i) ) * dx(j-1,i)) &
                 + (( hlf * (zr(k,j+1,i)-zr(k,j-1,i)) / dy(j  ,i) ) * dx(j  ,i)) )
            cA(5,k,j,i) =- qrt * ( &
                 (( hlf * (zr(k-1,j+1,i  )-zr(k-1,j-1,i  )) / dy(j,i) ) * dx(j,i)) + &
                 (( hlf * (zr(k,j,i  )-zr(k,j-2,i  )) / dy(j-1,i) ) * dx(j-1,i)) )     ! couples with k-1 j-1
  