anchor.diff:      end do
   
   ! wrap fields around for periodic simulations
  -    wrap_fields_3D(unew,nx,ny,layers)
  -    wrap_fields_3D(vnew,nx,ny,layers)
  -    wrap_fields_3D(hnew,nx,ny,layers)
  -    wrap_fields_2D(etanew,nx,ny)
  +    call wrap_fields_3D(unew,nx,ny,layers)
  +    call wrap_fields_3D(vnew,nx,ny,layers)
  +    call wrap_fields_3D(hnew,nx,ny,layers)
  +    call wrap_fields_2D(etanew,nx,ny)
   
   
   !    Accumulate average fields
  
positive.code:     integer nx,ny, i,j, maxits, n
      double precision dt
      double precision rjac, eps
      double precision rhs(nx,ny), 
      double precision res(nx,ny)
      double precision norm, norm0, norm_old
      double precision relax_param
  
negative.code:     double precision depth(0:nx+1,0:ny+1)
      double precision H0 ! default depth in no file specified
  !   Pressure solver variables
      double precision a(5,0:nx+1,0:ny+1)
      double precision phi(0:nx+1,0:ny+1), phiold(0:nx+1,0:ny+1)
  !
  !     Bernoulli potential and relative vorticity 
  