	program sw
c
c	triple deck solver 
c	modified for Stewartson  & Williams free interaction
c       in a compressible flow
c
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      parameter (iy=50,iyp=iy+1,iz=3*(iy+1)+1,ix=120)


	real *8 u(0:ix,0:iy),du(0:iyp)
	real *8 w(0:ix,0:iy),dw(0:iyp)
	real *8 psi(0:ix,0:iy),dpsi(0:iyp)
	real *8 v(0:ix,0:iy)
	real *8 p(0:ix),px(0:ix),dpx
	real *8 a(0:ix),ad(0:ix)
	real *8 a0(0:ix)
	real *8 cof(iz,iz),rhs(iz)
	real *8 x(0:ix),hx(0:ix),hy(0:iy),y(0:iy)
	real *8 dxi,xm,z,zm,third,hxi
	real *8 hs,h2,alpha,fac,fac1,xt,rho,nu,var
	real *8 tolerance,hx0
	integer ipx,ipy,ix,iy,i,j,k,msize,msizep,row,nits

	nits=5

	tolerance=1.e-6
	ipx=1
	ipy=1
	iterations=10
c
c	set up y mesh
c
	hy0=.05
	facy=1.07
	y(0)=0
	hy(1)=hy0
	y(1)=y(0)+hy0
	do j=2,iy
		hy(j)=facy*hy(j-1)
		y(j)=y(j-1)+hy(j)
	enddo
	zm=y(iy)
c
c	set up x mesh
c
	hx0=0.1
	do i=0,ix
	   hx(i)=hx0
	   x(i)=i*hx0
	enddo
	xm=x(ix)

	write(7,7)(x(i),hx(i),i=0,ix,ipx)
	write(8,7)(y(j),hy(j),j=0,iy,ipy)
7	format(2e15.5)
	
	write(6,5)ix,iy,zm,xm,hx0,hy0
5	format('Interactive bl solver',/, 'ix = ',i4,
     +   ' iy  = ',i4,
     +   /,'zm =',f8.0, 'xm = ', f12.4,/,
     +     '  hx0 = ',f10.4, ' hy0 = ',f10.4)
c
c	set up initial guess for a(x),p(x)
c
	write(6,1)
1	format('Initialising a(x), p(x)')
	do i=0,ix
		ad(i)=0.
		a(i)=0.
		p(i)=0.
		px(i)=0.
	end do
c
c	introduce initial perturbation
c
	p(0)=0.004
	
	write(6,3)
3	format('Setting up plate start solution')
	do i=0,ix
		u(i,0)=0.
		psi(i,0)=0.
		w(i,0)=1.
		do j=1,iy
			z=y(j)
			psi(i,j)=(z*z/2+z*a(i))
			u(i,j)=(z+a(i))
			w(i,j)=1.
		end do
	end do
c
c	now have initial guess for all variables,
c	output intial values
c
67	format('vdata',/,(f12.3,f15.6))
66	format('udata',/,(f12.3,f15.6))
65	format('data',/, (f12.3,f15.6))
c--------------------------------------------------------------------
c	main iteration                                              |
c--------------------------------------------------------------------
	do k=1,iterations
	   do i=0,ix
		u(i,iy)=(zm+a(i))
	   enddo
	   do i=1,ix
c--------------------------------------------------------------------
c       iterate along plate rows - u(0)=0                           |
c--------------------------------------------------------------------
	do newton=1,nits
		do l=1,iz
		do m=1,iz
			cof(l,m)=0.
		enddo
		enddo
	row=1
c	set first two rows
		cof(row,row)=1
		rhs(row)=0
		row=row+1
		cof(row,row)=1
		rhs(row)=0
