anchor.diff:  use diff
   use save
   use actuator
  +use boundary
   implicit none
   
   integer,parameter :: seed = 86456
  -real,parameter :: pi = 3.14159
   integer :: ix, iy, it, ii
   
   
  -real :: gamma
  +real :: gamma, pow
   
   real, allocatable,dimension(:,:) :: u,v,p,ua,va,ud,vd,uvx,uvy,u2x,v2y,rhs
   
  
positive.code: call ddx(v,1.0)
  u = u + uvel
  p = p*0.0
  call initialize_actuator(10)
  
  do it=1,nts
  
      !! Calculate the uvx and uvy with upwinding
      gamma = min(1.2*dt*max(maxval(abs(u))/h,maxval(abs(v))/h),1.0);
      ua = u; call yavg(ua); ua = cshift(ua,shift=-1,dim=2);
      ud = u; call ddy(ud,2.0); ud = cshift(ud,shift=-1,dim=2);
      va = v; call xavg(va); va = cshift(va,shift=-1,dim=1);
      vd = v; call ddx(vd,2.0); vd = cshift(vd,shift=-1,dim=1);
  
      uvx = (ua * va) - (gamma * abs(ua) * vd)
      call ddx(uvx,h);
  
      uvy = (ua * va) - (gamma * abs(va) * ud)
      call ddy(uvy,h);
  
      !! Calculate the u2x and v2y terms
      ua = u; call xavg(ua);
  
negative.code:     v = v - dt*(uvx+v2y)
  
      !! Pressure
      call apply_boundary(u,v,1,it)
  
      ud = u; call ddx(ud,hx)
      vd = v; call ddy(vd,hy)
      rhs = ud + vd
      do ii=1,niter
          p = 0.25*(myiter(p,hm) - hx*hy*rhs)
      end do
      ud = cshift(p,shift=-1,dim=1); call ddx(ud,hx)
      vd = cshift(p,shift=-1,dim=2); call ddy(vd,hy)
      u = u - ud
      v = v - vd
  
      !! actuator disk
      call update_actuator(u)
      call apply_actuator(u,dt,hx,pow)
  
      !! driving pressure gradient
      u = u + dt*dpg
  
      if (mod(it,10) == 0) then
          call save_vel(it,u,v)
          call save_pow(it,pow)
      end if 
  
  end do
  