anchor.diff:  program cfd
  +! Run the 2d cfd simulation
   
   use parameters
   use diff
  
positive.code: p = p*0.0
  call initialize_actuator()
  
  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);
  
negative.code:     u = u - dt*(uvy+u2x)
      v = v - dt*(uvx+v2y)
  
      !! Iteratively solve for pressure
      call apply_boundary(u,v,1,it)
  
      ud = u; call ddx(ud,hx)
      vd = v; call ddy(vd,hy)
      rhs = ud + vd
      iterate_pressure: do ii=1,niter
          p = 0.25*(myiter(p,hm) - hx*hy*rhs)
      end do iterate_pressure
      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
  
      !! Apply actuator force
      call update_actuator(u)
      call apply_actuator(u,dt,hx,pow)
  
      !! driving pressure gradient
      u = u + dt*dpg
  
      !! Save the information if requested
      if (mod(it,write_freq) == 0) then
          call save_vel(it,u,v)
          call save_pow(it,pow)
      end if 
  
  end do time_step
  
  
  
  