	program tdeckn
c
c	Triple deck solver using Newton iteration
c
c Modified 5 March 2021
c
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      parameter (iy=64,iyp=iy+1,iz=3*(iy+1)+1,ixm=200,ixp=200,
     +             ixm2=2*ixm,ixp2=2*ixp)
	common/wakedata/h0(0:2000),dh0(0:2000),ddh0(0:2000),hwake,nwake


	real *8 u(-ixm:ixp,0:iy),du(0:iyp)
	real *8 w(-ixm:ixp,0:iy),dw(0:iyp)
	real *8 psi(-ixm:ixp,0:iy),dpsi(0:iyp)
	real *8 v(-ixm:ixp,0:iy)
	real *8 p(-ixm2:ixp2),px(-ixm2:ixp2),dpx
	real *8 a(-ixm2:ixp2),ad(-ixm2:ixp2)
	real *8 a0(-ixm:ixp)
	real *8 cof(iz,iz),rhs(iz)
	real *8 x(-ixm2:ixp2),hx(-ixm2:ixp2),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 mu0,gamma0,lambda0,zeta0,tolerance,hx0
	integer ixu,ixd,ipx,ipy,iy,i,j,k,msize,msizep,row,nits

c	set up initial guess

	real *8 gamma
	mu0    =  1.132138751
	gamma0 =  0.244547191
	lambda0=  0.87890145
	zeta0=mu0/lambda0
	third=1./3.
	nits=5
	a1=.3321
	xscale=a1**(1.25)
	pscale=a1**(-.5)
	ascale=a1**.75

	tolerance=1.e-6
	ixu=200
	ixd=200
	ipx=5
	ipy=1
	msize=ixm2
	msizep=ixp2

	iterations=30
c
c	set up y mesh
c
	hy0=.05
	facy=1.05
	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.02
	facu=1.05
	facd=1.05
	x(0)=0.
	x(1)=hx0
	hx(1)=hx0
	hx(0)=hx0
	x(-1)=-hx0
	do i=2,msizep
		if(i.le.40)then
		hx(i)=facd*hx(i-1)
		else

		hx(i)=1.01*hx(i-1)
		endif
		x(i)=x(i-1)+hx(i)
	enddo
	do i=2,msize
		if(i.le.40)then
		hx(-i+1)=facu*hx(-i+2)
		else
		hx(-i+1)=1.01*hx(-i+2)
		endif

		x(-i)=x(-i+1)-hx(-i+1)

	enddo
	xm=x(ixd)

	write(7,7)(x(i),hx(i),i=-ixu,ixd,ipx)
	write(8,7)(y(j),hy(j),j=0,iy,ipy)
7	format(2e15.5)

	write(6,5)ixu,ixd,iy,zm,xm,hx0,hy0,msize,msizep
5	format('Interactive bl solver',/, 'ixu = ',i4,
     +   ' ixd = ',i4, ' iy  = ',i4,
     +   /,'zm =',f8.0, 'xm = ', f12.4,/,
     +     '  hx0 = ',f10.4, ' hy0 = ',f10.4,2i10)
        write(6,10)k,var,x(-ixu/2),x(0),x(ixd/4),
     +      x(ixd/2),x(3*ixd/4),x(ixd),p(0)

	gamma=zeta0/(3.*3.**.5)
        alpha=2*lambda0*lambda0*lambda0*gamma0

c
c	set up initial guess for a(x),p(x)
c
	write(6,1)
1	format('Initialising a(x), p(x)')
	do i=-msize,-1
		ad(i)=0.
		a(i)=0.
		xt=(-x(i))**third
		p(i)=-2*gamma/(xt*xt)
		xc=-(x(i)-0.5*hx(i))
		xt=xc**third
		px(i)=-2*2*gamma/(3*xc*xt*xt)
	end do
		ad(0)=0
		a(0)=0
		p(0)=0
		xc=0.5*hx(0)
		xt=xc**third
		px(0)=-2*2*gamma/(3*xc*xt*xt)
	do i=1,msizep
		xt=(x(i))**third
		ad(i)=zeta0/(3*xt*xt)
		a(i)=zeta0*xt
		p(i)= gamma/(xt*xt)
		xc=x(i)-0.5*hx(i)
		xt=xc**third
		px(i)=2*gamma/(3*xc*xt*xt)
	end do
