anchor.diff:  
     integer nx, ny, layers
     integer i, j, k
  -  double precision h(0:nx, 0:ny, layers)
     double precision u(0:nx+1, 0:ny+1, layers)
     double precision v(0:nx+1, 0:ny+1, layers)
     double precision zeta(0:nx+1, 0:ny+1, layers)
  
positive.code: subroutine calc_eta_star(ub, vb, eta, etastar, freesurfFac, nx, ny, dx, dy, dt)
    implicit none
  
    integer nx, ny, layers
    integer i, j, k
    double precision eta(0:nx+1, 0:ny+1)
    double precision etastar(0:nx+1, 0:ny+1)
    double precision ub(0:nx+1, 0:ny+1)
  
negative.code:   double precision :: H0 ! default depth in no file specified
    ! Pressure solver variables
    double precision :: a(5, nx, ny)
  
    ! Bernoulli potential and relative vorticity
    double precision :: b(0:nx+1, 0:ny+1, layers)
  