anchor.diff:  p = p*0.0
   call initialize_actuator()
   
  -do it=1,nts
  +time_step: do it=1,nts
   
       call apply_boundary(u,v,1,it)
  +
       !! Calculate the uvx and uvy with upwinding
       gamma = min(1.2*dt*max(maxval(abs(u))/hx,maxval(abs(v))/hy),1.0);
       ua = u; call yavg(ua); ua = cshift(ua,shift=-1,dim=2);
  
positive.code:     u = u - dt*(uvy+u2x)
      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
  
  
  
  
negative.code: program cfd
  ! Run the 2d cfd simulation
  
  use parameters
  use diff
  