c
c	now have asymptotic form for a,p
c
	write(6,2)
2	format('Calling Hilbert transform')
	call hilb(ad,msize,msizep,p,ixu,ixd,hx,x)
c
c	calculate how numerical hilbert transform sees a(x) from p(x)
c
	a(-ixu)=0
	do i=-ixu+1,ixd
		a(i)=a(i-1)+hx(i)*(ad(i-1)+ad(i))/2
	enddo
	write(6,3)
3	format('Setting up plate start solution')
	do i=-ixu,0
		u(i,0)=0.
		psi(i,0)=0.
		w(i,0)=alpha
		do j=1,iy
			z=y(j)
			psi(i,j)=alpha*(z*z/2+z*a(i))
			u(i,j)=alpha*(z+a(i))
			w(i,j)=alpha
		end do
	end do
	sm=zm/(xm**third)
	write(6,4)
4	format('Setting up wake start solution')
c
c	intialise wake solution
c
	call wake0
	do i=1,ixd
		xt=x(i)**third
		do j=0,iy
			z=y(j)
			s=z/xt
			if(s.le.sm)then
		 	   call wake(s,ho,dho,ddho)
			   psi(i,j)=xt*xt*ho
			   u  (i,j)=xt*dho
			   w(i,j)=ddho
           		else
			   u(i,j)=alpha*(z+a(i))
			   psi(i,j)=alpha*(z*z/2+z*a(i))
			   w(i,j)=alpha
			endif
		end do
	end do
c
c	now have initial guess for all variables,
c	output intial values
c
	write(13,21)((psi(i,j),j=0,iy,ipy),i=-ixu,ixd,ipx)
	write(14,21)((u(i,j),j=0,iy,ipy),i=-ixu,ixd,ipx)
	write(15,65)(xscale*x(i),ascale*a(i),  i=-msize,msizep)
	write(16,65)(xscale*x(i),pscale*p(i),  i=-msize,msizep)
	write(17,65)(x(i),ad(i), i=-msize,msizep)
	write(18,65)(x(i)-0.5*hx(i),px(i),i=-msize+1,msizep)
	write(19,66)(y(j),u(-ixu,j),j=0,iy)
	write(19,66)(y(j),u(0,j),j=0,iy)
	write(19,66)(y(j),u(ixd/4,j),j=0,iy)
	write(19,66)(y(j),u(ixd/2,j),j=0,iy)
	write(19,66)(y(j),u(3*ixd/4,j),j=0,iy)
	write(19,66)(y(j),u(ixd,j),j=0,iy)
	write(20,67)(x(i),w(i,0),i=-ixu,0),(x(i),u(i,0),i=1,ixd)
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=-ixu,ixd
		u(i,iy)=alpha*(zm+a(i))
	   enddo
	   do i=-ixu+1,0
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            iterate along wake                                   |
c------------------------------------------------------------------
	   do i=1,ixd
	do newton=1,10
		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)=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
	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,71)error/(3*iy)
71	format('wake ',e12.5)
	dpx=rhs(iz)
	px(i)=px(i)+dpx
	enddo
	enddo
c
c	get pressures from gradients, gradients at midpoints
c

	do i=-ixu+1,ixd
		p(i)=p(i-1)+hx(i)*px(i)
	enddo
c
c ----------------	use Hilbert transform to get a' from p(x)
c
	   call hilb(ad,msize,msizep,p,ixu,ixd,hx,x)
c
c   calculate a(x) from a'(x)
c

	a(-ixu)=0.

	omega=.05
	a0(-ixu)=0.
	do i=-ixu+1,ixd
		a0(i)=a0(i-1)+hx(i)*(ad(i-1)+ad(i))/2
		a(i)=(1-omega)*a(i)+omega*a0(i)
	enddo


        write(6,10)k,a(-ixu/2),a(0),a(ixd/4),
     +      a(ixd/2),a(3*ixd/4),a(ixd),p(0),px(1)
