	program channel
c
c	implement newton method to solve for interactive channel flow
c
        implicit real*8 (a-h,o-z)
        implicit integer (i-n)
        parameter (iy=64,iyp=iy+1,iz=3*(iy+1)+1,ixm=128,ixp=128,iyg=128)
	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(-ixm:ixp),px(-ixm:ixp),dpx
	real *8 a(-ixm:ixp),ad(-ixm:ixp)
	real *8 cof(iz,iz),rhs(iz)
	real *8 x(-ixm:ixp),hx(-ixm:ixp),hy(0:iy),y(0:iy)
	real *8 xs(-ixm:ixp,0:iy),ys(-ixm:ixp,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
c
c	need matrices to hold physical coordinates and velocities
c
	real *8 yg(-ixm:ixp,0:iyg),dyg(-ixm:ixp)
	real *8 yl(-ixm:ixp,0:iy),yu(-ixm:ixp,0:iy)
	real *8 uo(-ixm:ixp,0:iyg),psio(-ixm:ixp,0:iyg)
	real *8 ucomp(-ixm:ixp,0:iyg),uinner(-ixm:ixp,0:iy),uloc(0:iy)
	real *8 upper(-ixm:ixp,0:iy),yloc(0:iy),flux(-ixm:ixp)
c
c	yg= y coordinates in global domain
c	yl= y coordinates in lower wall layer }
c	yu= y coordinates in upper wall layer } in global coordinates, only known when R given
c	ys= Y coordinates
c
	integer ixu,ixd,ipx,ipy,iy,i,j,k,msize,msizep,row
c
c	specify flux through channel
c
	fluxin=1.
	fluxcorrection=fluxin*12
	fluxcorrection=1.
c	set up initial guess
	tolerance=1.e-6
	ixu=ixm
	ixd=ixp
	ipx=1
	ipy=1
	msize=ixm
	msizep=ixp
	iterations=200
c
c	set up mesh
c
	hx0=1./(ixu+ixd)
	hy0=0.05
	facy=1.04
	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)

	facu=1.0
	facd=1.0
	x(0)=0.
	x(1)=hx0
	hx(1)=hx0
	hx(0)=hx0
	x(-1)=-hx0
	do i=2,msizep
		hx(i)=facd*hx(i-1)
		x(i)=x(i-1)+hx(i)
	enddo
	do i=2,msize
		hx(-i+1)=facu*hx(-i+2)
		x(-i)=x(-i+1)-hx(-i+1)

	enddo
	xm=x(ixd)
	write(7,7)(x(i),hx(i),i=-ixu,ixd)
	write(8,7)(y(j),hy(j),j=0,iy)
7	format(2e15.5)
	
	write(6,5)ixu,ixd,iy,zm,xm,hx0,hy0,msize,msizep
5	format('Interactive channel 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)

 

c	set up initial guess for a(x),p(x)
	write(6,1)
1	format('Initialising a(x)')
c
c++++++++++++++++++++++++++++++++++
c                                 +
	write(6,*)' Input wall size factor'
	read(5,*)sfac
c                                 +
c++++++++++++++++++++++++++++++++++
c
	do i=-ixu,ixd
		t=(x(i)-x(-ixu))/(x(ixd)-x(-ixu))
		do j=0,iy
			xs(i,j)=x(i)
			ys(i,j)=-sfac*(1-cos(2*3.14159*t))/2+y(j)
		enddo

c
c	note that a(X) here is B(X)=A(X)+f(X)=f/2=-A if g=0
c
		a(i)=ys(i,0)/2
		p(i)=0
		px(i)=0
	end do
	write(6,3)
3	format('Setting up plate start solution')
	do i=-ixu,ixp
		u(i,0)=0.
		psi(i,0)=0.
		w(i,0)=fluxcorrection/2
		do j=1,iy
			z=y(j)
			psi(i,j)=fluxcorrection*(z*z/2+z*a(i))/2
			u(i,j)=fluxcorrection*(z+a(i))/2
			w(i,j)=fluxcorrection/2
		end do
	end do

c	now have initial guess for all variables,
c	output intial values
c
	write(9,21) ((xs(i,j),j=0,iy,ipy),i=-ixu,ixd,ipx)
	write(10,21)((ys(i,j),j=0,iy,ipy),i=-ixu,ixd,ipx)
	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)(x(i),a(i),  i=-ixu,ixd)
	write(16,65)(x(i),p(i),  i=-ixu,ixd)
	write(17,65)(x(i),ad(i), i=-ixu,ixd)
	write(18,65)(x(i)-0.5*hx(i),px(i),i=-ixu+1,ixd)
	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,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+1,ixp
c--------------------------------------------------------------------
c       iterate along plate rows - u(0)=0                           |
c--------------------------------------------------------------------
	   do newton=1,5
		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
	do i=-ixu+1,ixd
		p(i)=p(i-1)+hx(i)*px(i)
	enddo


