anchor.diff:  
     ! Free surface (eta)
     double precision, dimension(:,:), allocatable :: eta
  -  double precision, dimension(:,:), allocatable :: etastar
     double precision, dimension(:,:), allocatable :: etanew
     ! for saving average fields
     double precision, dimension(:,:), allocatable :: etaav
  
positive.code:   allocate(dudtold(0:nx+1, 0:ny+1, layers))
    allocate(dudtveryold(0:nx+1, 0:ny+1, layers))
    allocate(unew(0:nx+1, 0:ny+1, layers))
    allocate(ub(nx+1, ny))
    allocate(uhalf(0:nx+1, 0:ny+1, layers))
    allocate(uav(0:nx+1, 0:ny+1, layers))
  
  
negative.code:   double precision, intent(in)    :: g_vec(layers)
    integer,          intent(in)    :: nx, ny, layers, n
  
    ! barotropic velocity components (for pressure solver)
    double precision :: ub(nx+1, ny)
    double precision :: vb(nx, ny+1)
    double precision :: etastar(0:nx+1, 0:ny+1)
  