10      format(i6,(10e12.3))

	enddo

	write(22,21)((psi(i,j),j=0,iy,ipy),i=-ixu,ixd,ipx)
	write(23,21)((u(i,j),j=0,iy,ipy),i=-ixu,ixd,ipx)
	write(24,21)((w(i,j),j=0,iy,ipy),i=-ixu,ixd,ipx)
	write(25,65)(xscale*x(i),      ascale*a(i),  i=-ixu,ixd)
	write(26,65)(xscale*x(i),pscale*p(i),  i=-ixu,ixd)
	write(27,65)(x(i),      ad(i), i=-ixu,ixd)
	write(28,65)(x(i)-0.5*hx(i),px(i),i=-msize+1,msizep)
	write(29,66)(y(j),u(-ixu,j),j=0,iy)
	write(29,66)(y(j),u(0,j),j=0,iy)
	write(29,66)(y(j),u(ixd/4,j),j=0,iy)
	write(29,66)(y(j),u(ixd/2,j),j=0,iy)
	write(29,66)(y(j),u(3*ixd/4,j),j=0,iy)
	write(29,66)(y(j),u(ixd,j),j=0,iy)
	write(30,67)(x(i),w(i,0),i=-ixu,0),(x(i),u(i,0),i=1,ixd)

21	format(8e15.4)
c
c	calculate part of R^(-7/8) drag component in integral in (3.50), page 90
c note that this integral from x=x(-ixu) to x=0 leaves out the long tail
c which decays as (-X)^-(4/3) upstream on the plate which will integrate
c from -infinity to -X=x(ixu) to give an additional part proportional
c to x(-ixu)^(-1/3).
c
c My suggestion is to run the program with two different values
c for x(-ixu) and use the values of the finite length part of the integral,
c together with a factor of size x^(-1/3) to back out overall value for
c integral over -infinity to zero in order to estimate the R^(-7/8) term.
c
c I am grateful to Prof Peter Schaffarczyk for raising this with me
c
	dragintegral=0.
	do i=-ixu,-1
		dragintegral=dragintegral+(w(i+1,0)+w(i,0)-2*alpha)*(x(i+1)-x(i))/2
	enddo
	write(6,25)dragintegral
25	format('Calculated integral from x(-ixu) to 0 for second drag term', f15.6)
	stop
	end
c
c-------------------------------------------------------------------------------hilb-------------
c
	subroutine hilb(ad,msize,msizep,p,ixu,ixd,hx,x)
	real *8 p(-msize:msizep),ad(-msize:msizep)
	real *8 x(-msize:msizep),hx(-msize:msizep)
	integer msize,msizep,ixu,ixd
	real *8 mu0,gamma0,lambda0,zeta0
	real *8 phx,k1,sum,b,t1,t2,dx
	mu0    =  1.132138751
	gamma0 =  0.244547191
	lambda0=  0.87890145
	zeta0=mu0/lambda0
	phx=1./3.14159
	do i=-ixu,ixd
