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:     ! h2av = 0.0
  
    end if
  
  
  
    return
  end subroutine maybe_dump_output
  