!	calculate the transportation from neighbour cell
	subroutine transport(c,itype, tanalfa, tanbeta, u1,u3,u4,u2,dx,dy,dt,f1,f3,f4,f2)
	
	implicit double precision (a-h,o-z)
	cotalfa=1.0/tanalfa
	cotbeta=1.0/tanbeta
	f1=0
	f3=0
	f4=0
	f2=0
	if(itype.eq.1)then
		s1=0
		s3=sqrt(2.0*c*tanalfa)
		s4=0
		s2=sqrt(2.0*c*tanalfa)
		if(u1.gt.0)then
			if(u1*dt.le.(1.0-s2)*dy)then
				f1=0
			else
				temp1=u1*dt-(1.0-s2)*dy
				f1=0.5*temp1*temp1*cotbeta
			endif
		endif
		if(u2.gt.0)then
			if(u2*dt.le.s3*dx)then
				f2=c*dx*dy
			else
				f2=0.5*u2*dt*(2.0-u2*dt/(s3*dx))*s2*dy
			endif
		endif
		if(u3.lt.0)then
			if((abs(u3)*dt.ge.s2*dy))then
				f3=c*dx*dy
			else
				f3=0.5*abs(u3)*dt*(2.0-abs(u3)*dt/(s2*dy))*s3*dx
			endif
		endif
		if(u4.lt.0)then
			if(abs(u4)*dt.le.(1.0-s3)*dx)then
				f4=0
			else
				temp1=abs(u4)*dt-(1.0-s3)*dx
				f4=0.5*temp1*temp1*tanbeta
			endif
		endif
	else if(itype.eq.2)then
		s1=0
		s3=1.0
		s4=c-0.5*tanalfa
		s2=c+0.5*tanalfa
		if(u1.gt.0)then
			if(u1*dt.le.(1.0-s2)*dy)then
				f1=0
			else if(u1*dt.le.(1.0-s4)*dy)then
				temp1=u1*dt-(1.0-s2)*dy
				f1=0.5*temp1*temp1*cotbeta
			else
				f1=u1*dt*dx-(1.0-c)*dx*dy
			endif

		endif
		if(u2.gt.0)then
			f2=u2*dt*(s2*dy-0.5*u2*dt*tanbeta)
		endif
		if(u3.lt.0)then
			if((abs(u3)*dt.le.s4*dy))then
				f3=abs(u3)*dt*dx
			else if((abs(u3)*dt.le.s2*dy))then
				temp1=abs(u3)*dt-s4*dy
				f3=abs(u3)*dt*dx-0.5*temp1*temp1*cotbeta
			else
				f3=c*dx*dy
			endif
		endif
		if(u4.lt.0)then
			f4=abs(u4)*dt*(s4*dy+0.5*abs(u4)*dt*tanbeta)
		endif
	else if(itype.eq.3)then
		s1=c-0.5*tanalfa
		s3=c+0.5*tanalfa
		s4=0
		s2=1.0
		if(u1.gt.0)then
			f1=u1*dt*(s1*dx+0.5*u1*dt*cotbeta)
		endif
		if(u2.gt.0)then
			if(u2*dt.le.s1*dx)then
				f2=u2*dt*dy
			else if(u2*dt.le.s3*dx)then
				temp1=u2*dt-s1*dx
				f2=u2*dt*dy-0.5*temp1*temp1*tanbeta
			else
				f2=c*dx*dy
			endif
		endif
		if(u3.lt.0)then
			f3=abs(u3)*dt*(s3*dy+0.5*abs(u3)*dt*cotbeta)
		endif
		if(u4.lt.0)then
			if((abs(u4)*dt.le.(1.0-s3)*dx))then
				f4=0
			else if((abs(u4)*dt.le.(1.0-s1)*dx))then
				temp1=abs(u4)*dt-(1.0-s3)*dx
				f4=0.5*temp1*temp1*tanbeta
			else
				f4=abs(u4)*dt*dy-(1.0-c)*dx*dy
			endif
		endif
	else if(itype.eq.4)then
		s1=1.0-sqrt(2.0*(1.0-c)*cotalfa)
		s3=1.0
		s4=1.0-sqrt(2.0*(1.0-c)*tanalfa)
		s2=1.0
		if(u1.gt.0)then
			if(u1*dt.ge.(1.0-s4)*dy)then
				f1=u1*dt*dx-(1.0-c)*dx*dy
			else
				f1=u1*dt*(s1*dx+0.5*u1*dt*cotbeta)
			endif
		endif
		if(u2.gt.0)then
			if(u2*dt.le.s1*dx)then
				f2=u2*dt*dy
			else
				temp1=u2*dt-s1*dx
				f2=u2*dt*dy-0.5*temp1*temp1*tanbeta
			endif
		endif
		if(u3.lt.0)then
			if(abs(u3)*dt.le.s4*dy)then
				f3=abs(u3)*dt*dx
			else
				temp1=abs(u3)*dt-s4*dy
				f3=abs(u3)*dt*dx-0.5*temp1*temp1*cotbeta
			endif
		endif
		if(u4.lt.0)then
			if((abs(u4)*dt.ge.(1.0-s1)*dx))then
				f4=abs(u4)*dt*dy-(1.0-c)*dx*dy
			else
				f4=abs(u4)*dt*(s4*dy+0.5*abs(u4)*dt*tanbeta)
			endif
		endif
	else
		write(*,*) 'error in subroutine transport'
		stop
	endif
	end

