c     ********************************************************
c     **  two-dimensional incompressible viscous flow       **
c     **  in general coordinate                             **
c     **                                                    **
c     **  the 3rd order finite difference scheme            **
c     **                                                    **
c     **  for o-type grid                                   **
c     **                                                    **
c     **  euler forward scheme (semi implicit)              **
c     **                                                    **
c     **  marker and cell method                            **
c     **                                                    **
c     **                                                    **
c     **  unit  1:grid data read                            **
c     **  unit  2:pre-results read                          **
c     **  unit  3:results of instantanious write            **
c     **  unit  7:cd & cl write                             **
c     **  unit 10:history write                             **
c     **                                                    **
c     **  2004/06/28                計算に燃えてる皆様へ    **
c     ********************************************************
c
c
c     ******************
c     ** main routine **
c     ******************
c
      open(2,file='r2.dat')  ! read result
      open(3,file='r1.dat')  ! write result
      open(4,file='a1.dat')  ! time average
c
c
         call init
c
         call input
c
         call vtrans
c
         call solve
c
         call output
c
c
      stop
      end
c
c     ********************************************
c     ** calculate flow by time marching method **
c     ********************************************
c
      subroutine solve
c
      implicit double precision(a-h,o-z)
      parameter(mx=401,my=201)
      common /tmcnd/ time,nstep,iskip
      common /reyno/ re,rei,dt
c
c     ***** time marching loop *****
c
      do 100 istep=1,nstep
c
      if(mod(istep,iskip).eq.0) then
         ichk=0
      else
         ichk=1
      end if
c
         time=time+dt
c
      if(ichk.eq.0) then
      write(*,600) time
  600 format(' time=',f8.4)
      end if
c
c
         call boundv
c
         call navi(ichk)
c
         call boundp
c
         call poiss(ichk)
c
         call cdcl(ichk)
c
         call judge(errs)
c
         if(istep.gt.100 .and. errs.lt.1.d-6) goto 90
c
  100 continue
c
   90 continue
c
c
      return
      end
c
c     *****************************
c     ** subroutine initial data **
c     *****************************
c
      subroutine init
c
      implicit double precision(a-h,o-z)
      common /divide/ nxp1,nyp1,nx,ny,nxm1,nym1
      common /conv1/ eps1,omega1,iter1
      common /conv2/ eps2,omega2,iter2
      common /tmcnd/ time,nstep,iskip
      common /reyno/ re,rei,dt
      common /conti/ icon
c
c     *** iteration parameter ***
c
         iter1=300
         iter2= 50
c
c     *** relaxation parameter ***
c
         omega1=1.9d0
         omega2=1.1d0
c
c     *** reynolds number ***
c
         re =1.d+3
         rei=1.d0/re
c
c     *** convergence parameter ***
c
         eps1=1.d-3
         eps2=1.d-8
c
c     *** pre-result reading yes=1 no=0 ***
c
         icon=0
c
c     *** time step ***
c
         dt=0.002d0
c
         nstep=10000
c
         iskip=10
c
c
      return
      end
c
c     ***************************
c     ** subroutine input data **
c     ***************************
c
      subroutine input
c
      implicit double precision(a-h,o-z)
      parameter(mx=401,my=201)
      common /divide/ nxp1,nyp1,nx,ny,nxm1,nym1
      common /coord/ x(mx,my),y(mx,my)
      common /result/ u(0:mx+1,0:my+1),v(0:mx+1,0:my+1),p(mx,my)
      common /tmcnd/ time,nstep,iskip
      common /conti/ icon
c
c     *** grid data reading ***
c
c      read(1,*) nx,ny
c      read(1,*)((x(i,j),y(i,j),i=1,nx),j=1,ny)
c
  500 format(2i10)
  510 format(2f15.10)
c
c      close(1)
c
c     *** 代数的に生成 ***
c
      pi=4.d0*datan(1.d0)
c
c     *** init ***
c     格子数を変更すると発散することがありますが
c     rtを調��して下さい。べき乗なので，慎重に変更して下さい
c
        nx=181
        ny=121
c
        th=2.d0*pi/dfloat(nx-1)
        rt=1.053d0   ! これで格子分布を調��します
c
c     *** compute coordinate around circular cylinder ***
c
      do 100 j=1,ny
      do 100 i=1,nx
