program thom
c	
c	implementation of Thom's stream function-vorticity iteration(sort of)
c
	parameter (ix1=-500,ix2=2000,iy=400)
	complex *8 z(ix1:ix2,0:iy),t(ix1:ix2,0:iy),zp
	real *8 x(ix1:ix2,0:iy),y(ix1:ix2,0:iy)
	real *8 psi(ix1:ix2,0:iy),zeta(ix1:ix2,0:iy)
        real *8 q(ix1:ix2,0:iy),cp(ix1:ix2,0:iy),cpinf(ix1:ix2)
	real *8 th(ix1:ix2),v(ix1:ix2,0:iy)
	real *8 newpsi,inertia
	character *30 file_root,ext,temp
	character *30 x_file,y_file,psi_file,speed_file,vort_file
	character *30 cp_file,cp_cyl_file,log_file
	integer x_unit,y_unit,psi_unit,vort_unit,speed_uni
	integer cp_unit,cp_cyl_unit,log_unit,reynolds
	pi=3.14159265358980
c-------------------------------------------------------------
c   read starting data
c
c       mesh is m1:m2 in xi direction, 0:n in eta direction
c
c	om is relation parameter, if used
c       kiters is number of Reynolds numbers to be iterated
c       iterations is base number
c       r is final Reynolds number
c-------------------------------------------------------------
	read(5,*)m1
	read(5,*)m2
	read(5,*)n
	read(5,*)om
	read(5,*)dxi
	read(5,*)kiters
	read(5,*)iterations
	read(5,*)r
	read(5,*)file_root
	call termstr(file_root,3)
	write(6,90)file_root
90	format('file root =:',a30)
	mp=1.000001/dxi
	deta=dxi
	deta2=deta*deta
c-------------------------------------------------------------
c	generate t,z coordinate fields
c-------------------------------------------------------------
	do i=m1,m2
		do j=0,n
			t(i,j)=cmplx(i*dxi,j*deta)
		if(i.ge. 0)then
			z(i,j)=t(i,j)+sqrt(t(i,j)*t(i,j)-1.)
		else
			z(i,j)=t(i,j)-sqrt(t(i,j)*t(i,j)-1.)
		endif
			x(i,j)=real(z(i,j))
			y(i,j)=aimag(z(i,j))
			q(i,j)=(0.5*abs(1.-1./(z(i,j)*z(i,j))))**2
			zeta(i,j)=0.
			psi(i,j)=2*j*deta
		enddo
		if(y(i,0).lt.0.)y(i,0)=-y(i,0)
	enddo
	th(-mp)=0.
	do i=-mp+1,mp-1
		th(i)=acos(-i*dxi)
	enddo
	th(mp)=pi
c	file_root="test\0"
	x_file="\0"
	call stradd(file_root,"x\0",x_file)
	x_unit=20
	y_file="\0"
	call stradd(file_root,"y\0",y_file)
	y_unit=21
	write(6,*)x_file,y_file
	open(unit=x_unit,file=x_file)
	write(x_unit,20)((x(i,j),i=m1,m2),j=0,n)
	close(unit=x_unit)
	open(unit=y_unit,file=y_file)
	write(y_unit,20)((y(i,j),i=m1,m2),j=0,n)
	close(unit=y_unit)
c-------------------------------------------------------------
c       open log file
	log_unit=19
	call stradd(file_root,"log\0",log_file)
	open(unit=log_unit,file=log_file)
	write(log_unit,5)m1,m2,n,om,dxi,kiters,iterations,r
5	format('Thom method log file',/,
     +         'm1 = ',i6,' m2 = ',i6,' n = ',i6,
     +         'omega      = ',f10.3,/,
     +         'dxi        = ',f10.4,/,
     +         'kiters     = ',i10,/,
     +         'iterations = ',i10,/,
     +         'final r    = ',f10.0)
c-------------------------------------------------------------
	lp=iterations/500
	lscreen=iterations/500
	lpindex=1
	if(lp.eq.0)lp=1
	if(lscreen.eq.0)lscreen=1
c-------------------------------------------------------------
c main loop for changing reynolds number
c-------------------------------------------------------------
	do kit=1,kiters
		r1=(r*kit)/kiters
		om1=(1+kiters-kit)*om
		if(om1.gt.1.)om1=1.
		iterations1=5000+((iterations-5000)*kit)/kiters	