!	program main
!	solve fluid volume function equation using Young's VOF method
	program main
	
	parameter(isize=100,pi=3.14159265)
	implicit double precision (a-h,o-z)
	double precision c0(0:isize,0:isize),c1(0:isize,0:isize)
	double precision u(0:isize,0:isize),v(0:isize,0:isize)
	double precision x(0:isize),y(0:isize)
	
	open(unit=1,file='d:/you1_0a.plt',status='unknown')
	open(unit=2,file='d:/you1_0b.plt',status='unknown')
	
	write(*,*) 'Input T:'
	read(*,*) supt
	
	h=1.0/isize
	dt=0.1*h
	dx=h
	dy=h
	emikro=1.0e-10
	emk2=0
	emk2=0
	x0=0.5
	y0=0.3
	r1=0.2
	r2=0.15
	
	do i=0,isize
		x(i)=dx*i
		y(i)=dy*i
	enddo
	do i=0,isize
		do j=0,isize
			u(i,j)=-pi*cos(pi*(x(i)-0.5))*sin(pi*(y(j)-0.5))
			v(i,j)=pi*sin(pi*(x(i)-0.5))*cos(pi*(y(j)-0.5))
		enddo
	enddo
	do i=1,isize
		do j=0,isize
!			if(dist((x(i),y(j)),(x0,y0)).le.r1)then
			if(sqrt((x(i)-x0)*(x(i)-x0)+(y(j)-y0)*(y(j)-y0)).le.r1)then
				c0(i,j)=1.0
			else
				c0(i,j)=0
			endif
		enddo
	enddo
	do ll=1,2
		t=0
		it=0
		do while(t.lt.supt)
			do i=1,isize
				do j=1,isize
					c1(i,j)=c0(i,j)
				enddo
			enddo
			do i=1,isize-1
				do j=1,isize-1
					ft=0
					fb=0
					fl=0
					fr=0
					c=c0(i,j)
					ut=v(i,j)
					ub=v(i,j-1)
					ul=u(i-1,j)
					ur=u(i,j)
					if (c.ge.1.0-emikro)then
						if(ut.gt.0)then
							ft=ut*dt*dx*c
						endif
						if(ub.lt.0)then
							fb=abs(ub)*dt*dx*c
						endif
						if(ul.lt.0)then
							fl=abs(ul)*dt*dy*c
						endif
						if(ur.gt.0)then
							fr=ur*dt*dy*c
						endif
					else if(c.gt.0)then
						rnx=(c0(i+1,j+1)+2.0*c0(i+1,j)+c0(i+1,j-1)-c0(i-1,j+1)-2.0*c0(i-1,j)-c0(i-1,j-1))/dx
						rny=(c0(i+1,j+1)+2.0*c0(i,j+1)+c0(i-1,j+1)-c0(i+1,j-1)-2.0*c0(i,j-1)-c0(i-1,j-1))/dy
						if(abs(rny).le.emk2)then
							if(rnx.ge.0)then
								u1=ut
								u2=ur
								u3=ub
								u4=ul
							else
								u1=ut
								u2=-ul
								u3=ub
								u4=-ur
							endif
							f1=0
							f2=0
							f3=0
							f4=0
							if(u1.gt.0)then
								f1=abs(u1)*dt*c*dx
							endif
							if(u3.lt.0)then
								f3=abs(u3)*dt*c*dx
							endif
							if(u2.gt.0)then
								if(abs(u2)*dt.le.c*dx)then
									f2=abs(u2)*dt*dy
								else
									f2=c*dx*dy
								endif
							endif
							if(u4.lt.0)then
								if(abs(u4)*dt.le.(1.0-c)*dx)then
									f4=0
								else
									f4=(abs(u4)*dt-(1.0-c)*dx)*dy
								endif
							endif
							if(rnx.ge.0)then
								ft=f1
								fb=f3
								fr=f2
								fl=f4
							else
								ft=f1
								fb=f3
								fr=f4
								fl=f2
							endif
						else if(abs(rnx).le.emk2)then
							if(rny.ge.0)then
								u1=ut
								u2=ur
								u3=ub
								u4=ul
							else
								u1=-ub
								u2=ur
								u3=-ut
								u4=ul
							endif
							f1=0
							f2=0
							f3=0
							f4=0
							if(u1.gt.0)then
								if(u1*dt.le.c*dy)then
									f1=u1*dt*dx
								else
									f1=c*dx*dy
								endif
							endif
							if(u3.lt.0)then
								if(abs(u3)*dt.le.(1.0-c)*dy)then
									f3=0
								else
									f3=(abs(u3)*dt-(1.0-c)*dy)*dx
								endif
							endif
							if(u2.gt.0)then
								f2=u2*dt*c*dy
							endif
							if(u4.lt.0)then
								f4=abs(u4)*dt*c*dy
							endif
							if(rnx.ge.0)then
								ft=f1
								fb=f3
								fr=f2
								fl=f4
							else
								ft=f3
								fb=f1
								fr=f2
								fl=f4
							endif
						else
							if(rnx.gt.0.and.rny.lt.0)then
								u1=ut
								u2=ur
								u3=ub
								u4=ul
							else if(rnx.lt.0.and.rny.gt.0)then
								u1=-ub
								u2=-ul
								u3=-ut
								u4=-ur
							else if(rnx.lt.0.and.rny.lt.0)then
								u1=ut
								u2=-ul
								u3=ub
								u4=-ur
							else
								u1=-ub
								u2=ur
								u3=-ut
								u4=ul
							endif
							rnx1=abs(rnx)
							rny1=-abs(rny)
							tanbeta=-rnx1/rny1
							cotbeta=1.0/tanbeta
							tanalfa=dx/dy*tanbeta
							cotalfa=1.0/tanalfa
							if(tanalfa.le.1.0)then
								if(c.le.0.5*tanalfa)then
									itype=1
								else if(c.le.1.0-0.5*tanalfa)then
									itype=2
								else
									itype=4
								endif
							else
								if(c.le.0.5*cotalfa)then
									itype=1
								else if(c.le.1.0-0.5*cotalfa)then
									itype=3
								else
									itype=4
								endif
							endif
							call transport(c,itype, tanalfa, tanbeta, u1,u3,u4,u2,dx,dy,dt,f1,f3,f4,f2)
							if(rnx.gt.0.and.rny.lt.0)then
								ft=f1
								fb=f3
								fl=f4
								fr=f2
							else if(rnx.lt.0.and.rny.gt.0)then
								ft=f3
								fb=f1
								fl=f2
								fr=f4
							else if(rnx.lt.0.and.rny.lt.0)then
								ft=f1
								fb=f3
								fl=f2
								fr=f4
							else
								ft=f3
								fb=f1
								fl=f4
								fr=f2
							endif
						endif
					endif
					c1(i,j+1)=c1(i,j+1)+ft/dx/dy
					c1(i,j-1)=c1(i,j-1)+fb/dx/dy
					c1(i+1,j)=c1(i+1,j)+fr/dx/dy
					c1(i,j)=c1(i,j)-(ft+fb+fl+fr)/dx/dy
				enddo
			enddo
			do i=1,isize
				do j=1,isize
					c0(i,j)=max(min(c1(i,j),1.0),0.0)
				enddo
			enddo
			t=t+1
			it=it+1
		enddo
		write(ll,*) 'TITLE="EXANPLE:3D GEOMETRIES"'
		write(ll,*) 'VARIABLES="X","Y","Z"'
		write(ll,*) 'ZONE T="Floor",I=',isize,'J=',isize,'F=POINT'
		do i=1,isize
			do j=1,isize
				write(ll,102) x(i),y(j),c0(i,j)
			enddo
		enddo
		do i=1,isize
			do j=1,isize
				u(i,j)=pi*cos(pi*(x(i)-0.5))*sin(pi*(y(j)-0.5))
				v(i,j)=-pi*sin(pi*(x(i)-0.5))*cos(pi*(y(j)-0.5))
			enddo
		enddo
	enddo
102	format(1x,e20.10,e20.10,e20.10)
	close(1)
	close(2)
	end