c
        r=.5d0+0.00316d0*(1.d0-rt**dble(j-1))/(1.d0-rt)
        x(i,j)=r*dcos(th*dble(i-1))
        y(i,j)=r*dsin(th*dble(i-1))
c
  100 continue
c
c
c     *** number of divided ***
c
          nxp1= nx+1
          nyp1= ny+1
          nxm1= nx-1
          nym1= ny-1
c
c     *****  impulsive start  *****
c
      do 110 j=1,ny
      do 110 i=1,nx
c
         p(i,j) =0.d0
         u(i,j) =0.d0
         v(i,j) =0.d0
c
  110 continue
c
      if(icon.eq.1) then
c
c     *** pre-results data ***
c
      read(2,600) nx,ny,time,dmy
  600 format(2i10,f10.5,f10.0)
c
      read(2,602)((u(i,j),v(i,j),p(i,j),i=1,nx),j=1,ny)
  602 format(2f10.6,e15.7)
c
      close(2)
c
         call boundp
         call boundv
c
      end if
c
c
      return
      end
c
c     ****************************
c     ** output result to file  **
c     ****************************
c
      subroutine output
c
      implicit double precision(a-h,o-z)
      parameter(mx=401,my=201)
      common /tmcnd/ time,nstep,iskip
      common /reyno/ re,rei,dt
      common /coord/ x(mx,my),y(mx,my)
      common /divide/ nxp1,nyp1,nx,ny,nxm1,nym1
      common /result/ u(0:mx+1,0:my+1),v(0:mx+1,0:my+1),p(mx,my)
      common /mean/ au(mx,my),av(mx,my),ap(mx,my)
      dimension cp(mx),the(0:mx)
c
c
      write(3,600) nx,ny,time,re
  600 format(2i10,f10.5,f10.0)
c
      write(3,610)((u(i,j),v(i,j),p(i,j),i=1,nx),j=1,ny)
  610 format(2f10.6,e15.7)
c
      write(3,615)((x(i,j),y(i,j),i=1,nx),j=1,ny)
  615 format(2f10.5)
c
c
      close(3)
c
c     *** time average ***
c
      do 100 j=1,ny
      do 100 i=1,nx
c
         au(i,j)=au(i,j)/dble(nstep)
         av(i,j)=av(i,j)/dble(nstep)
         ap(i,j)=ap(i,j)/dble(nstep)
c
  100 continue
c
c
      write(4,620) nx,ny,nstep,re
  620 format(3i10,f12.4)
c
      write(4,625)((au(i,j),av(i,j),ap(i,j),i=1,nx),j=1,ny)
  625 format(3e15.7)
c
c     *** pressure coefficient ***
c
        pi=4.d0*datan(1.d0)
c
        dth=2.d0*pi/dble(nx-1)
        the(0)=-dth
c
      do 110 i=1,nx
c
        the(i)=the(i-1)+dth
c
        pref=ap(nx/2+1,ny)
  110   cp(i)=2.d0*(ap(i,1)-pref)
c
c
c      write(10,630)(x(i,1),cp(i),i=1,nx)
      write(10,630)(the(i)*180.d0/pi,cp(i),i=1,nx)
  630 format(2f10.5)
c
c
      return
      end
c
c
c     ***************************
c     ** judge of steady state **
c     ***************************
c
      subroutine judge(errs)
c
      implicit double precision(a-h,o-z)
      parameter(mx=401,my=201)
      common /divide/ nxp1,nyp1,nx,ny,nxm1,nym1
      common /result/ u(0:mx+1,0:my+1),v(0:mx+1,0:my+1),p(mx,my)
      common /uold/ uo(mx,my),vo(mx,my)
      common /tmcnd/ time,nstep,iskip
c
c     *** root mean square error ***
c
         errs=0.d0
c
      do 110 j=2,nym1
      do 110 i=2,nx
c
         rhsu=(u(i,j)-uo(i,j))**2
         rhsv=(v(i,j)-vo(i,j))**2
         errs=errs+rhsu+rhsv
c
  110 continue
c
         errs=dsqrt(errs/dfloat((nx-1)*(ny-2)))
c
c      write(*,*)' ********************************'
c      write(*,*)'     '
c      write(*,600) errs
c      write(*,*)'     '
c      write(*,*)' ********************************'
c
  600 format('   total energy check=',f8.6)