c
c ----------------	impose periodicity
c
	do j=0,iy
		psi(-ixu,j)=psi(ixd,j)
		u(-ixu,j)=u(ixd,j)
		w(-ixu,j)=w(ixd,j)
	enddo
	wmin=100.
	do i=-ixu,ixd
		if(w(i,0).lt.wmin)wmin=w(i,0)
	enddo

        write(6,10)k,wmin,w(-ixu/2,0),w(0,0),w(ixd/4,0),
     +      w(ixd/2,0),w(3*ixd/4,0),w(ixd,0),p(0),px(0)
10      format(i6,(10e12.3))
c-----------------------------------------------------------------------------------
c                                                             end of iteration loop
c-----------------------------------------------------------------------------------
	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)(x(i),      a(i),  i=-ixu,ixd)
	write(26,65)(x(i),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,ixd)

21	format(8e15.4)
	do kk=1,10

c
c	loop up to 10 times
c
		write(6,45)
45	format('looping various composite expansions')
		write(6,46)
46	format('input length of furrows')
		read(5,*)flength
		if(flength.lt.0.)go to 101
		epsilon=1./flength
		write(6,47)
47	format('input depth of furrow')
		read(5,*)fdepth
		fac=fdepth/(2*sfac)
		ain=log(fac)/log(epsilon)
		eps=epsilon**ain
		write(6,48)epsilon,ain,eps
48	format('epsilon ',f10.4, ' ain ',f15.4,' eps ',f15.4)
		re=flength**(1+3*ain)
		write(6,49)re
49	format('Reynolds number = ',f15.3)

c
c	generate composite solution
c
	write(6,50)
50	format('Generating composite expansion')
	eps2=eps*eps
c
c	doing case of flat top wall, sinuous lower wall
c
	do i=-ixu,ixd
		dyg(i)=(1.-eps*ys(i,0))/iyg
		do j=0,iyg
			yg(i,j)=eps*ys(i,0)+j*dyg(i)
			if(yg(i,j).ge. 0.)then
		 	  vo=yg(i,j)*(1.-yg(i,j))/2
			  dvo=0.5-yg(i,j)
c	A(X)=-a(X)
			  uo(i,j)=fluxcorrection*(vo-eps*a(i)*dvo)
			else
			  uo(i,j)=fluxcorrection*(yg(i,j)-eps*a(i))/2
			endif
		enddo
		do j=0,iy
			yl(i,j)=eps*ys(i,j)
			yu(i,j)=1.-eps*y(j)
			zu=1-yu(i,j)
			uinner(i,j)=eps*u(i,j)-
     +                          fluxcorrection*(yl(i,j)/2-eps*a(i)/2)
			upper(i,j)=eps*u(i,j)-
     +                          fluxcorrection*(zu/2+eps*a(i)/2)
		enddo
	enddo
	write(35,66)(yg(0,j),uo(0,j),j=0,iyg)
	write(35,66)(yl(0,j),eps*u(0,j),j=0,iy)
	write(35,66)(yu(0,j),eps*u(0,j),j=0,iy)
	do i=-ixu,ixd
	do j=0,iy
		yloc(j)=yl(i,j)
		uloc(j)=uinner(i,j)
	enddo
	do j=0,iyg
	ucomp(i,j)=uo(i,j)+ainterp(yloc,uloc,iy,yg(i,j))
	enddo
	do j=0,iy
		yloc(j)=yu(i,j)
		uloc(j)=upper(i,j)
	enddo
	do j=0,iyg
	ucomp(i,j)=ucomp(i,j)+ainterp(yloc,uloc,iy,yg(i,j))
	enddo

	
	flux(i)=0.
	do j=1,iyg
		flux(i)=flux(i)+(ucomp(i,j)+ucomp(i,j-1))*(yg(i,j)-yg(i,j-1))/2
	enddo
	enddo
	write(36,66)(xs(i,0)/epsilon,yg(i,0),i=-ixu,ixd)
	write(36,66)(xs(i,0)/epsilon,yg(i,iyg),i=-ixu,ixd)

	do i=-ixu,ixd,ixu/2
	write(36,66)(xs(i,0)/epsilon+5*ucomp(i,j),yg(i,j),j=0,iyg)
	enddo
c	write(6,88)(xs(i,0)/epsilon,flux(i),i=-ixu,ixd)
88	format(' flux = ',(2f12.5))
	write(37,66)(yg(0,j),ucomp(0,j),j=0,iyg)
	enddo
101	continue	
	stop
	end

	real *8 function ainterp(y,f,iy,p)
	real *8 y(0:iy),f(0:iy),p,fac
	integer j,iy
c
c	interpolate data given at (y,f) at point p
c
	ainterp=0.
	
	do j=1,iy
		if((y(j)-p)*(y(j-1)-p).le.0.)then
			fac=(p-y(j-1))/(y(j)-y(j-1))
			ainterp=(1-fac)*f(j-1)+fac*f(j)
		endif
	enddo
	return
	end
c-----------------------------------------------------------------------------------
c                                                                matrix inverter   |
c-----------------------------------------------------------------------------------
  	subroutine invert(a,n,np,b)
	implicit real *8(a-h,o-z)
	real*8 a(np,np),b(np)
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
c****************************************************
c   YOU MUST PROVIDE THIS ROUTINE                ****
c   for instance from Numerical Recipies         ****
c****************************************************
c
	return
	end