c	loop through iy sets
		do j=1,iy
		        wc=(w(i,j)+w(i-1,j)+
     +                      w(i,j-1)+w(i-1,j-1))/4
			uc=(u(i,j)+u(i-1,j)+
     +                      u(i,j-1)+u(i-1,j-1))/4
			phiy=psi(i,j)+  psi(i-1,j)-
     +                       psi(i,j-1)-psi(i-1,j-1)
			uy=    u(i,j)+    u(i-1,j)-
     +                         u(i,j-1)-  u(i-1,j-1)
			wy=    w(i,j)+    w(i-1,j)-
     +                         w(i,j-1)-  w(i-1,j-1)
		        phix=(psi(i,  j)+psi(i,  j-1)-
     +                        psi(i-1,j)-psi(i-1,j-1))/(2*hx(i))
		        ux=(    u(i,  j)+  u(i,  j-1)-
     +                          u(i-1,j)-  u(i-1,j-1))/(2*hx(i))

			row=row+1
			cof(row,row-2)=-1
			cof(row,row-1)=-2*hy(j)/4
			cof(row,row+1)=1
			cof(row,row+2)=-2*hy(j)/4
			rhs(row)=2*hy(j)*uc-phiy

			row=row+1
			cof(row,row-2)=-1
			cof(row,row-1)=-2*hy(j)/4
			cof(row,row+1)=1
			cof(row,row+2)=-2*hy(j)/4
			rhs(row)=2*hy(j)*wc-uy

			row=row+1
			cof(row,row-4)=   2*hy(j)*wc/(2*hx(i))
			cof(row,row-3)=  -2*hy(j)*ux/4-2*hy(j)*uc/(2*hx(i))
			cof(row,row-2)=-1+2*hy(j)*phix/4
			cof(row,row-1)=   2*hy(j)*wc/(2*hx(i))
			cof(row,row)  =  -2*hy(j)*ux/4-2*hy(j)*uc/(2*hx(i))
			cof(row,row+1)= 1+2*hy(j)*phix/4
			cof(row,iz)   =  -2*hy(j)
			rhs(row)=2*hy(j)*(px(i)
     +					+ uc*ux - wc*phix) -wy
		enddo

		row=row+1
		cof(row,row-1)=1.
		rhs(row)=0

		row=row+1
		cof(row,row-1)=1.
		rhs(row)=0
c
c	cof set up
c
	call invert(cof,iz,iz,rhs)
	error=0
	do j=0,iy		
		dpsi(j)=rhs(3*j+1)
		psi(i,j)=psi(i,j)+dpsi(j)
		du(j)=rhs(3*j+2)
		u(i,j)=u(i,j)+du(j)
		dw(j)=rhs(3*j+3)
		w(i,j)=w(i,j)+dw(j)
		error=error+abs(dpsi(j))/psi(ixd,iy)
     +                     +abs(du(j))/u(ixd,iy)
     +                     +abs(dw(j))/w(ixd,iy)

	enddo
c	write(6,72)error/(3*iy)
72	format('plate',e12.3)
	dpx=rhs(iz)
	px(i)=px(i)+dpx

	enddo
        enddo
c
c	get pressures from gradients, gradients at midpoints
c
	
	do i=1,ix
		p(i)=p(i-1)+hx(i)*px(i)
		a(i)=a(i-1)-0.5*hx(i)*(p(i)+p(i-1))
	enddo
c
c   calculate a(x) from a'(x)
c


        write(6,10)k,a(0),a(ix/4),a(ix/2)
10      format(i6,(10e12.3))

	enddo

	write(22,21)((psi(i,j),j=0,iy,ipy),i=0,ix,ipx)
	write(23,21)((u(i,j),j=0,iy,ipy),i=0,ix,ipx)
	write(24,21)((w(i,j),j=0,iy,ipy),i=0,ix,ipx)
	write(25,65)(x(i),      a(i),  i=0,ix)
	write(31,65)(x(i),      w(i,0),  i=0,ix)
	write(26,65)(x(i),      p(i),  i=0,ix)
	write(27,65)(x(i),      ad(i), i=0,ix)
	write(28,65)(x(i)-0.5*hx(i),px(i),i=0,ix)
	

21	format(8e15.4)
	stop
	end


  	subroutine invert(a,n,np,b)
	implicit real *8(a-h,o-z)
	real*8 a(np,np),b(np)
c---------------------------------------------------------------------
c
c	YOU MUST PROVIDE A ROUTINE FOR INVERSION HERE
c
c	for instace from Numerical Recipes
c
c---------------------------------------------------------------------
c
c	routine to solve ax=b
c
c	maximum dimension of a is np x np
c
c	part being used is n x n
c
c	solution returned in b
c


	return
	end