c
c
      return
      end
c
c     ************************************
c     ** independent variable translate **
c     ************************************
c
      subroutine vtrans
c
      implicit double precision(a-h,o-z)
      parameter(mx=401,my=201)
      common /divide/ nxp1,nyp1,nx,ny,nxm1,nym1
      common /coord/ x(mx,my),y(mx,my)
      common /coeff/ c1(mx,my),c2(mx,my),c3(mx,my),c4(mx,my),c5(mx,my)
      common /trans/ xx(mx,my),xe(mx,my),yx(mx,my),ye(mx,my)
      common /jacobi/ rjaci(mx,my)
c
c
      do 110 j=2,nym1
c
         jp1=j+1
         jm1=j-1
c
      do 110 i=2,nx
c
         ip1=i+1
         im1=i-1
c
      if(i.eq.nx) ip1=2
c
c     *****  dx/dxi *****
c
      xx(i,j)=.5d0*(x(ip1,j)-x(im1,j))
c
c     *****  dx/deta  *****
c
      xe(i,j)=.5d0*(x(i,jp1)-x(i,jm1))
c
c     *****  dy/dxi  *****
c
      yx(i,j)=.5d0*(y(ip1,j)-y(im1,j))
c
c     *****  dy/deta  *****
c
      ye(i,j)=.5d0*(y(i,jp1)-y(i,jm1))
c
c     *****  jacobian  *****
c
      rjac=xx(i,j)*ye(i,j)-xe(i,j)*yx(i,j)
      rjaci(i,j)=1.d0/rjac
c
c     *** 2nd order derivative ***
c
      xxx=x(ip1,j)-2.d0*x(i,j)+x(im1,j)
      xee=x(i,jp1)-2.d0*x(i,j)+x(i,jm1)
      yxx=y(ip1,j)-2.d0*y(i,j)+y(im1,j)
      yee=y(i,jp1)-2.d0*y(i,j)+y(i,jm1)
      xxe=0.25d0*(x(ip1,jp1)-x(im1,jp1)-x(ip1,jm1)+x(im1,jm1))
      yxe=0.25d0*(y(ip1,jp1)-y(im1,jp1)-y(ip1,jm1)+y(im1,jm1))
c
      alpha=xe(i,j)*xe(i,j)+ye(i,j)*ye(i,j)
      beta =xx(i,j)*xe(i,j)+yx(i,j)*ye(i,j)
      gamma=xx(i,j)*xx(i,j)+yx(i,j)*yx(i,j)
c
      rjacs=rjaci(i,j)*rjaci(i,j)
      rjacc=rjacs*rjaci(i,j)
c
      dx=alpha*xxx-2.d0*beta*xxe+gamma*xee
      dy=alpha*yxx-2.d0*beta*yxe+gamma*yee
c
c     *****  coefficient of laplacian  *****
c
      c1(i,j) = alpha*rjacs
      c2(i,j) =-2.d0*beta*rjacs
      c3(i,j) = gamma*rjacs
      c4(i,j) = (xe(i,j)*dy-ye(i,j)*dx)*rjacc
      c5(i,j) = (yx(i,j)*dx-xx(i,j)*dy)*rjacc
c
  110 continue
c
c     *** metric on the surface wall ***
c
      do 120 i=2,nx
c
         ip1=i+1
c
      if(i.eq.nx) ip1=2
c
        xx(i,1 )=.5d0*(x(ip1,1)-x(i-1,1))
        yx(i,1 )=.5d0*(y(ip1,1)-y(i-1,1))
        xe(i,1 )=.5d0*(-x(i,3   )+4.d0*x(i,2   )-3.d0*x(i,1 ))
        ye(i,1 )=.5d0*(-y(i,3   )+4.d0*y(i,2   )-3.d0*y(i,1 ))
c
        rjac=xx(i,1)*ye(i,1)-xe(i,1)*yx(i,1)
        rjaci(i,1)=1.d0/rjac
c
  120 continue
c
c
      return
      end
c
c     *************************************
c     ** substituting boundary condition **
c     *************************************
c
      subroutine boundv
c
      implicit double precision(a-h,o-z)
      parameter(mx=401,my=201)
      common /divide/ nxp1,nyp1,nx,ny,nxm1,nym1
      common /result/ u(0:mx+1,0:my+1),v(0:mx+1,0:my+1),p(mx,my)