c				pressure at grid points
		sum1=0.
		sum=hx(i+1)*p(i+1)/(x(i)-x(i+1))
		do j=i+1,2*ixd-1
			sum1=sum1+0.5*(hx(j)+hx(j+1))*p(j)/(x(i)-x(j))
		enddo
		sum1=sum1+0.5*p(2*ixd)*hx(2*ixd)/(x(i)-x(2*ixd))
		sum1=sum1*phx
		b=x(2*ixd)
		t2=tail(x(i),b,100)*zeta0/(3**1.5)
		j=-2*ixu
		sum2=0.5*hx(j+1)*p(j)/(x(i)-x(j))
		do j=-2*ixu+1,i-1
			sum2=sum2+0.5*(hx(j)+hx(j+1))*p(j)/(x(i)-x(j))
		enddo
		sum2=sum2*phx
		b=-x(-2*ixu)
		t1=tail(-x(i),b,100)*2*zeta0/(3**1.5)
		ad(i)=-(sum1+sum2+t1+t2)

	enddo
	return
	end

 	real*8 function tail(x,b,m)
	integer i,k,m
	real *8 b1,bx,x,b
	b1=3.*b**(-2./3.)/3.14159
	k=3*m+2
	bx=x/b
	tail=1./k
	do i=1,m
		k=k-3
		tail=1./k+bx*tail
	enddo
	tail=-b1*tail
	return
	end
	subroutine wake(s,ho,dho,ddho)
	implicit real*8(a-h,o-z)
	implicit integer (i-n)
	common/wakedata/h0(0:2000),dh0(0:2000),ddh0(0:2000),h,n
	real *8 a0,b0,df,e
	real *8 fa,fb,fc,fd,dfa,dfb,dfc,dfd,ddfa,ddfb,ddfc,ddfd
	real *8 sp,sm,mu0,mu1,gamma0,gamma1,lambda0
	integer i,j,k,l,m,n

	lambda0=0.87890145
	z=lambda0*s

	h=0.01
	n=1000
	if(z.gt.n*h)write(6,1)s,z
1	format('panic in wake', 2f15.6)
	fc=n*h
	gamma0=ddh0(n)/2
	fa=h0(n)/gamma0
	fb=fa**0.5
	mu0=fb-fc
	np=z/h
	dt=z-np*h
	dt1=1-dt
	ho=h0(np)*dt1+dt*h0(np+1)
	dho=dh0(np)*dt1+dt*dh0(np+1)
	ddho=ddh0(np)*dt1+dt*ddh0(np+1)
	ho=lambda0*ho
	dho=lambda0*lambda0*dho
	ddho=lambda0*lambda0*lambda0*ddho
	return
	end




	subroutine wake0
	implicit real*8(a-h,o-z)
	implicit integer (i-n)
	common/wakedata/h0(0:2000),dh0(0:2000),ddh0(0:2000),h,n
	real *8 a0,b0,df,e
	real *8 fa,fb,fc,fd,dfa,dfb,dfc,dfd,ddfa,ddfb,ddfc,ddfd
	real *8 sp,sm,mu0,mu1,gamma0,gamma1,lambda0
	integer i,j,k,l,m,n

	h=0.01
	n=1000
	h0(0)=0.
	dh0(0)=1.
	ddh0(0)=0.
	do i=1,n
		fa=   h*dh0(i-1)
		dfa=  h*ddh0(i-1)
		ddfa= h*(dh0(i-1)*dh0(i-1)-2*h0(i-1)*ddh0(i-1))/3

		fb=   h*(dh0(i-1)+dfa/2)
		dfb=  h*(ddh0(i-1)+ddfa/2)
		ddfb= h*((dh0(i-1)+dfa/2)*(dh0(i-1)+dfa/2)-
     +                      2*(h0(i-1)+fa/2)*(ddh0(i-1)+ddfa/2))/3

		fc=   h*(dh0(i-1)+dfb/2)
		dfc=  h*(ddh0(i-1)+ddfb/2)
		ddfc= h*((dh0(i-1)+dfb/2)*(dh0(i-1)+dfb/2)-
     +                      2*(h0(i-1)+fb/2)*(ddh0(i-1)+ddfb/2))/3

		fd=   h*(dh0(i-1)+dfc)
		dfd=  h*(ddh0(i-1)+ddfc)
		ddfd= h*((dh0(i-1)+dfc)*(dh0(i-1)+dfc)-
     +                      2*(h0(i-1)+fc)*(ddh0(i-1)+ddfc))/3


		h0(i)=h0(i-1)+(fa+2*fb+2*fc+fd)/6
		dh0(i)=dh0(i-1)+(dfa+2*dfb+2*dfc+dfd)/6
		ddh0(i)=ddh0(i-1)+(ddfa+2*ddfb+2*ddfc+ddfd)/6
	enddo
	return
	end





  	subroutine invert(a,n,np,b)
	implicit real *8(a-h,o-z)
	real*8 a(np,np),b(np)
	integer indx(500)
c
c----------------------------------------------------------------
c
c	You must provide matrix inversion here. For instance
c	using routine 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

