anchor.diff:      v = v - dt*(uvx+v2y)
   
       !! Pressure
  +    call apply_boundary(u,v,1,it)
   
  -    ud = u; call ddx(ud,h)
  -    vd = v; call ddy(vd,h)
  +    ud = u; call ddx(ud,hx)
  +    vd = v; call ddy(vd,hy)
       rhs = ud + vd
       do ii=1,niter
  -        p = 0.25*(myiter(p,h) - h**2*rhs)
  +        p = 0.25*(myiter(p,hm) - hx*hy*rhs)
       end do
  -    ud = cshift(p,shift=-1,dim=1); call ddx(ud,h)
  -    vd = cshift(p,shift=-1,dim=2); call ddy(vd,h)
  +    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,h)
  +    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
  
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:     vd = v; call ddy(vd,2.0);
  
      u2x = ua*ua - gamma*abs(ua)*ud
      call ddx(u2x,hx); u2x = cshift(u2x,shift=-1,dim=1)
  
      v2y = va*va - gamma*abs(va)*vd
      call ddy(v2y,hy); v2y = cshift(v2y,shift=-1,dim=2)
  
  
      !! update the u and v variables
  