c
c
c     **** uniform flow condition ****
c
      do 100 i=1,nxp1
c
         u(i,nyp1)= u(i,nym1)
         v(i,nyp1)=-v(i,nym1)
         u(i,ny)  =1.d0
         v(i,ny)  =0.d0
c
c     *** on the body surface ***
c
         u(i,1) = 0.d0
         v(i,1) = 0.d0
         u(i,0) =-u(i,2)
         v(i,0) =-v(i,2)
c
  100 continue
c
c     *** branch cut ***
c
      do 120 j=1,ny
c
         u(0,j)    = u(nxm1,j)
         v(0,j)    = v(nxm1,j)
         u(1,j)    = u(nx,j)
         v(1,j)    = v(nx,j)
         u(nxp1,j) = u(2,j)
         v(nxp1,j) = v(2,j)
         u(nx+2,j) = u(3,j)
  120    v(nx+2,j) = v(3,j)
c
c
      return
      end
c
c     *************************************
c     ** substituting boundary condition **
c     *************************************
c
      subroutine boundp
c
      implicit double precision(a-h,o-z)
      parameter(mx=401,my=201)
      common /divide/ nxp1,nyp1,nx,ny,nxm1,nym1
      common /result/ u(0:mx+1,0:my+1),v(0:mx+1,0:my+1),p(mx,my)
      common /trans/ xx(mx,my),xe(mx,my),yx(mx,my),ye(mx,my)
      common /coeff/ c1(mx,my),c2(mx,my),c3(mx,my),c4(mx,my),c5(mx,my)
      common /reyno/ re,rei,dt
c
c     *** uniform flow ***
c     *** on the body surface ***
c
      do 100 i=2,nx
c
      vx=c1(i,2)*(u(i+1,2)-2.d0*u(i,2)+u(i-1,2))
     &  +c2(i,2)*(u(i+1,3)-u(i-1,3)-u(i+1,1)+u(i-1,1))*.25d0
     &  +c3(i,2)*(u(i,3)-2.d0*u(i,2)+u(i,1))
     & +(c4(i,2)*(u(i+1,2)-u(i-1,2))
     &  +c5(i,2)*(u(i  ,3)-u(i  ,1)))*.5d0
c
      vy=c1(i,2)*(v(i+1,2)-2.d0*v(i,2)+v(i-1,2))
     &  +c2(i,2)*(v(i+1,3)-v(i-1,3)-v(i+1,1)+v(i-1,1))*.25d0
     &  +c3(i,2)*(v(i,3)-2.d0*v(i,2)+v(i,1))
     & +(c4(i,2)*(v(i+1,2)-v(i-1,2))
     &  +c5(i,2)*(v(i  ,3)-v(i  ,1)))*.5d0
c
       p(i,1)=((4.d0*p(i,2)-p(i,3))
     &      -2.d0*(xe(i,1)*vx+ye(i,1)*vy)*rei)/3.d0
c
  100  p(i,ny)=(4.d0*p(i,nym1)-p(i,ny-2))/3.d0
c
c
      do 110 j=1,ny
c
         p(1   ,j)=p(nx,j)
         p(nxp1,j)=p(2 ,j)
c
  110 continue

c
c
      return
      end
c
c     **********************************
c     ** poisson equation of pressure **
c     **********************************
c
      subroutine poiss(ll)
c
      implicit double precision(a-h,o-z)
      parameter(mx=401,my=201)
      common /conv1/ eps1,omega1,iter1
      common /divide/ nxp1,nyp1,nx,ny,nxm1,nym1
      common /reyno/ re,rei,dt
      common /result/ u(0:mx+1,0:my+1),v(0:mx+1,0:my+1),p(mx,my)
      common /coeff/ c1(mx,my),c2(mx,my),c3(mx,my),c4(mx,my),c5(mx,my)
      common /trans/ xx(mx,my),xe(mx,my),yx(mx,my),ye(mx,my)
      common /jacobi/ rjaci(mx,my)
      dimension rhs(mx,my)
c
c     ***** right hand side of poisson equation *****
c
      do 100 j=2,nym1
c
         jp1=j+1
         jm1=j-1
c
      do 100 i=2,nx
