Main Page
src/balansstep.f

c
c     Source code by Sergey Karabasov
c
      subroutine balansstep(mx,my,mbc,meqn,xx,yy,dx,dy,q,ddt,ccfl,fx,fy)
      implicit double precision (a-h,o-z)

      dimension xx(1-mbc:mx+mbc),yy(1-mbc:my+mbc)
      dimension q(meqn,1-mbc:mx+mbc,1-mbc:my+mbc)
      dimension fx(meqn,1-mbc:mx+mbc,1-mbc:my+mbc)
      dimension fy(meqn,1-mbc:mx+mbc,1-mbc:my+mbc)

      COMMON/CELL/ ROCELL(2001,2001),UCELL(2001,2001),VCELL(2001,2001),
     &     ETCELL(2001,2001),EICELL(2001,2001)
      COMMON/GX/UX(2001,2001),VX(2001,2001),ROX(2001,2001),PX(2001,2001)
      COMMON/GY/UY(2001,2001),VY(2001,2001),ROY(2001,2001),PY(2001,2001)
      COMMON/XY/X(2001),Y(2001)

      COMMON/D1X/PN(2001),UN(2001),RON(2001),VN(2001)

      COMMON/FLUXES/FLUX1_X(2001,2001),FLUX1_Y(2001,2001)
     &     ,FLUX2_X(2001,2001),FLUX2_Y(2001,2001)
     &     ,FLUX3_X(2001,2001),FLUX3_Y(2001,2001)
     &     ,FLUX4_X(2001,2001),FLUX4_Y(2001,2001)
 
      COMMON/BAL/NX,NY,CFL,EPS,GAM,DT

c     Things to do:
c     - Don't touch ghost cells, reformulate loops with 1-mbc:mx+mbc
c     - Reformulate solver with q(meqn), fx(meqn), fy(meqn)
c     - Clean up usage of loop bounds
c     - Use parameter lists instead of common blocks, at least for q, f
c     - specify constants as double precision
c     - don't format codes with tabs (this is an extension to the F77 standard)
c
      NX = mx+2*mbc
      NY = my+2*mbc
      DT = ddt
      CFL = 0.d0

      do i=1-mbc,mx+mbc 
         X(i+mbc)=xx(i)
      enddo
         
      do j=1-mbc,my+mbc 
         Y(j+mbc)=yy(j)
      enddo

      do 10 j=1-mbc,my+mbc
         do 10 i=1-mbc,mx+mbc 
            
            II = i+mbc
            JJ = j+mbc
            ROCELL(II,JJ) = q(1,i,j)
            UCELL(II,JJ)  = q(2,i,j)/q(1,i,j)
            VCELL(II,JJ)  = q(3,i,j)/q(1,i,j)
            ETCELL(II,JJ) = q(4,i,j)/q(1,i,j)
            EICELL(II,JJ) = ETCELL(II,JJ)-0.5*UCELL(II,JJ)**2-
     &           0.5*VCELL(II,JJ)**2

            ROX(II,JJ) = q(5,i,j)
            UX(II,JJ)  = q(6,i,j)
            VX(II,JJ)  = q(7,i,j)
            PX(II,JJ)  = q(8,i,j)

            ROY(II,JJ) = q(9,i,j)
            UY(II,JJ)  = q(10,i,j)
            VY(II,JJ)  = q(11,i,j)
            PY(II,JJ)  = q(12,i,j)

 10   continue

      CALL BALANS
      CALL CHARX
      CALL CHARY
      CALL BALANS

      do 20 j=1-mbc,my+mbc
         do 20 i=1-mbc,mx+mbc 
            
            II = i+mbc
            JJ = j+mbc

            q(1,i,j) = ROCELL(II,JJ)
            q(2,i,j) = ROCELL(II,JJ)*UCELL(II,JJ)
            q(3,i,j) = ROCELL(II,JJ)*VCELL(II,JJ)
            q(4,i,j) = ROCELL(II,JJ)*ETCELL(II,JJ)

            q(5,i,j) = ROX(II,JJ)
            q(6,i,j) = UX(II,JJ)
            q(7,i,j) = VX(II,JJ)
            q(8,i,j) = PX(II,JJ)

            q(9,i,j)  = ROY(II,JJ) 
            q(10,i,j) = UY(II,JJ)
            q(11,i,j) = VY(II,JJ) 
            q(12,i,j) = PY(II,JJ) 

            fx(1,i,j) = FLUX1_X(II,JJ)
            fx(2,i,j) = FLUX2_X(II,JJ)
            fx(3,i,j) = FLUX3_X(II,JJ)
            fx(4,i,j) = FLUX4_X(II,JJ)

            fy(1,i,j) = FLUX1_Y(II,JJ)
            fy(2,i,j) = FLUX2_Y(II,JJ)
            fy(3,i,j) = FLUX3_Y(II,JJ)
            fy(4,i,j) = FLUX4_Y(II,JJ)

 20   continue

      ccfl = CFL 

      return
      end