anchor.diff:      u = u - dt*(uvy+u2x)
       v = v - dt*(uvx+v2y)
   
  -    !! Pressure
  +    !! 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
  -    do ii=1,niter
  +    iterate_pressure: do ii=1,niter
           p = 0.25*(myiter(p,hm) - hx*hy*rhs)
  -    end do
  +    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
   
  -    !! actuator disk
  +    !! Apply actuator force
       call update_actuator(u)
       call apply_actuator(u,dt,hx,pow)
   
       !! driving pressure gradient
       u = u + dt*dpg
   
  -    if (mod(it,10) == 0) then
  +    !! 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
  +end do time_step
   
   
   
  
positive.code: integer :: ix, iy, it, ii
  
  
  real :: gamma, pow
  
  real, allocatable,dimension(:,:) :: u,v,p,ua,va,ud,vd,uvx,uvy,u2x,v2y,rhs
  
  call save_parameters()
  
negative.code: p = p*0.0
  call initialize_actuator()
  
  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);
  