tranx1.src

#include "zeus2d.def"
c=======================================================================
c//////////////////////////  SUBROUTINE TRANX1  \\\\\\\\\\\\\\\\\\\\\\\\
c
      subroutine tranx1(mflx1,s3)
c
c  PURPOSE:  Transports all zone centered variables in the 1-direction
c  only.  Currently transported are: density
c                                    energy density
c                                    3-momentum
c                                    3-magnetic field.
c                                    radiation energy density
c  The consistent transport algorithm is used, including the effects of
c  grid compression.  The transported fluxes are thus given by the mass
c  fluxes times the time centered area of the control volume faces times
c  the interpolated variable. Interpolations are performed in X1INTZC.
c
c  INPUT ARGUMENTS:
c   s3 = momentum density in 3-direction
c
c  OUTPUT ARGUMENTS:
c   mflx1 = mass flux in 1-direction
c   s3 = "half"-updated momentum density in 3-direction
c
c  EXTERNALS: X1INTZC
c
c  LOCALS:
c-----------------------------------------------------------------------
      implicit NONE
#include "param.h"
#include "root.h"
#include "grid.h"
#include "field.h"
#include "scratch.h"
      REAL    mflx1(in,jn),s3(in,jn)
c
      integer i,j
      REAL atwid(in),vel1(in),td(in),eod1(in),pr(in),esc
     &    ,dtwid1(in),etwid1(in)
     &    ,dflx1(in) ,eflx1(in)
      equivalence (atwid,wi0) , (vel1,wi1) , (eod1,eflx1,wi2)
     &           ,(dtwid1,wi3)
     &           ,(etwid1,wi4) , (td,dflx1,wi5) , (pr,wi13)
#ifdef ROTATE
      REAL tv3(in),v3twid1(in),s3flx1(in)
      equivalence (v3twid1,wi7) , (tv3,s3flx1,wi8)
#endif
#ifdef MHD
      REAL b3od1(in),b3twid1(in),b3flx1(in)
      equivalence (b3od1,b3flx1,wi9) , (b3twid1,wi10)
#endif
#ifdef RAD
      REAL erod1(in),ertwid1(in),erflx1(in),ersc,sasum
      external sasum
      equivalence (erod1,erflx1,wi11),(ertwid1,wi12)
#endif
      external x1intzc
c\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////////////
c=======================================================================
c  Check for 1-D problem in 2-direction
c  Compute time-centered area factor
c
      if (nx1z .le. 1) return
      do 10 i=is,ie+1
        atwid(i) = g2ah(i)*g31ah(i)
10    continue
      e sc = 1.0
#ifdef RAD
      ersc = sasum(in*jn,er,1)/float(nx1z*nx2z)
      e sc = sasum(in*jn,e ,1)/float(nx1z*nx2z)
#endif
c
c  Interpolate quantities to zone faces in the 1-direction. 
c  vel is the relative fluid velocity at interpolation points.
c
      do 100 j=js,je
        do 20 i=ii(j),iop1(j)
          vel1(i) = v1(i,j) - vg1(i)
20      continue
        do 30 i=iim2(j),iop2(j)
          td   (i) = d (i,j)
          eod1 (i) = e (i,j)/d(i,j)/esc
          if (gamma-1.0.eq.0.0) then
            pr (i) = ciso**2*d(i,j)
          else
            pr (i) = (gamma-1.0)*e(i,j)
          endif
#ifdef ROTATE
          tv3  (i) = v3(i,j)
#endif
#ifdef MHD
          b3od1(i) = b3(i,j)/d(i,j)
#endif
#ifdef RAD
          erod1(i) = er(i,j)/d(i,j)/ersc
#endif
30      continue
c
        call x1intzc(td   ,vel1,pr,j,iordd ,istpd ,d twid1)
        call x1intzc(e od1,vel1,pr,j,iorde ,   0  ,e twid1)
#ifdef ROTATE
        call x1intzc(tv3  ,vel1,pr,j,iords3,   0  ,v3twid1)
#endif
#ifdef MHD
        call x1intzc(b3od1,vel1,pr,j,iordb3,   0  ,b3twid1)
#endif
#ifdef RAD
        call x1intzc(erod1,vel1,pr,j,iorder,   0  ,ertwid1)
#endif
c
c  Construct fluxes at interfaces, including the mass flux which is
c  passed to MOMX1
c
        do 40 i=ii(j),iop1(j)
          m flx1(i,j) = vel1(i)*dt*dtwid1(i)
          d flx1(i) = mflx1(i,j)*atwid(i)
          e flx1(i) = dflx1(i  )*e twid1(i)*esc
#ifdef ROTATE
          s3flx1(i) = dflx1(i  )*v3twid1(i)*g31a(i)*g32b(j)
#endif
#ifdef MHD
          b3flx1(i) = mflx1(i,j)*b3twid1(i)*g2ah(i)
#endif
#ifdef RAD
          erflx1(i) = dflx1(i  )*ertwid1(i)*ersc
#endif
40      continue
c
c  Perform advection using fluxes.  Note timestep dt is hidden in the
c  fluxes via the mass fluxes.
c
        do 50 i=ii(j),io(j)
          d (i,j) = (d (i,j)*dvl1a(i)-(d flx1(i+1)-d flx1(i)))/dvl1an(i)
          e (i,j) = (e (i,j)*dvl1a(i)-(e flx1(i+1)-e flx1(i)))/dvl1an(i)
#ifdef ROTATE
          s3(i,j) = (s3(i,j)*dvl1a(i)-(s3flx1(i+1)-s3flx1(i)))/dvl1an(i)
#endif
#ifdef MHD
          b3(i,j) = (b3(i,j)*g2b(i)*dx1a(i)-(b3flx1(i+1)-b3flx1(i)))
     &                   /(g2bn(i)*dx1an(i))
#endif
#ifdef RAD
          er(i,j) = (er(i,j)*dvl1a(i)-(erflx1(i+1)-erflx1(i)))/dvl1an(i)
#endif
50      continue
100   continue
c
      return
      end