c
c =====================================================
subroutine ic(maxmx,maxmy,meqn,mbc,mx,my,x,y,dx,dy,q)
c =====================================================
c
c Copyright (C) 2003-2007 California Institute of Technology
c Ralf Deiterding, ralf@amroc.net
c
implicit none
c
include "cuser.i"
c
integer i,j,meqn,maxmx,maxmy,my,mx,mbc,n,m,steps
double precision r,ur,p,phi,dx,dy,y,x,
& q,xp,yp,qs,dxs,dys
dimension q(meqn, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc)
dimension x(1-mbc:maxmx+mbc),y(1-mbc:maxmy+mbc)
dimension qs(4)
c
steps = 4
do 160 i = 1-mbc, mx+mbc
do 160 j = 1-mbc, my+mbc
q(1,i,j) = 0.d0
q(2,i,j) = 0.d0
q(3,i,j) = 0.d0
q(4,i,j) = 0.d0
dxs = dx/(2*steps+1)
dys = dy/(2*steps+1)
do 170 n = -steps, steps
do 170 m = -steps, steps
xp = x(i)+dxs*n
yp = y(j)+dys*m
r = ( (xp-xlc)**2 + (yp-ylc)**2 )**0.5d0
phi = datan( (yp-ylc) / (xp-xlc) )
if( xp-xlc .lt. 0.d0 ) phi=phi+pi
if (r .lt. rvortex) then
if(r .lt. rvortex/2.d0) then
ur = ua*2.d0*r/rvortex
p = (r/rvortex)**2+1.d0-2.d0*dlog(2.d0)
else
ur = ua*2.d0*(1.d0-r/rvortex)
p = (r/rvortex)**2+3.d0-4.d0*r/rvortex+
& 2.d0*dlog(r/rvortex)
end if
qs(2) = -sin(phi)*ur
qs(3) = cos(phi)*ur
else
p=0.d0
qs(2)=0.d0
qs(3)=0.d0
endif
qs(1) = rho0
qs(2) = qs(1)*qs(2)
qs(3) = qs(1)*qs(3)
p = p0 + 2.d0*rho0*ua**2*p
qs(4) = p/gamma1+0.5d0*(qs(2)**2+qs(3)**2)/qs(1)
q(1,i,j) = q(1,i,j) + qs(1)*dxs*dys
q(2,i,j) = q(2,i,j) + qs(2)*dxs*dys
q(3,i,j) = q(3,i,j) + qs(3)*dxs*dys
q(4,i,j) = q(4,i,j) + qs(4)*dxs*dys
170 continue
q(1,i,j) = q(1,i,j) / (dx*dy)
q(2,i,j) = q(2,i,j) / (dx*dy)
q(3,i,j) = q(3,i,j) / (dx*dy)
q(4,i,j) = q(4,i,j) / (dx*dy)
160 continue
return
end
c
c