	program freeint
c
c	free interaction in a parallel channel
c
c	have to set the initial value of a
c
c	internal channel problem
c	
c
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      parameter (iy=50,iyp=iy+1,iz=2*3*(iy+1)+3,ix=100,iyg=128)
c
c	variables on upper wall
c
	real *8 uu(0:ix,0:iy),duu(0:iyp)
	real *8 wu(0:ix,0:iy),dwu(0:iyp)
	real *8 psiu(0:ix,0:iy),dpsiu(0:iyp)
	real *8 vu(0:ix,0:iy)
	real *8 pu(0:ix),pxu(0:ix),dpxu
c
c	variables on lower wall
c
	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 f(0:ix),g(0:ix)
	real *8 a(-3:ix+2),ad(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 xs(0:ix,0:iy),ys(0:ix,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
	

c
c	need matrices to hold physical coordinates and velocities
c
	real *8 yg(0:ix,0:iyg),dyg(0:ix)
	real *8 yl(0:ix,0:iy),yu(0:ix,0:iy)
	real *8 uo(0:ix,0:iyg),psio(0:ix,0:iyg)
	real *8 ucomp(0:ix,0:iyg),uinner(0:ix,0:iy),uloc(0:iy)
	real *8 upper(0:ix,0:iy),yloc(0:iy),flux(0:ix)
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 ix,ipx,ipy,iy,i,j,k,msize,msizep,row
c	set up initial guess

	c=5.73
	tolerance=1.e-6
	ipx=1
	ipy=1

	iterations=1
c
c	set up mesh
c
	hy0=.01
	hx0=0.013
	facy=1.1
	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)

	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)
	write(8,7)(y(j),hy(j),j=0,iy)
7	format(2e15.5)
	
	write(6,5)ix,iy,zm,xm,hx0,hy0,msize,msizep
5	format('Interactive channel solver',/, 'ix = ',i4,' iy  = ',i4,
     +   /,'zm =',f8.0, 'xm = ', f12.4,/,
     +     '  hx0 = ',f10.4, ' hy0 = ',f10.4,2i10)
        write(6,10)k,var,x(0),x(ix/4),
     +      x(ix/2),x(3*ix/4),x(ix),p(0)

 

c	set up initial guess for a(x),p(x)
	write(6,1)
1	format('Initialising a(x)')
c
	do i=0,ix
		do j=0,iy
			xs(i,j)=x(i)
			ys(i,j)=y(j)
		enddo
		a(i)=0.
		p(i)=0
		px(i)=0
		pu(i)=0
		pxu(i)=0
	end do
	a(0)=-0.0001

	a(-1)=a(0)*exp(-c*hx0)
	a(-2)=a(-1)*exp(-c*hx0)
	a(-3)=a(-2)*exp(-c*hx0)
	a(ix+1)=0.
	a(ix+2)=0.

	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./2

		uu(i,0)=0.
		psiu(i,0)=0.
		wu(i,0)=1./2

		do j=1,iy
			z=y(j)
			psi(i,j)=z*z/4
			u(i,j)=z/2
			w(i,j)=1./2
			psiu(i,j)=z*z/4
			uu(i,j)=z/2
			wu(i,j)=1./2
		end do
	end do
c
c	set pressure perturbation
c
	p(0)=0.0

c	now have initial guess for all variables,
c	output intial values
c
	write(9,21) ((xs(i,j),j=0,iy,ipy),i=0,ix,ipx)
	write(10,21)((ys(i,j),j=0,iy,ipy),i=0,ix,ipx)
	write(13,21)((psi(i,j),j=0,iy,ipy),i=0,ix,ipx)
	write(14,21)((u(i,j),j=0,iy,ipy),i=0,ix,ipx)
	write(15,65)(x(i),a(i),  i=0,ix)
	write(16,65)(x(i),p(i),  i=0,ix)
	write(17,65)(x(i),ad(i), i=0,ix)
	write(18,65)(x(i)-0.5*hx(i),px(i),i=0+1,ix)
	write(19,66)(y(j),u(0,j),j=0,iy)
	write(19,66)(y(j),u(0,j),j=0,iy)
	write(19,66)(y(j),u(ix/4,j),j=0,iy)
	write(19,66)(y(j),u(ix/2,j),j=0,iy)
	write(19,66)(y(j),u(3*ix/4,j),j=0,iy)
	write(19,66)(y(j),u(ix,j),j=0,iy)
	write(20,67)(x(i),w(i,0),i=0,ix)
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=1,ix
		do j=1,iy
			psi(i,j)=psi(i-1,j)
			u(i,j)=u(i-1,j)
			w(i,j)=w(i-1,j)
			psiu(i,j)=psiu(i-1,j)
			uu(i,j)=uu(i-1,j)
			wu(i,j)=wu(i-1,j)
		enddo
c--------------------------------------------------------------------
c       iterate along channel rows - u(0)=0                         |
c--------------------------------------------------------------------
	   do newton=1,5
		do l=1,iz
		do m=1,iz
			cof(l,m)=0.
		enddo
		enddo
c
c	lower layer
c
		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)   =  -2*hy(j)
			rhs(row)=2*hy(j)*(px(i)
     +					+ uc*ux - wc*phix) -wy
		enddo

		row=row+1
		cof(row,row-1)=-1.
		cof(row,iz)=0.5
		rhs(row)=u(i,iy)-(y(iy)+a(i))/2

		row=row+1
		cof(row,row-1)=1.
		rhs(row)=0
