c 2d upwind finite volume solver (2nd order in space and 1st order in time) implicit real*8 (a-h,o-z) parameter(nx=201,ny=201) common /aray1/ q(0:nx,0:ny,4) common /aray2/ d(0:nx,0:ny,4) common /aray3/ qlx(0:nx,0:ny,4),qrx(0:nx,0:ny,4), & qle(0:nx,0:ny,4),qre(0:nx,0:ny,4) common /aray4/ fl(nx,ny,4) common /aray5/ rc(0:nx,0:ny),uc(0:nx,0:ny),vc(0:nx,0:ny), & pc(0:nx,0:ny),cc(0:nx,0:ny) common /coor1/ x(nx,ny),y(nx,ny) common /coor2/ imax,jmax common /metr1/ axx(nx,ny),axy(nx,ny),axs(nx,ny) common /metr2/ aex(nx,ny),aey(nx,ny),aes(nx,ny) common /metr3/ vol(nx,ny) common /cntl0/ istart common /cntl1/ nstep,nout,mdd common /cntl2/ cflmx common /cntl3/ dt common /cntl4/ ibl,ibr,ibb,ibt common /cntl5/ nac common /parm1/ xminf common /parm2/ rinf,uinf,vinf,pinf,cinf common /parm3/ gamm,gam1 common /parm4/ ramp common /stat1/ time common /stat2/ loop,locl common /file1/ file_field,file_tecplt character*50 file_field,file_tecplt c initialization sequence call parm ! read run-time parameters call grid ! make a grid call metr ! metrices for cell interfaces call disp ! shows chosen run-time parameters call init ! set initial conditions if(istart.ne.0) then call dred ! read field data (if istart .ne. 0) endif c main sequence locl =0 1000 continue locl =locl+1 loop =loop+1 call rbdy ! set boundary conditions for MUSCL step if(nac.eq.1) then call slp_1st ! set interface values (1st order) else call slp_2nd ! set interface values (2nd order) endif call sbdy ! set boundary conditions for flux calculation call flux_xi ! numerical flux for xi-interfaces call flux_eta ! numerical flux for eta-interface call dtsz ! determine time step size call intg ! integrate equations c post sequence m =mod(locl,nout) if(m.eq.0) then call wrtd call tecp endif if(locl.lt.nstep) go to 1000 c exit stop end subroutine parm c ************************************************************ c * run_time parameters * c ************************************************************ implicit real*8 (a-h,o-z) parameter(nx=201,ny=201) common /cntl0/ istart common /cntl1/ nstep,nout,mdd common /cntl2/ cflmx common /cntl3/ dt common /cntl4/ ibl,ibr,ibb,ibt common /cntl5/ nac common /parm1/ xminf common /parm3/ gamm,gam1 common /parm4/ ramp common /file1/ file_field,file_tecplt character*50 file_field,file_tecplt c set run-time parameters write(6,*) '?input : file name of field data ...' read(5,500) file_field write(6,*) '?input : file name of tecplot data ...' read(5,500) file_tecplt 500 format(a50) write(6,*) '?select : start from initial(0) / disk(1) ...' read(5,*) istart write(6,*) '?input : number of steps (nstep) ...' read(5,*) nstep write(6,*) '?input : interval to store data (nout) ...' read(5,*) nout write(6,*) '?input : mdd (interval for residual est.) ...' read(5,*) mdd write(6,*) '?input : nac (1=1st order, 2=2nd order) ...' read(5,*) nac write(6,*) '?input : maximum allowable cfl limit ...' read(5,*) cflmx write(6,*) '?input : mach (freestream Mach number) ...' read(5,*) xminf write(6,*) '?input : gamm (Cp/Cv) ...' read(5,*) gamm gam1 =gamm-1.0d0 write(6,*) '?choose : left BC (1=sld, 2=amb, 3=otf, 4=ext)' read(5,*) ibl write(6,*) '?choose : right BC (1=sld, 2=amb, 3=otf, 4=ext)' read(5,*) ibr write(6,*) '?choose : bottom BC (1=sld, 2=amb, 3=otf, 4=ext)' read(5,*) ibb write(6,*) '?choose : top BC (1=sld, 2=amb, 3=otf, 4=ext)' read(5,*) ibt write(6,*) '?input : ramp angle (deg) ...' read(5,*) ramp c exit return end subroutine grid c ************************************************************ c * make a grid system * c ************************************************************ implicit real*8 (a-h,o-z) parameter(nx=201,ny=201) common /coor1/ x(nx,ny),y(nx,ny) common /coor2/ imax,jmax common /parm4/ ramp c parameters imax =101 ! number of grid points in xi-direction jmax =101 ! number of grid points in eta-direction xl = 1.0d0 ! horizontal length of the computational domain yl = 1.0d0 ! vertical length of the computational domain c initial mesh dx =xl/dfloat(imax-1) dy =yl/dfloat(jmax-1) do j=1,jmax do i=1,imax x(i,j) =dx*dfloat(i-1) y(i,j) =dy*dfloat(j-1) end do end do c make a ramp ib =(1+imax)/4 pai =dacos(-1.0d0) the =ramp*pai/180.0d0 do i=ib,imax y(i,1) =y(i,1)+dtan(the)*(x(i,1)-x(ib,1)) end do c regeneration do i=ib,imax yt =y(i,jmax)-y(i,1) dy =yt/dfloat(jmax-1) do j=1,jmax y(i,j) =y(i,1)+dy*float(j-1) end do end do c termination return end subroutine metr c ************************************************************ c * evaluation of metrices * c ************************************************************ implicit real*8 (a-h,o-z) parameter(nx=201,ny=201) common /coor1/ x(nx,ny),y(nx,ny) common /coor2/ imax,jmax common /metr1/ axx(nx,ny),axy(nx,ny),axs(nx,ny) common /metr2/ aex(nx,ny),aey(nx,ny),aes(nx,ny) common /metr3/ vol(nx,ny) c axx,axy & axs (xi-const interfaces) do j=1,jmax-1 do i=1,imax axxt = y( i,j+1)-y( i, j) axyt =-(x( i,j+1)-x( i, j)) axs(i,j) =dsqrt(axxt*axxt+axyt*axyt) axx(i,j) =axxt/axs(i,j) axy(i,j) =axyt/axs(i,j) end do end do c aex,aey & aes (eta-const interfaces) do j=1,jmax do i=1,imax-1 aext =-(y(i+1, j)-y( i, j)) aeyt = x(i+1, j)-x( i, j) aes(i,j) =dsqrt(aext*aext+aeyt*aeyt) aex(i,j) =aext/aes(i,j) aey(i,j) =aeyt/aes(i,j) end do end do c volume do j=1,jmax-1 do i=1,imax-1 dx13 =x(i+1,j+1)-x(i,j) dy13 =y(i+1,j+1)-y(i,j) vol(i,j) =0.5d0*(dx13*(axx(i,j)*axs(i,j)+aex(i,j)*aes(i,j)) & +dy13*(axy(i,j)*axs(i,j)+aey(i,j)*aes(i,j))) end do end do c termination return end subroutine disp c ************************************************************ c * displays current run_time parameters * c ************************************************************ implicit real*8 (a-h,o-z) parameter(nx=201,ny=201) common /coor2/ imax,jmax common /cntl0/ istart common /cntl1/ nstep,nout,mdd common /cntl2/ cflmx common /cntl3/ dt common /cntl5/ nac common /parm1/ xminf common /parm3/ gamm,gam1 common /parm4/ ramp common /stat1/ time common /stat2/ loop,locl if(istart.eq.0) then write(6,600) 600 format(///1h ,'disp) calculation starts from initial state.') else write(6,610) loop,time 610 format(///1h ,'disp) calculation starts from the disk file.'/ & 1h ,' loop = ',i8,' time = ',e12.5) endif write(6,620) imax,jmax 620 format(//1h ,'disp) coordinate parameters.'/ & 1h ,' imax = ',i8,' jmax = ',i8) if(nac.eq.1) then write(6,625) 625 format(//1h ,'disp) first order in space ') else write(6,626) 626 format(//1h ,'disp) second order in space ') endif write(6,630) cflmx 630 format(//1h ,'disp) step size control.'/ & 1h ,' cflmx = ',e12.5) write(6,640) nstep,nout,mdd 640 format(//1h ,'disp) run_time parameters(1).'/ & 1h ,' nstep = ',i8/ & 1h ,' nout = ',i8/ & 1h ,' mdd = ',i8) write(6,650) xminf,gamm,ramp 650 format(//1h ,'disp) run_time parameters(2).'/ & 1h ,' xminf = ',e12.5/ & 1h ,' gamm = ',e12.5/ & 1h ,' ramp = ',e12.5//) c termination return end subroutine init c ************************************************************ c * set initial condition * c ************************************************************ implicit real*8 (a-h,o-z) parameter(nx=201,ny=201) common /aray1/ q(0:nx,0:ny,4) common /coor2/ imax,jmax common /parm1/ xminf common /parm2/ rinf,uinf,vinf,pinf,cinf common /parm3/ gamm,gam1 common /stat1/ time common /stat2/ loop,locl c characterization of freestream gas rinf =1.0d0 uinf =xminf vinf =0.0d0 pinf =1.0d0/gamm cinf =1.0d0 einf =pinf/gam1+0.5d0*rinf*xminf*xminf ruinf =rinf*uinf rvinf =rinf*vinf c set initial condition do j=1,jmax-1 do i=1,imax-1 q(i,j,1) =rinf q(i,j,2) =ruinf q(i,j,3) =rvinf q(i,j,4) =einf end do end do c set time parameters loop =0 time =0.0 c termination return end subroutine dred c ************************************************************ c * read disk data (if istart .ne. 0) * c ************************************************************ implicit real*8 (a-h,o-z) parameter(nx=201,ny=201) common /aray1/ q(0:nx,0:ny,4) common /coor2/ imax,jmax common /stat1/ time common /stat2/ loop,locl common /file1/ file_field,file_tecplt character*50 file_field,file_tecplt dimension icom(20),pcom(20) c read data open(15,file=file_field,form='formatted') read(15,*) (icom(i),i=1,20) read(15,*) (pcom(i),i=1,20) read(15,*) ((q(i,j,1),q(i,j,2),q(i,j,3),q(i,j,4), & i=1,imax),j=1,jmax) close(15) c set parameters time =pcom(1) loop =icom(3) write(6,600) loop,time 600 format(/1x,'dred:disp) loop/time are reset to followings.'/ & 1x,' loop = ',i13/ & 1x,' time = ',e13.6) c termination return end subroutine dtsz c ************************************************************ c * step_size estimation * c ************************************************************ implicit real*8 (a-h,o-z) parameter(nx=201,ny=201) common /aray1/ q(0:nx,0:ny,4) common /aray5/ rc(0:nx,0:ny),uc(0:nx,0:ny),vc(0:nx,0:ny), & pc(0:nx,0:ny),cc(0:nx,0:ny) common /coor2/ imax,jmax common /metr1/ axx(nx,ny),axy(nx,ny),axs(nx,ny) common /metr2/ aex(nx,ny),aey(nx,ny),aes(nx,ny) common /metr3/ vol(nx,ny) common /cntl2/ cflmx common /cntl3/ dt common /parm3/ gamm,gam1 c step size estimation dt =1.0e30 do j=1,jmax-1 do i=1,imax-1 cvxmp =dmax1(0.5d0*(cc(i-1,j)-uc(i-1,j)*axx(i ,j) & -vc(i-1,j)*axy(i ,j) & +cc(i ,j)-uc(i ,j)*axx(i ,j) & -vc(i ,j)*axy(i ,j)), & 0.0d0) cvxpp =dmax1(0.5d0*(cc(i+1,j)+uc(i+1,j)*axx(i+1,j) & +vc(i+1,j)*axy(i+1,j) & +cc(i ,j)+uc(i ,j)*axx(i+1,j) & +vc(i ,j)*axy(i+1,j)), & 0.0d0) cvemp =dmax1(0.5d0*(cc(i,j-1)-uc(i,j-1)*aex(i,j ) & -vc(i,j-1)*aey(i,j ) & +cc(i,j )-uc(i,j )*aex(i,j ) & -vc(i,j )*aey(i,j )), & 0.0d0) cvepp =dmax1(0.5d0*(cc(i,j+1)+uc(i,j+1)*aex(i,j+1) & +vc(i,j+1)*aey(i,j+1) & +cc(i,j )+uc(i,j )*aex(i,j+1) & +vc(i,j )*aey(i,j+1)), & 0.0d0) dt =vol(i,j)/(cvxmp*axs(i,j)+cvxpp*axs(i+1,j) & +cvemp*aes(i,j)+cvepp*aes(i,j+1)) end do end do dt =cflmx*dt c termination return end subroutine rbdy c ************************************************************ c * set ghost cells * c ************************************************************ implicit real*8 (a-h,o-z) parameter(nx=201,ny=201) common /aray1/ q(0:nx,0:ny,4) common /aray5/ rc(0:nx,0:ny),uc(0:nx,0:ny),vc(0:nx,0:ny), & pc(0:nx,0:ny),cc(0:nx,0:ny) common /coor2/ imax,jmax common /metr1/ axx(nx,ny),axy(nx,ny),axs(nx,ny) common /metr2/ aex(nx,ny),aey(nx,ny),aes(nx,ny) common /cntl4/ ibl,ibr,ibb,ibt common /parm2/ rinf,uinf,vinf,pinf,cinf common /parm3/ gamm,gam1 c convert conservative variables into primitive variables do j=1,jmax-1 do i=1,imax-1 rc(i,j) =q(i,j,1) uc(i,j) =q(i,j,2)/q(i,j,1) vc(i,j) =q(i,j,3)/q(i,j,1) pc(i,j) =gam1*(q(i,j,4)-0.5d0*(uc(i,j)*q(i,j,2) & +vc(i,j)*q(i,j,3))) cc(i,j) =dsqrt(gamm*pc(i,j)/rc(i,j)) end do end do c left boundary (xi) if(ibl.eq.1) then do j=1,jmax-1 qn =2.0d0*(axx(1,j)*uc(1,j)+axy(1,j)*vc(1,j)) rc(0,j) =rc(1,j) uc(0,j) =uc(1,j)-qn*axx(1,j) vc(0,j) =vc(1,j)-qn*axy(1,j) pc(0,j) =pc(1,j) cc(0,j) =cc(1,j) end do elseif(ibl.eq.2) then do j=1,jmax-1 rc(0,j) =rinf uc(0,j) =uinf vc(0,j) =vinf pc(0,j) =pinf cc(0,j) =cinf end do elseif(ibl.eq.3) then do j=1,jmax-1 rc(0,j) =rc(1,j) uc(0,j) =uc(1,j) vc(0,j) =vc(1,j) pc(0,j) =pinf cc(0,j) =dsqrt(gamm*pinf/rc(0,j)) end do else do j=1,jmax-1 rc(0,j) =rc(1,j) uc(0,j) =uc(1,j) vc(0,j) =vc(1,j) pc(0,j) =pc(1,j) cc(0,j) =cc(1,j) end do endif c right boundary (xi) if(ibr.eq.1) then do j=1,jmax-1 qn =2.0d0*(axx(imax,j)*uc(imax-1,j) & +axy(imax,j)*vc(imax-1,j)) rc(imax,j) =rc(imax-1,j) uc(imax,j) =uc(imax-1,j)-qn*axx(imax,j) vc(imax,j) =vc(imax-1,j)-qn*axy(imax,j) pc(imax,j) =pc(imax-1,j) cc(imax,j) =cc(imax-1,j) end do elseif(ibr.eq.2) then do j=1,jmax-1 rc(imax,j) =rinf uc(imax,j) =uinf vc(imax,j) =vinf pc(imax,j) =pinf cc(imax,j) =cinf end do elseif(ibr.eq.3) then do j=1,jmax-1 rc(imax,j) =rc(imax-1,j) uc(imax,j) =uc(imax-1,j) vc(imax,j) =vc(imax-1,j) pc(imax,j) =pinf cc(imax,j) =dsqrt(gamm*pinf/rc(imax,j)) end do else do j=1,jmax-1 rc(imax,j) =rc(imax-1,j) uc(imax,j) =uc(imax-1,j) vc(imax,j) =vc(imax-1,j) pc(imax,j) =pc(imax-1,j) cc(imax,j) =cc(imax-1,j) end do endif c bottom boundary (eta) if(ibb.eq.1) then do i=1,imax-1 qn =2.0d0*(aex(i,1)*uc(i,1)+aey(i,1)*vc(i,1)) rc(i,0) =rc(i,1) uc(i,0) =uc(i,1)-qn*aex(i,1) vc(i,0) =vc(i,1)-qn*aey(i,1) pc(i,0) =pc(i,1) cc(i,0) =cc(i,1) end do elseif(ibb.eq.2) then do i=1,imax-1 rc(i,0) =rinf uc(i,0) =uinf vc(i,0) =vinf pc(i,0) =pinf cc(i,0) =cinf end do elseif(ibb.eq.3) then do i=1,imax-1 rc(i,0) =rc(i,1) uc(i,0) =uc(i,1) vc(i,0) =vc(i,1) pc(i,0) =pinf cc(i,0) =dsqrt(gamm*pinf/rc(i,0)) end do else do i=1,imax-1 rc(i,0) =rc(i,1) uc(i,0) =uc(i,1) vc(i,0) =vc(i,1) pc(i,0) =pc(i,1) cc(i,0) =cc(i,1) end do endif c top boundary (eta) if(ibt.eq.1) then do i=1,imax-1 qn =2.0d0*(aex(i,jmax)*uc(i,jmax-1) & +aey(i,jmax)*vc(i,jmax-1)) rc(i,jmax) =rc(i,jmax-1) uc(i,jmax) =uc(i,jmax-1)-qn*aex(i,jmax) vc(i,jmax) =vc(i,jmax-1)-qn*aey(i,jmax) pc(i,jmax) =pc(i,jmax-1) cc(i,jmax) =cc(i,jmax-1) end do elseif(ibt.eq.2) then do i=1,imax-1 rc(i,jmax) =rinf uc(i,jmax) =uinf vc(i,jmax) =vinf pc(i,jmax) =pinf cc(i,jmax) =cinf end do elseif(ibt.eq.3) then do i=1,imax-1 rc(i,jmax) =rc(i,jmax-1) uc(i,jmax) =uc(i,jmax-1) vc(i,jmax) =vc(i,jmax-1) pc(i,jmax) =pinf cc(i,jmax) =dsqrt(gamm*pinf/rc(i,jmax)) end do else do i=1,imax-1 rc(i,jmax) =rc(i,jmax-1) uc(i,jmax) =uc(i,jmax-1) vc(i,jmax) =vc(i,jmax-1) pc(i,jmax) =pc(i,jmax-1) cc(i,jmax) =cc(i,jmax-1) end do endif c termination return end subroutine slp_1st c ************************************************************ c * set interface values (1st order in space) * c ************************************************************ implicit real*8 (a-h,o-z) parameter(nx=201,ny=201) common /aray1/ q(0:nx,0:ny,4) common /aray3/ qlx(0:nx,0:ny,4),qrx(0:nx,0:ny,4), & qle(0:nx,0:ny,4),qre(0:nx,0:ny,4) common /aray5/ rc(0:nx,0:ny),uc(0:nx,0:ny),vc(0:nx,0:ny), & pc(0:nx,0:ny),cc(0:nx,0:ny) common /coor2/ imax,jmax c xi do j=1,jmax-1 do i=1,imax-1 qlx(i+1,j,1) =rc(i,j) qlx(i+1,j,2) =uc(i,j) qlx(i+1,j,3) =vc(i,j) qlx(i+1,j,4) =pc(i,j) qrx(i ,j,1) =rc(i,j) qrx(i ,j,2) =uc(i,j) qrx(i ,j,3) =vc(i,j) qrx(i ,j,4) =pc(i,j) end do end do do j=1,jmax-1 qlx(1 ,j,1) =rc(0 ,j) qlx(1 ,j,2) =uc(0 ,j) qlx(1 ,j,3) =vc(0 ,j) qlx(1 ,j,4) =pc(0 ,j) qrx(imax,j,1) =rc(imax,j) qrx(imax,j,2) =uc(imax,j) qrx(imax,j,3) =vc(imax,j) qrx(imax,j,4) =pc(imax,j) end do c eta do j=1,jmax-1 do i=1,imax-1 qle(i,j+1,1) =rc(i,j) qle(i,j+1,2) =uc(i,j) qle(i,j+1,3) =vc(i,j) qle(i,j+1,4) =pc(i,j) qre(i,j ,1) =rc(i,j) qre(i,j ,2) =uc(i,j) qre(i,j ,3) =vc(i,j) qre(i,j ,4) =pc(i,j) end do end do do i=1,imax-1 qle(i,1 ,1) =rc(i,0 ) qle(i,1 ,2) =uc(i,0 ) qle(i,1 ,3) =vc(i,0 ) qle(i,1 ,4) =pc(i,0 ) qre(i,jmax,1) =rc(i,jmax) qre(i,jmax,2) =uc(i,jmax) qre(i,jmax,3) =vc(i,jmax) qre(i,jmax,4) =pc(i,jmax) end do c termination return end subroutine slp_2nd c ************************************************************ c * set interface values (2nd order in space) * c ************************************************************ implicit real*8 (a-h,o-z) parameter(nx=201,ny=201) common /aray1/ q(0:nx,0:ny,4) common /aray3/ qlx(0:nx,0:ny,4),qrx(0:nx,0:ny,4), & qle(0:nx,0:ny,4),qre(0:nx,0:ny,4) common /aray5/ rc(0:nx,0:ny),uc(0:nx,0:ny),vc(0:nx,0:ny), & pc(0:nx,0:ny),cc(0:nx,0:ny) common /coor1/ x(nx,ny),y(nx,ny) common /coor2/ imax,jmax common /metr1/ axx(nx,ny),axy(nx,ny),axs(nx,ny) common /metr2/ aex(nx,ny),aey(nx,ny),aes(nx,ny) common /metr3/ vol(nx,ny) common /parm3/ gamm,gam1 c estimation of gradient values do j=1,jmax-1 do i=1,imax-1 drip =rc(i+1,j )-rc(i ,j ) drim =rc(i ,j )-rc(i-1,j ) drjp =rc(i ,j+1)-rc(i ,j ) drjm =rc(i ,j )-rc(i ,j-1) duip =uc(i+1,j )-uc(i ,j ) duim =uc(i ,j )-uc(i-1,j ) dujp =uc(i ,j+1)-uc(i ,j ) dujm =uc(i ,j )-uc(i ,j-1) dvip =vc(i+1,j )-vc(i ,j ) dvim =vc(i ,j )-vc(i-1,j ) dvjp =vc(i ,j+1)-vc(i ,j ) dvjm =vc(i ,j )-vc(i ,j-1) dpip =pc(i+1,j )-pc(i ,j ) dpim =pc(i ,j )-pc(i-1,j ) dpjp =pc(i ,j+1)-pc(i ,j ) dpjm =pc(i ,j )-pc(i ,j-1) dri =(dsign(0.25d0,drip)+dsign(0.25d0,drim)) & *dmin1(dabs(drip), dabs(drim)) drj =(dsign(0.25d0,drjp)+dsign(0.25d0,drjm)) & *dmin1(dabs(drjp), dabs(drjm)) dui =(dsign(0.25d0,duip)+dsign(0.25d0,duim)) & *dmin1(dabs(duip), dabs(duim)) duj =(dsign(0.25d0,dujp)+dsign(0.25d0,dujm)) & *dmin1(dabs(dujp), dabs(dujm)) dvi =(dsign(0.25d0,dvip)+dsign(0.25d0,dvim)) & *dmin1(dabs(dvip), dabs(dvim)) dvj =(dsign(0.25d0,dvjp)+dsign(0.25d0,dvjm)) & *dmin1(dabs(dvjp), dabs(dvjm)) dpi =(dsign(0.25d0,dpip)+dsign(0.25d0,dpim)) & *dmin1(dabs(dpip), dabs(dpim)) dpj =(dsign(0.25d0,dpjp)+dsign(0.25d0,dpjm)) & *dmin1(dabs(dpjp), dabs(dpjm)) qlx(i+1,j,1) =rc(i,j)+dri qlx(i+1,j,2) =uc(i,j)+dui qlx(i+1,j,3) =vc(i,j)+dvi qlx(i+1,j,4) =pc(i,j)+dpi qrx(i ,j,1) =rc(i,j)-dri qrx(i ,j,2) =uc(i,j)-dui qrx(i ,j,3) =vc(i,j)-dvi qrx(i ,j,4) =pc(i,j)-dpi qle(i,j+1,1) =rc(i,j)+drj qle(i,j+1,2) =uc(i,j)+duj qle(i,j+1,3) =vc(i,j)+dvj qle(i,j+1,4) =pc(i,j)+dpj qre(i ,j,1) =rc(i,j)-drj qre(i ,j,2) =uc(i,j)-duj qre(i ,j,3) =vc(i,j)-dvj qre(i ,j,4) =pc(i,j)-dpj end do end do c termination return end subroutine sbdy c ************************************************************ c * set interface boundary values * c ************************************************************ implicit real*8 (a-h,o-z) parameter(nx=201,ny=201) common /aray1/ q(0:nx,0:ny,4) common /aray3/ qlx(0:nx,0:ny,4),qrx(0:nx,0:ny,4), & qle(0:nx,0:ny,4),qre(0:nx,0:ny,4) common /coor2/ imax,jmax common /metr1/ axx(nx,ny),axy(nx,ny),axs(nx,ny) common /metr2/ aex(nx,ny),aey(nx,ny),aes(nx,ny) common /cntl4/ ibl,ibr,ibb,ibt common /parm2/ rinf,uinf,vinf,pinf,cinf common /parm3/ gamm,gam1 c left boundary (xi) if(ibl.eq.1) then do j=1,jmax-1 qn =2.0d0*(axx(1,j)*qrx(1,j,2)+axy(1,j)*qrx(1,j,3)) qlx(1,j,1) =qrx(1,j,1) qlx(1,j,2) =qrx(1,j,2)-qn*axx(1,j) qlx(1,j,3) =qrx(1,j,3)-qn*axy(1,j) qlx(1,j,4) =qrx(1,j,4) end do elseif(ibl.eq.2) then do j=1,jmax-1 qlx(1,j,1) =rinf qlx(1,j,2) =uinf qlx(1,j,3) =vinf qlx(1,j,4) =pinf end do elseif(ibl.eq.3) then do j=1,jmax-1 qlx(1,j,1) =qrx(1,j,1) qlx(1,j,2) =qrx(1,j,2) qlx(1,j,3) =qrx(1,j,3) qlx(1,j,4) =pinf end do else do j=1,jmax-1 qlx(1,j,1) =qrx(1,j,1) qlx(1,j,2) =qrx(1,j,2) qlx(1,j,3) =qrx(1,j,3) qlx(1,j,4) =qrx(1,j,4) end do endif c right boundary (xi) if(ibr.eq.1) then do j=1,jmax-1 qn =2.0d0*(axx(imax,j)*qlx(imax,j,2) & +axy(imax,j)*qlx(imax,j,3)) qrx(imax,j,1) =qlx(imax,j,1) qrx(imax,j,2) =qlx(imax,j,2)-qn*axx(imax,j) qrx(imax,j,3) =qlx(imax,j,3)-qn*axy(imax,j) qrx(imax,j,4) =qlx(imax,j,4) end do elseif(ibr.eq.2) then do j=1,jmax-1 qrx(imax,j,1) =rinf qrx(imax,j,2) =uinf qrx(imax,j,3) =vinf qrx(imax,j,4) =pinf end do elseif(ibr.eq.3) then do j=1,jmax-1 qrx(imax,j,1) =qlx(imax,j,1) qrx(imax,j,2) =qlx(imax,j,2) qrx(imax,j,3) =qlx(imax,j,3) qrx(imax,j,4) =pinf end do else do j=1,jmax-1 qrx(imax,j,1) =qlx(imax,j,1) qrx(imax,j,2) =qlx(imax,j,2) qrx(imax,j,3) =qlx(imax,j,3) qrx(imax,j,4) =qlx(imax,j,4) end do endif c bottom boundary (eta) if(ibb.eq.1) then do i=1,imax-1 qn =2.0d0*(aex(i,1)*qre(i,1,2)+aey(i,1)*qre(i,1,3)) qle(i,1,1) =qre(i,1,1) qle(i,1,2) =qre(i,1,2)-qn*aex(i,1) qle(i,1,3) =qre(i,1,3)-qn*aey(i,1) qle(i,1,4) =qre(i,1,4) end do elseif(ibb.eq.2) then do i=1,imax-1 qle(i,1,1) =rinf qle(i,1,2) =uinf qle(i,1,3) =vinf qle(i,1,4) =pinf end do elseif(ibb.eq.3) then do i=1,imax-1 qle(i,1,1) =qre(i,1,1) qle(i,1,2) =qre(i,1,2) qle(i,1,3) =qre(i,1,3) qle(i,1,4) =pinf end do else do i=1,imax-1 qle(i,1,1) =qre(i,1,1) qle(i,1,2) =qre(i,1,2) qle(i,1,3) =qre(i,1,3) qle(i,1,4) =qre(i,1,4) end do endif c top boundary (eta) if(ibt.eq.1) then do i=1,imax-1 qn =2.0d0*(aex(i,jmax)*qle(i,jmax,2) & +aey(i,jmax)*qle(i,jmax,3)) qre(i,jmax,1) =qle(i,jmax,1) qre(i,jmax,2) =qle(i,jmax,2)-qn*aex(i,jmax) qre(i,jmax,3) =qle(i,jmax,3)-qn*aey(i,jmax) qre(i,jmax,4) =qle(i,jmax,4) end do elseif(ibt.eq.2) then do i=1,imax-1 qre(i,jmax,1) =rinf qre(i,jmax,2) =uinf qre(i,jmax,3) =vinf qre(i,jmax,4) =pinf end do elseif(ibt.eq.3) then do i=1,imax-1 qre(i,jmax,1) =qle(i,jmax,1) qre(i,jmax,2) =qle(i,jmax,2) qre(i,jmax,3) =qle(i,jmax,3) qre(i,jmax,4) =pinf end do else do i=1,imax-1 qre(i,jmax,1) =qle(i,jmax,1) qre(i,jmax,2) =qle(i,jmax,2) qre(i,jmax,3) =qle(i,jmax,3) qre(i,jmax,4) =qle(i,jmax,4) end do endif c termination return end subroutine flux_xi c ************************************************************ c * upwind flux for xi-const. interface * c ************************************************************ implicit real*8 (a-h,o-z) parameter(nx=201,ny=201) common /aray1/ q(0:nx,0:ny,4) common /aray2/ d(0:nx,0:ny,4) common /aray3/ qlx(0:nx,0:ny,4),qrx(0:nx,0:ny,4), & qle(0:nx,0:ny,4),qre(0:nx,0:ny,4) common /aray4/ fl(nx,ny,4) common /coor2/ imax,jmax common /metr1/ axx(nx,ny),axy(nx,ny),axs(nx,ny) common /parm3/ gamm,gam1 c numerical flux function for xi-const interfaces do j=1,jmax-1 do i=1,imax c left state rl =qlx(i,j,1) ul =qlx(i,j,2) vl =qlx(i,j,3) pl =qlx(i,j,4) unl =axx(i,j)*ul+axy(i,j)*vl cl2 =gamm*pl/rl cl =dsqrt(cl2) hl =cl2/gam1+0.5d0*(ul*ul+vl*vl) el =rl*hl-pl fl1 =rl*unl fl2 =fl1*ul+axx(i,j)*pl fl3 =fl1*vl+axy(i,j)*pl fl4 =fl1*hl c right state rr =qrx(i,j,1) ur =qrx(i,j,2) vr =qrx(i,j,3) pr =qrx(i,j,4) unr =axx(i,j)*ur+axy(i,j)*vr cr2 =gamm*pr/rr cr =dsqrt(cr2) hr =cr2/gam1+0.5d0*(ur*ur+vr*vr) er =rr*hr-pr fr1 =rr*unr fr2 =fr1*ur+axx(i,j)*pr fr3 =fr1*vr+axy(i,j)*pr fr4 =fr1*hr c dissipation coefficient dd =0.5d0*(dabs(unl)+cl+dabs(unr)+cr) c state differences drc =rr -rl dru =rr*ur-rl*ul drv =rr*vr-rl*vl dec =er -el c numerical flux fl(i,j,1) =0.5d0*(fr1+fl1-dd*drc)*axs(i,j) fl(i,j,2) =0.5d0*(fr2+fl2-dd*dru)*axs(i,j) fl(i,j,3) =0.5d0*(fr3+fl3-dd*drv)*axs(i,j) fl(i,j,4) =0.5d0*(fr4+fl4-dd*dec)*axs(i,j) end do end do c integration (xi contributing part) do j=1,jmax-1 do i=1,imax-1 d(i,j,1) =-(fl(i+1,j,1)-fl(i,j,1)) d(i,j,2) =-(fl(i+1,j,2)-fl(i,j,2)) d(i,j,3) =-(fl(i+1,j,3)-fl(i,j,3)) d(i,j,4) =-(fl(i+1,j,4)-fl(i,j,4)) end do end do c termination return end subroutine flux_eta c ************************************************************ c * upwind flux for eta-const. interface * c ************************************************************ implicit real*8 (a-h,o-z) parameter(nx=201,ny=201) common /aray1/ q(0:nx,0:ny,4) common /aray2/ d(0:nx,0:ny,4) common /aray3/ qlx(0:nx,0:ny,4),qrx(0:nx,0:ny,4), & qle(0:nx,0:ny,4),qre(0:nx,0:ny,4) common /aray4/ fl(nx,ny,4) common /coor2/ imax,jmax common /metr2/ aex(nx,ny),aey(nx,ny),aes(nx,ny) common /parm3/ gamm,gam1 c numerical flux function for xi-const interfaces do j=1,jmax do i=1,imax-1 c left state rl =qle(i,j,1) ul =qle(i,j,2) vl =qle(i,j,3) pl =qle(i,j,4) unl =aex(i,j)*ul+aey(i,j)*vl cl2 =gamm*pl/rl cl =dsqrt(cl2) hl =cl2/gam1+0.5d0*(ul*ul+vl*vl) el =rl*hl-pl fl1 =rl*unl fl2 =fl1*ul+aex(i,j)*pl fl3 =fl1*vl+aey(i,j)*pl fl4 =fl1*hl c right state rr =qre(i,j,1) ur =qre(i,j,2) vr =qre(i,j,3) pr =qre(i,j,4) unr =aex(i,j)*ur+aey(i,j)*vr cr2 =gamm*pr/rr cr =dsqrt(cr2) hr =cr2/gam1+0.5d0*(ur*ur+vr*vr) er =rr*hr-pr fr1 =rr*unr fr2 =fr1*ur+aex(i,j)*pr fr3 =fr1*vr+aey(i,j)*pr fr4 =fr1*hr c dissipation coefficient dd =0.5d0*(dabs(unl)+cl+dabs(unr)+cr) c state differences drc =rr -rl dru =rr*ur-rl*ul drv =rr*vr-rl*vl dec =er -el c numerical flux fl(i,j,1) =0.5d0*(fr1+fl1-dd*drc)*aes(i,j) fl(i,j,2) =0.5d0*(fr2+fl2-dd*dru)*aes(i,j) fl(i,j,3) =0.5d0*(fr3+fl3-dd*drv)*aes(i,j) fl(i,j,4) =0.5d0*(fr4+fl4-dd*dec)*aes(i,j) end do end do c integration (eta contributing part) do j=1,jmax-1 do i=1,imax-1 d(i,j,1) =d(i,j,1)-(fl(i,j+1,1)-fl(i,j,1)) d(i,j,2) =d(i,j,2)-(fl(i,j+1,2)-fl(i,j,2)) d(i,j,3) =d(i,j,3)-(fl(i,j+1,3)-fl(i,j,3)) d(i,j,4) =d(i,j,4)-(fl(i,j+1,4)-fl(i,j,4)) end do end do c termination return end subroutine intg c ************************************************************ c * integration * c ************************************************************ implicit real*8 (a-h,o-z) parameter(nx=201,ny=201) common /aray1/ q(0:nx,0:ny,4) common /aray2/ d(0:nx,0:ny,4) common /coor2/ imax,jmax common /metr3/ vol(nx,ny) common /cntl1/ nstep,nout,mdd common /cntl2/ cflmx common /cntl3/ dt common /stat1/ time common /stat2/ loop,locl common /parm3/ gamm,gam1 c set new time time =time+dt c residual estimation iwr =mod(loop,mdd) if(iwr.eq.0) then drmx =-1.0d30 rms = 0.0d0 do j=1,jmax-1 do i=1,imax-1 drot =dabs(d(i,j,1)/vol(i,j)) rms =rms+drot*drot if(drot.gt.drmx) then drmx =drot irmx =i jrmx =j endif end do end do rms =dsqrt(rms/(dfloat(imax-1)*dfloat(jmax-1))) write(6,600) loop,time,dt,drmx,irmx,jrmx,rms 600 format(1x,i8,1x,e10.3,1x,e10.3,1x,e10.3,1x,2i8,1x,e10.3) endif c advance a time step do j=1,jmax-1 do i=1,imax-1 dt1 =dt/vol(i,j) q(i,j,1) =q(i,j,1)+dt1*d(i,j,1) q(i,j,2) =q(i,j,2)+dt1*d(i,j,2) q(i,j,3) =q(i,j,3)+dt1*d(i,j,3) q(i,j,4) =q(i,j,4)+dt1*d(i,j,4) end do end do c termination return end subroutine wrtd c ************************************************************ c * write field data * c ************************************************************ implicit real*8 (a-h,o-z) parameter(nx=201,ny=201) common /aray1/ q(0:nx,0:ny,4) common /coor2/ imax,jmax common /stat1/ time common /stat2/ loop,locl common /parm3/ gamm,gam1 common /file1/ file_field,file_tecplt character*50 file_field,file_tecplt dimension icom(20),pcom(20) data icom/20*0/pcom/20*0.0/ c set variables of the headder icom(1) =imax icom(2) =jmax icom(3) =loop pcom(1) =time pcom(2) =gamm c write data open(15,file=file_field,form='formatted') write(15,*) (icom(i),i=1,20) write(15,*) (pcom(i),i=1,20) write(15,*) ((q(i,j,1),q(i,j,2),q(i,j,3),q(i,j,4),i=1,imax), & j=1,jmax) close(15) c termination return end subroutine tecp c ************************************************************** c * sample program for tecplot data * c ************************************************************** implicit real*8 (a-h,o-z) parameter(nx=201,ny=201) common /aray1/ q(0:nx,0:ny,4) common /aray2/ d(0:nx,0:ny,4) common /coor1/ x(nx,ny),y(nx,ny) common /coor2/ imax,jmax common /parm3/ gamm,gam1 common /file1/ file_field,file_tecplt character*50 file_field,file_tecplt dimension w(0:nx,0:ny) c open file open(20,file=file_tecplt,form='formatted',status='unknown') c write headder write(20,501) 501 format(35hVARIABLES="X","Y","RHO","P","S","M") c transfer data (cell center to cell node) do j=0,jmax do i=0,imax w(i,j) =0.0 end do end do do j=1,jmax-1 do i=1,imax-1 w(i,j) =1.0 end do end do do j=1,jmax do i=1,imax ddv =1.0/(w(i-1,j-1)+w(i-1,j)+w(i,j-1)+w(i,j)) cmm =ddv* w(i-1,j-1) cmp =ddv* w(i-1,j ) cpm =ddv* w(i ,j-1) cpp =ddv* w(i ,j ) d(i,j,1) =cmm*q(i-1,j-1,1)+cmp*q(i-1,j ,1) & +cpm*q(i ,j-1,1)+cpp*q(i ,j ,1) d(i,j,2) =cmm*q(i-1,j-1,2)+cmp*q(i-1,j ,2) & +cpm*q(i ,j-1,2)+cpp*q(i ,j ,2) d(i,j,3) =cmm*q(i-1,j-1,3)+cmp*q(i-1,j ,3) & +cpm*q(i ,j-1,3)+cpp*q(i ,j ,3) d(i,j,4) =cmm*q(i-1,j-1,4)+cmp*q(i-1,j ,4) & +cpm*q(i ,j-1,4)+cpp*q(i ,j ,4) end do end do c write data write(20,502) imax,jmax 502 format(7hZONE I=,i4,3h,J=,i4,1x,7hF=POINT) do j=1,jmax do i=1,imax xg =x(i,j) yg =y(i,j) rov =1.0/d(i,j,1) pc =gam1*(d(i,j,4)-0.5*rov*(d(i,j,2)**2+d(i,j,3)**2)) cc2 =gamm*pc*rov xm =dsqrt((d(i,j,2)**2+d(i,j,3)**2)*rov**2/cc2) sc =pc*rov**gamm write(20,503) xg,yg,d(i,j,1),pc,sc,xm end do end do 503 format(1x,8(e15.7,1h,)) c plot boundary write(20,504) 504 format(33hGEOMETRY T=LINE, CS=GRID, F=POINT) write(20,*) '1' write(20,*) 2*imax+2*jmax-3 do j=1,jmax-1 write(20,505) x(1,j),y(1,j) end do do i=1,imax-1 write(20,505) x(i,jmax),y(i,jmax) end do do j=jmax,1,-1 write(20,505) x(imax,j),y(imax,j) end do do i=imax-1,1,-1 write(20,505) x(i,1),y(i,1) end do 505 format(2e20.12) close(20) c termination return end