c-------------------------------------------------------------
c      set up output files
c-------------------------------------------------------------
	reynolds=r1
	ext="\0"
	psi_file="\0"
	psi_unit=22
	vort_file="\0"
	vort_unit=23
	speed_file="\0"
	speed_unit=24
	cp_file="\0"
	cp_unit=25
	cp_cyl_file="\0"
	cp_cyl_unit=26
	temp="\0"
	call writestr(ext,reynolds)
	call stradd(file_root,ext,temp)
	call stradd(temp,"psi\0",psi_file)
	open(unit=psi_unit,file=psi_file)
	call stradd(temp,"vort\0",vort_file)
	open(unit=vort_unit,file=vort_file)
	call stradd(temp,"speed\0",speed_file)
	open(unit=speed_unit,file=speed_file)
	call stradd(temp,"cp\0",cp_file)
	open(unit=cp_unit,file=cp_file)
	call stradd(temp,"cpcyl\0",cp_cyl_file)
	open(unit=cp_cyl_unit,file=cp_cyl_file)
	
c-------------------------------------------------------------
c        iterate solution at the Reynolds  number r1
c	 assumes dxi=deta
c-------------------------------------------------------------
	 do k=1,iterations
	  do i=m1+1,m2-1
	    do j=1,n-1
c-------------------------------------------------------------
c      calculate cell centred values, labelled a,b,c,d clockwise
c      from NE corner
c-------------------------------------------------------------
		i1=i
		j1=j
		qbar=(q(i1,j1)+q(i1+1,j1)+
     +                 q(i1,j1+1)+q(i1+1,j1+1))/4
		wbar=(zeta(i1,j1)+zeta(i1+1,j1)+
     +                  zeta(i1,j1+1)+zeta(i1+1,j1+1))/4
		psibar=(psi(i1,j1)+psi(i1+1,j1)+
     +                  psi(i1,j1+1)+psi(i1+1,j1+1))/4
		inertia=r1*((psi(i1+1,j1+1)-psi(i1,j1))*
     +                      (zeta(i1+1,j1)-zeta(i1,j1+1))+
     +                      (psi(i1+1,j1)-psi(i1,j1+1))*
     +                      (zeta(i1,j1)-zeta(i1+1,j1+1)))/8  
	        wa=wbar-inertia/2
		psia=psibar+deta2*wbar/(8*qbar)
c-----------------
c	SE corner
c-----------------
		i1=i
		j1=j-1
		qbar=(q(i1,j1)+q(i1+1,j1)+
     +                 q(i1,j1+1)+q(i1+1,j1+1))/4
		wbar=(zeta(i1,j1)+zeta(i1+1,j1)+
     +                  zeta(i1,j1+1)+zeta(i1+1,j1+1))/4
		psibar=(psi(i1,j1)+psi(i1+1,j1)+
     +                  psi(i1,j1+1)+psi(i1+1,j1+1))/4
		inertia=r1*((psi(i1+1,j1+1)-psi(i1,j1))*
     +                      (zeta(i1+1,j1)-zeta(i1,j1+1))+
     +                      (psi(i1+1,j1)-psi(i1,j1+1))*
     +                      (zeta(i1,j1)-zeta(i1+1,j1+1)))/8  
	        wb=wbar-inertia/2
		psib=psibar+deta2*wbar/(8*qbar)
c-----------------
c  	SW corner
c-----------------
		i1=i-1
		j1=j-1
		qbar=(q(i1,j1)+q(i1+1,j1)+
     +                 q(i1,j1+1)+q(i1+1,j1+1))/4
		wbar=(zeta(i1,j1)+zeta(i1+1,j1)+
     +                  zeta(i1,j1+1)+zeta(i1+1,j1+1))/4
		psibar=(psi(i1,j1)+psi(i1+1,j1)+
     +                  psi(i1,j1+1)+psi(i1+1,j1+1))/4
		inertia=r1*((psi(i1+1,j1+1)-psi(i1,j1))*
     +                      (zeta(i1+1,j1)-zeta(i1,j1+1))+
     +                      (psi(i1+1,j1)-psi(i1,j1+1))*
     +                      (zeta(i1,j1)-zeta(i1+1,j1+1)))/8  
	        wc=wbar-inertia/2
		psic=psibar+deta2*wbar/(8*qbar)
c-----------------
c        NW corner
c-----------------
		i1=i-1
		j1=j
		qbar=(q(i1,j1)+q(i1+1,j1)+
     +                 q(i1,j1+1)+q(i1+1,j1+1))/4
		wbar=(zeta(i1,j1)+zeta(i1+1,j1)+
     +                  zeta(i1,j1+1)+zeta(i1+1,j1+1))/4
		psibar=(psi(i1,j1)+psi(i1+1,j1)+
     +                  psi(i1,j1+1)+psi(i1+1,j1+1))/4
		inertia=r1*((psi(i1+1,j1+1)-psi(i1,j1))*
     +                      (zeta(i1+1,j1)-zeta(i1,j1+1))+
     +                      (psi(i1+1,j1)-psi(i1,j1+1))*
     +                      (zeta(i1,j1)-zeta(i1+1,j1+1)))/8  
	        wd=wbar-inertia/2
		psid=psibar+deta2*wbar/(8*qbar)
