anchor.diff:    call evaluate_zeta(zeta, u, v, nx, ny, layers, dx, dy)
     if (debug_level .ge. 4) then
       call write_output_3d(zeta, nx, ny, layers, 1, 1, &
  -      n, 'snap.zeta.')
  +      n, 'output/snap.zeta.')
     end if
   
     ! Calculate dhdt, dudt, dvdt at current time step
  
positive.code:     call evaluate_b_RedGrav(b, h, u, v, nx, ny, layers, g_vec)
      if (debug_level .ge. 4) then
        call write_output_3d(b, nx, ny, layers, 0, 0, &
          n, 'snap.BP.')
      end if
    else
      call evaluate_b_iso(b, h, u, v, nx, ny, layers, g_vec, depth)
      if (debug_level .ge. 4) then
        call write_output_3d(b, nx, ny, layers, 0, 0, &
          n, 'snap.BP.')
      end if
    end if
  
  
negative.code:   ! print *, maxval(abs(etastar))
    if (debug_level .ge. 4) then
      call write_output_2d(etastar, nx, ny, 0, 0, &
        n, 'output/snap.eta_star.')
    end if
  
    ! Prevent barotropic signals from bouncing around outside the
  