c
         ip1=i+1
         im1=i-1
c
         ux=.5d0*(u(i+1,j)-u(i-1,j))
         ue=.5d0*(u(i,jp1)-u(i,jm1))
         vx=.5d0*(v(i+1,j)-v(i-1,j))
         ve=.5d0*(v(i,jp1)-v(i,jm1))
c
         dudx=ux*ye(i,j)-ue*yx(i,j)
         dvdy=ve*xx(i,j)-vx*xe(i,j)
         dudy=ue*xx(i,j)-ux*xe(i,j)
         dvdx=vx*ye(i,j)-ve*yx(i,j)
c
        dd=dudx**2+dvdy**2+2.d0*dudy*dvdx
       dcon=dd*rjaci(i,j)*rjaci(i,j)
c
c      *** cal. divergence ***
c
         div=(dudx+dvdy)*rjaci(i,j)
c
        rhs(i,j)=dcon-div/dt
c
  100 continue
c
c     ***** sor iteration loop *****
c
      do 500 loop=1,iter1
c
         err=0.d0
c
c     **** cal. internal points ****
c
      do 120 j=2,nym1
c
         jp1=j+1
         jm1=j-1
c
      do 120 i=2,nx
c
         ip1=i+1
         im1=i-1
c
c     ***** finite difference equation for poisson eq. *****
c
      dp=c1(i,j)*(p(ip1,j)+p(im1,j))
     &  +c2(i,j)*(p(ip1,jp1)-p(im1,jp1)-p(ip1,jm1)+p(im1,jm1))*.25d0
     &  +c3(i,j)*(p(i,jp1)+p(i,jm1))
     & +(c4(i,j)*(p(ip1,j)-p(im1,j))
     &  +c5(i,j)*(p(i,jp1)-p(i,jm1)))*.5d0+rhs(i,j)
c
c     ***** residual *****
c
        rp=dp/(2.d0*(c1(i,j)+c3(i,j)))-p(i,j)
c
        p(i,j)=p(i,j)+omega1*rp
c
        err=err+rp*rp
c
  120 continue
c
         call boundp
c
         err=dsqrt(err/dfloat((nxm1)*(ny-2)))
         if(err.lt.eps1) goto 99
c
  500 continue
c
   99 continue
c
      if(ll.eq.0) then
      write(*,610) loop,err
  610 format(' poisson iteration & error =',i5,3x,e10.4)
      end if
c
c
      return
      end
c
c     ***********************************
c     ** solved navier-stokes equation **
c     ***********************************
c
      subroutine navi(ll)
c
      implicit double precision(a-h,o-z)
      parameter(mx=401,my=201)
      common /divide/ nxp1,nyp1,nx,ny,nxm1,nym1
      common /reyno/ re,rei,dt
      common /conv2/ eps2,omega2,iter2
      common /result/ u(0:mx+1,0:my+1),v(0:mx+1,0:my+1),p(mx,my)
      common /coeff/ c1(mx,my),c2(mx,my),c3(mx,my),c4(mx,my),c5(mx,my)
      common /trans/ xx(mx,my),xe(mx,my),yx(mx,my),ye(mx,my)
      common /jacobi/ rjaci(mx,my)
      common /uold/ uo(mx,my),vo(mx,my)
      common /mean/ au(mx,my),av(mx,my),ap(mx,my)
      dimension coef(mx,my),dmyu(mx),dmyv(mx)
c
c     ***** coefficient for sor method *****
c
      do 100 j=2,nym1
      do 100 i=2,nx
c
         uo(i,j)=u(i,j)
         vo(i,j)=v(i,j)
c
         vis =2.d0*rei*(c1(i,j)+c3(i,j))
c
         coef(i,j)=1.d0/(1.d0+dt*vis)
c
  100 continue
c
c     ***** sor iteration loop *****
c
      do 150 loop=1,iter2
c
         err=0.d0
c
         call boundv
c
      do 110 j=2,nym1
c
         jp1=j+1
         jp2=j+2
         jm1=j-1
         jm2=j-2
c
      do 120 i=2,nx
c
         ip1=i+1
         ip2=i+2
         im1=i-1
         im2=i-2