c-----------------
c   centre
c-----------------
		wbar=(wa+wb+wc+wd)/4
		psibar=(psia+psib+psic+psid)/4
		inertia=r1*((psia-psic)*(wb-wd)+
     +                      (psib-psid)*(wc-wa))/8  
	        zeta(i,j)=(1-om)*zeta(i,j)+om*(wbar-inertia/2)		
		psi(i,j)=psibar+deta2*wbar/(8*q(i,j))
	    enddo
	  enddo
	  do i=-mp,mp
		zeta(i,0)=-2*psi(i,1)*q(i,0)/deta2
	  enddo
c------------------------------------------------------------
c        output data to screen if k is multiple of lscreen
c------------------------------------------------------------
	 if((k-(k/lscreen)*lscreen).eq.0)then
c------------------------------------------------------------
c	determine length, width of any vortex
c------------------------------------------------------------
	ilength=0
	do i=m2-1,mp+1,-1
		sp=psi(i+1,1)-psi(i+1,0)
		sm=psi(i-1,1)-psi(i-1,0)

		if(sp*sm.lt.0.)then
			if(ilength.eq.0)ilength=i
		endif
	enddo
	vortex_length=x(ilength,0)
	iwidth=0
	jwidth=0
	do i=0,m2
		do j=1,n-1
			sp=psi(i,j+1)
			sm=psi(i,j-1)
			if(sm*sp.lt.0.)then
				if(j.gt.jwidth)then
				  jwidth=j
				  iwidth=i
			       endif
			endif
		enddo
	enddo
	vortex_width=y(iwidth,jwidth)
	iangle=0
	do i=-mp+1,mp-1
		sp=psi(i+1,1)-psi(i+1,0)
		sm=psi(i-1,1)-psi(i-1,0)
		if(sp*sm.lt.0.)then
			if(iangle.eq.0)iangle=i
		endif
	enddo
	angle=180*acos(-iangle*dxi)/3.14159
c------------------------------------------------------------
c       length & width calculated
c------------------------------------------------------------
	  write(6,10)r1,k,zeta(-mp/2,0),zeta(0,0),zeta(mp/2,0),
     +     vortex_width,vortex_length,angle
10	  format(f8.0,i8,3f8.3,3f15.6)

	endif
c------------------------------------------------------------
c       output check value to file
c------------------------------------------------------------
	  if((k-(k/lp)*lp).eq.0)then
		write(11,18)lpindex,zeta(mp/2,0)
		lpindex=lpindex+1
18	        format(i10,f15.6)
	  endif
c------------------------------------------------------------
	enddo
c------------------------------------------------------------
c       end of iteration loop
c------------------------------------------------------------
c
c	determine length of any vortex
c------------------------------------------------------------
	ilength=0
	do i=m2-1,mp+1,-1
		sp=psi(i+1,1)-psi(i+1,0)
		sm=psi(i-1,1)-psi(i-1,0)

		if(sp*sm.lt.0.)then
			if(ilength.eq.0)ilength=i
		endif
	enddo
	vortex_length=x(ilength,0)
	iwidth=0
	jwidth=0
	do i=0,m2
		do j=1,n-1
			sp=psi(i,j+1)
			sm=psi(i,j-1)
			if(sm*sp.lt.0.)then
				if(j.gt.jwidth)then
				  jwidth=j
				  iwidth=i
			       endif
			endif
		enddo
	enddo
	vortex_width=y(iwidth,jwidth)
	iangle=0
	do i=-mp+1,mp-1
		sp=psi(i+1,1)-psi(i+1,0)
		sm=psi(i-1,1)-psi(i-1,0)
		if(sp*sm.lt.0.)then
			if(iangle.eq.0)iangle=i
		endif
	enddo
	angle=180*acos(-iangle*dxi)/3.14159
c-------------------------------------------------------------
c	calculate pressure
c-------------------------------------------------------------
	cp(m1,0)=0.
	cpinf(m1)=0.
c	calculate pressure on surface
	sum=0.
	do i=m1+1,m2
		u=sqrt(q(i,0))*psi(i,1)/deta
		if((i.le.mp).and.(i.ge.-mp))u=0.
		cpinf(i)=1-u*u
		sum=sum+0.5*(zeta(i,1)-zeta(i,0)+
     +                       zeta(i-1,1)-zeta(i-1,0))
		cp(i,0)=1-u*u-2*sum/r1
	enddo