c
c	now do upper layer
c

		row=row+1
c		set first two rows
		cof(row,row-1)=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=(wu(i,j)+wu(i-1,j)+
     +                      wu(i,j-1)+wu(i-1,j-1))/4
			uc=(uu(i,j)+uu(i-1,j)+
     +                      uu(i,j-1)+uu(i-1,j-1))/4
			phiy=psiu(i,j)+  psiu(i-1,j)-
     +                       psiu(i,j-1)-psiu(i-1,j-1)
			uy=    uu(i,j)+    uu(i-1,j)-
     +                         uu(i,j-1)-  uu(i-1,j-1)
			wy=    wu(i,j)+    wu(i-1,j)-
     +                         wu(i,j-1)-  wu(i-1,j-1)
		        phix=(psiu(i,  j)+psiu(i,  j-1)-
     +                        psiu(i-1,j)-psiu(i-1,j-1))/(2*hx(i))
		        ux=(    uu(i,  j)+  uu(i,  j-1)-
     +                          uu(i-1,j)-  uu(i-1,j-1))/(2*hx(i))

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

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

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

		row=row+1
		cof(row,row-2)=-1.
		cof(row,iz)=-0.5
		rhs(row)=uu(i,iy)-(y(iy)-a(i))/2

		row=row+1
		cof(row,row-2)=1.
		rhs(row)=0
c
c	pressure relationship: row should be equal to iz
c
		addd=(a(i)-3*a(i-1)+3*a(i-2)-a(i-3))/(hx(i)**3)
		row=row+1
		cof(row,row-2)=1.
		cof(row,row-1)=-1.
		cof(row,row)= 1./(120.*hx(i)**3)
		
		rhs(row)=pxu(i)-px(i)-addd/120
		

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(ix,iy)
     +          	       +abs(du(j))/u(ix,iy)
     +           	       +abs(dw(j))/w(ix,iy)
		enddo
		do j=0,iy
			key=3*(iy+1)+3*j
			dpsiu(j)=rhs(key+1)
			psiu(i,j)=psiu(i,j)+dpsiu(j)
			duu(j)=rhs(key+2)
			uu(i,j)=uu(i,j)+duu(j)
			dwu(j)=rhs(key+3)
			wu(i,j)=wu(i,j)+dwu(j)
			error=error+abs(dpsiu(j))/psi(ix,iy)
     +          	       +abs(duu(j))/u(ix,iy)
     +           	       +abs(dwu(j))/w(ix,iy)

		enddo
c		write(6,72)error/(3*iy)
72		format('plate',e12.3)
		dpx=rhs(iz-2)
		px(i)=px(i)+dpx
		dpxu=rhs(iz-1)
		pxu(i)=pxu(i)+dpxu
		a(i)=a(i)+rhs(iz)
	   enddo
        enddo
c
c ----------------	impose periodicity
c
	a(ix+1)=a(ix)
	a(ix+2)=a(ix)
c
c --------------	estimate a'' at left boundary
c
	i=0
	add=(a(i+1)-2*a(i)+a(i-1))/(hx(i+1)**2)
c
c	calculate pressure on upper wall relative to lower wall
c	
	pu(0)=p(0)+add/120
	do i=1,ix
		p(i)=p(i-1)+hx(i)*px(i)
		pu(i)=pu(i-1)+hx(i)*pxu(i)
	enddo
        suma=0.
        do i=0,ix
                suma=suma+a(i)*a(i)
        enddo
        wlmax=w(0,0)
        wlmin=w(0,0)
        wumax=wu(0,0)
        wumin=wu(0,0)
        do i=1,ix
                if(w(i,0).gt.wlmax)wlmax=w(i,0)
                if(w(i,0).lt.wlmin)wlmin=w(i,0)
                if(wu(i,0).gt.wumax)wumax=wu(i,0)
                if(wu(i,0).lt.wumin)wumin=wu(i,0)
        enddo

        write(6,10)k,wlmin,wlmax,wumin,wumax,p(0),pu(0),a(0),suma

10      format(i6,(10e12.3))
c-----------------------------------------------------------------------------------
c                                                             end of iteration loop
c-----------------------------------------------------------------------------------
	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(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)

	write(29,66)(y(j),u(0,j),j=0,iy)
	write(29,66)(y(j),u(0,j),j=0,iy)
	write(29,66)(y(j),u(ix/4,j),j=0,iy)
	write(29,66)(y(j),u(ix/2,j),j=0,iy)
	write(29,66)(y(j),u(3*ix/4,j),j=0,iy)
	write(29,66)(y(j),u(ix,j),j=0,iy)

	write(30,67)(x(i),w(i,0),i=0,ix)
	write(50,67)(x(i),wu(i,0),i=0,ix)
	write(46,67)(x(i),pu(i),i=0,ix)

	write(42,21)((psiu(i,j),j=0,iy,ipy),i=0,ix,ipx)
	write(43,21)((uu(i,j),j=0,iy,ipy),i=0,ix,ipx)
	write(44,21)((wu(i,j),j=0,iy,ipy),i=0,ix,ipx)
	



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=0,ix
		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=0,ix
	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=0,ix)
	write(36,66)(xs(i,0)/epsilon,yg(i,iyg),i=0,ix)

	do i=0,ix,ix/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=0,ix)
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 needed here (E.G. from Numerical Recipes)
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
