anchor.diff:      vd = v; call ddy(vd,2.0);
   
       u2x = ua*ua - gamma*abs(ua)*ud
  -    call ddx(u2x,h); u2x = cshift(u2x,shift=-1,dim=1)
  +    call ddx(u2x,hx); u2x = cshift(u2x,shift=-1,dim=1)
   
       v2y = va*va - gamma*abs(va)*vd
  -    call ddy(v2y,h); v2y = cshift(v2y,shift=-1,dim=2)
  +    call ddy(v2y,hy); v2y = cshift(v2y,shift=-1,dim=2)
   
   
       !! update the u and v variables
  
positive.code:     v = v - dt*(uvx+v2y)
  
      !! Pressure
  
      ud = u; call ddx(ud,h)
      vd = v; call ddy(vd,h)
      rhs = ud + vd
      do ii=1,niter
          p = 0.25*(myiter(p,h) - h**2*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)
      u = u - ud
      v = v - vd
  
      call update_actuator(u)
      call apply_actuator(u,dt,h)
  
      if (mod(it,10) == 0) then
          call save_vel(it,u,v)
      end if 
  
  end do
  
negative.code: use diff
  use save
  use actuator
  use boundary
  implicit none
  
  integer,parameter :: seed = 86456
  integer :: ix, iy, it, ii
  
  
  real :: gamma, pow
  
  real, allocatable,dimension(:,:) :: u,v,p,ua,va,ud,vd,uvx,uvy,u2x,v2y,rhs
  
  