c
c     *****  contravariant velocity component  *****
c
      cu=(uo(i,j)*ye(i,j)-vo(i,j)*xe(i,j))*rjaci(i,j)
      cv=(vo(i,j)*xx(i,j)-uo(i,j)*yx(i,j))*rjaci(i,j)
c
c     *****  convective term  (kawamura scheme)  *****
c
      u1=.5d0*(u(ip1,j)-u(im1,j))
      u3=.5d0*(u(i,jp1)-u(i,jm1))
c
      v1=.5d0*(v(ip1,j)-v(im1,j))
      v3=.5d0*(v(i,jp1)-v(i,jm1))
c
      cx=cu*u1+cv*u3
      cy=cu*v1+cv*v3
c
c     *****  pressure term  *****
c
      px=(ye(i,j)*(p(ip1,j)-p(im1,j))-yx(i,j)*(p(i,jp1)-p(i,jm1)))*.5d0
      py=(xx(i,j)*(p(i,jp1)-p(i,jm1))-xe(i,j)*(p(ip1,j)-p(im1,j)))*.5d0
c
      cpx=-cx-px*rjaci(i,j)
      cpy=-cy-py*rjaci(i,j)
c
c     *****  viscous term  *****
c
      vx=c1(i,j)*(u(ip1,j)+u(im1,j))
     &  +c2(i,j)*(u(ip1,jp1)-u(im1,jp1)-u(ip1,jm1)+u(im1,jm1))*.25d0
     &  +c3(i,j)*(u(i,jp1)+u(i,jm1))
     &  +c4(i,j)*(u(ip1,j)-u(im1,j))*.5d0
     &  +c5(i,j)*(u(i,jp1)-u(i,jm1))*.5d0
c
      vy=c1(i,j)*(v(ip1,j)+v(im1,j))
     &  +c2(i,j)*(v(ip1,jp1)-v(im1,jp1)-v(ip1,jm1)+v(im1,jm1))*.25d0
     &  +c3(i,j)*(v(i,jp1)+v(i,jm1))
     &  +c4(i,j)*(v(ip1,j)-v(im1,j))*.5d0
     &  +c5(i,j)*(v(i,jp1)-v(i,jm1))*.5d0
c
c
         drx=cpx+vx*rei
         dry=cpy+vy*rei
c
         du=coef(i,j)*(uo(i,j)+dt*drx)-u(i,j)
         dv=coef(i,j)*(vo(i,j)+dt*dry)-v(i,j)
c
        err=err+du*du+dv*dv
c
         u(i,j)=u(i,j)+omega2*du
         v(i,j)=v(i,j)+omega2*dv
c
  120 continue
c
  110 continue
c
         err=dsqrt(err/nxm1/(ny-2))
      if(err.lt.eps2) goto 99
c
  150 continue
c
   99 continue
c
      if(ll.eq.0) then
      write(*,610) loop,err
  610 format('navier-stokes loop & error =',i5,3x,e10.4)
      end if
c
c     *** for time average ***
c
      do 200 j=1,ny
      do 200 i=1,nx
c
         au(i,j)=au(i,j)+u(i,j)
         av(i,j)=av(i,j)+v(i,j)
         ap(i,j)=ap(i,j)+p(i,j)
c
  200 continue
c
c
      return
      end
c
c     **********************************
c     ** cal. drag & lift coefficient **
c     **********************************
c
      subroutine cdcl(ll)
c
      implicit double precision(a-h,o-z)
      parameter(mx=401,my=201)
      common /divide/ nxp1,nyp1,nx,ny,nxm1,nym1
      common /result/ u(0:mx+1,0:my+1),v(0:mx+1,0:my+1),p(mx,my)
      common /trans/ xx(mx,my),xe(mx,my),yx(mx,my),ye(mx,my)
      common /jacobi/ rjaci(mx,my)
      common /reyno/ re,rei,dt
c
         cd =0.d0
         cl =0.d0
c
      do 110 i=2,nx
c
        dx=xx(i,1)
        dy=yx(i,1)
c
        dp=p(i,1)
c
        cd=cd-2.d0*dp*dy
        cl=cl+2.d0*dp*dx
c
  110 continue
c
      if(ll.eq.0) then
      write(*,*)'********************************'
      write(*,600) cd,cl
      write(*,*)'********************************'
      write(*,*)'     '
      end if
c
  600 format(' cd & cl =',2f10.4)
c
c
      return
      end