c	calculate pressure in field
	do i=m1+1,m2-1
	sum1=0.
	sum2=0.
	  do j=1,n-1
		ua=0.5*(psi(i,j+1)-psi(i,j-1))/deta
		ub=0.5*(psi(i-1,j)-psi(i+1,j))/deta
		u2=q(i,j)*(ua*ua+ub*ub)
		sum1=sum1+(zeta(i+1,j)-zeta(i-1,j)+
     +                       zeta(i+1,j-1)-zeta(i-1,j-1))/4
		sum2=sum2+(zeta(i,j)+zeta(i,j-1))*(psi(i,j)-psi(i,j-1))/2
		cp(i,j)=cp(i,0)-u2+2*sum1/r1-2*sum2
	  enddo
	  cp(i,n)=cp(i,n-1)
	enddo
	do j=1,n
		cp(m1,j)=0.
		cp(m2,j)=cp(m2-1,j)
	enddo
c	do j=1,n-1
c	cp(m1,j)=0.
c	sum1=0.
c	sum2=0.
c	do i=m1+1,m2-1
c		ua=0.5*(psi(i,j+1)-psi(i,j-1))/deta
c		ub=0.5*(psi(i-1,j)-psi(i+1,j))/deta
c		u2=q(i,j)*(ua*ua+ub*ub)
c		sum1=sum1+(zeta(i,j+1)-zeta(i,j-1)+
c     +                       zeta(i-1,j+1)-zeta(i-1,j-1))/4
c		sum2=sum2+(zeta(i,j)+zeta(i-1,j))*(psi(i,j)-psi(i-1,j))/2
c		cp(i,j)=1-u2-2*sum1/r1-2*sum2
c	enddo
c	cp(m2,j)=cp(m2-1,j)
c	enddo
	do i=m1,m2
		cp(i,n)=cp(i,n-1)
	enddo
	write(45,20)((cp(i,j),i=m1,m2),j=0,n)
	write(cp_unit,25)(i*dxi,cp(i,0),i=m1,m2)
	write(cp_cyl_unit,25)(180*acos(-i*dxi)/3.14159,cp(i,0),i=-mp,mp)
c------------------------------------------------------------
c	calculate pressure and viscous drag coefficients
c------------------------------------------------------------
	cpp=0.
	cpv=0.
	do i=-mp+1,mp
		ang=(th(i)+th(i-1))/2
		pressure=(cp(i,0)+cp(i-1,0))/2
		dth=th(i)-th(i-1)
		cpp=cpp+pressure*cos(ang)*dth
		cpv=cpv-2*zeta(i,0)*sin(ang)*dth/r1
	enddo
	cd=cpp+cpv
	write(6,60)vortex_width,vortex_length,angle,cpp,cpv,cd
	write(log_unit,60)vortex_width,vortex_length,angle,cpp,cpv,cd
	
c-------------------------------------------------------------
c	calculate speed
c-------------------------------------------------------------
	do i=m1+1,m2-1
	       us=(psi(i,1)-psi(i,0))/deta
	       vs=-0.5*(psi(i+1,0)-psi(i-1,0))/dxi
		v(i,0)=sqrt(q(i,0)*(us*us+vs*vs))
		do j=1,n-1
	           us=0.5*(psi(i,j+1)-psi(i,j-1))/deta
	           vs=-0.5*(psi(i+1,j)-psi(i-1,j))/dxi
		   v(i,j)=sqrt(q(i,j)*(us*us+vs*vs))
		enddo
	        v(i,n)=v(i,n-1)
	enddo
	do j=1,n
		v(m1,j)=v(m1+1,j)
		v(m2,j)=v(m2-1,j)
	enddo
	write(psi_unit  ,20)((psi(i,j),  i=m1,m2),j=0,n)
	write(vort_unit ,20)((zeta(i,j), i=m1,m2),j=0,n)
	write(speed_unit,20)((v(i,j),i=m1,m2),j=0,n)
	close(unit=psi_unit)
	close(unit=vort_unit)
	close(unit=speed_unit)
	close(unit=cp_unit)
	close(unit=cp_cyl_unit)

c--------------------------------------------------------------
c end of reynolds number loop
c--------------------------------------------------------------
        enddo
	close(unit=log_unit)
	kcontour=2
20	format(8f12.4)
25	format('data',/,(2f15.6))
30	format('Cpressure=',f15.6,' Cviscous=',f15.6,' Cd=',f15.6)
40	format(8f12.4)
60      format('width = ',f10.4,' length = ',f10.4, 
     +                          ' angle = ',f10.4,3f8.4)

	


	stop
	end