c***************************************************************************
c
c (c) Joan Masso: 28 April 1994: quick and dirty 2d transport equation
c Extended on 19 Apr 1995 to use fortran 90 instead of f77
c
c It evolves the equation:
c u,t + u,x + u,y = 0
c Using a maccormack scheme.
c The initial data is a cruddy gaussian.
c Boundaries are flat: copying the value of the neighbour
c
c***************************************************************************
c Initializes the field with a cruddy gaussian
subroutine initial(
$ u,ulb,uub,ushape,
$ time,dx,dy)
implicit none
integer ulb(2), uub(2), ushape(2)
real*8 u(ushape(1),ushape(2))
real*8 dx,dy
real*8 time
real*8 roffset
real*8 sigma
parameter (roffset=0.8)
parameter (sigma=.1)
integer nx,ny
integer sx,sy
integer i,j
real*8 x,y
nx = ushape(1)
ny = ushape(2)
sx = (uub(1)-ulb(1))/(nx-1)
sy = (uub(2)-ulb(2))/(ny-1)
do j = 1, ny
y = (ulb(2)+(j-1)*sy)*dy
do i=1,nx
x = (ulb(1)+(i-1)*sx)*dx
u(i,j) = exp(-((sqrt(x**2+y**2)-roffset-time)
$ /sigma)**2)
enddo
enddo
return
end
c
c
c ****************************************************************
c
c evolve the field.
subroutine evolveP(
* u,ulb,uub,ushape,
* upp,uplb,upub,upshape,
* dt,dx,dy);
implicit none
integer ulb(2), uub(2), ushape(2)
real*8 u(ushape(1),ushape(2))
integer uplb(2), upub(2), upshape(2)
real*8 upp(upshape(1),upshape(2))
real*8 dt,dx,dy
integer i,j
integer nx,ny
nx = ushape(1)
ny = ushape(2)
c Predictor step: backward derivatives. Note that we get a predicted
c value at the exterior.
u(2:,2:) = upp(2:,2:)
$ - dt/dx* ( upp(2:,2:) - upp(1:nx-1,2:) )
$ - dt/dy* ( upp(2:,2:) - upp(2:,1:ny-1) )
return
end
c ************************************************************
subroutine evolveC(
$ u,ulb,uub,ushape,
$ upp,uplb,upub,upshape,
$ dt,dx,dy);
implicit none
integer ulb(2), uub(2), ushape(2)
real*8 u(ushape(1),ushape(2))
integer uplb(2), upub(2), upshape(2)
real*8 upp(upshape(1),upshape(2))
real*8 dt,dx,dy
integer i,j
integer nx,ny
nx = ushape(1)
ny = ushape(2)
c Corrector step: forward derivatives. Average corrector with old value.
u(2:nx-1,2:ny-1) = (
$ u(2:nx-1,2:ny-1) + upp(2:nx-1,2:ny-1)
$ - dt/dx*(u(3:nx,2:ny-1)-u(2:nx-1,2:ny-1))
$ - dt/dy*(u(2:nx-1,3:ny)-u(2:nx-1,2:ny-1)))/2.0
return
end