fluxsw.F File Reference

#include "implicit_f.inc"
#include "com08_c.inc"
#include "units_c.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	fluxsw (nsflsw, neflsw, nnflsw, crflsw, flsw, x, v, elbuf)

Function/Subroutine Documentation

fluxsw()

subroutine fluxsw	(	integer	nsflsw,
		integer, dimension(*)	neflsw,
		integer, dimension(8,*)	nnflsw,
			crflsw,
			flsw,
			x,
			v,
			elbuf )

Definition at line 28 of file fluxsw.F.

C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER NSFLSW
      INTEGER NEFLSW(*), NNFLSW(8,*)
C     REAL
      my_real
     .   crflsw(6,*), flsw(9,*), x(3,*), v(3,*), elbuf(*)
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com08_c.inc"
#include      "units_c.inc"
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I2, IS, NEL, I1, I, IN, N1, N2, N3, N4, ND, IB2, IB3, IB4,
     .   IB5, IB6, IB10, IB11, IB12, J
C     REAL
      my_real
     .   crx, cry, crz, cnx, cny, cnz, srt, sv2, flvm, pre, ene, rho,
     .   vol, tke, tem, tde, x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4,
     .   z4, vx1, vy1, vz1, vx2, vy2, vz2, vx3, vy3, vz3, vx4, vy4, vz4,
     .   vmx, vmy, vmz, omx, omy, omz, sfx, sfy, sfz, vl2, vlm, sfm,
     .   flv, sw1, sw2
C-----------------------------------------------
      i2=0
      DO 500 is = 1, nsflsw
C
        flsw(1,is) = zero
        flsw(2,is) = zero
        flsw(3,is) = zero
        flsw(4,is) = zero
        flsw(5,is) = zero
        flsw(6,is) = zero
        flsw(7,is) = zero
        flsw(8,is) = zero
        flsw(9,is) = zero        
C
        nel = neflsw(is)
        i1=i2+1
        i2=i2+nel
        crx = crflsw(1,is)
        cry = crflsw(2,is)
        crz = crflsw(3,is)
        cnx = crflsw(4,is)
        cny = crflsw(5,is)
        cnz = crflsw(6,is)
        srt = zero
        sv2 = zero
        flvm = zero        
C
        DO 400 i = i1, i2
          in = nnflsw(1,i)
          n1 = nnflsw(2,i)
          n2 = nnflsw(3,i)
          n3 = nnflsw(4,i)
          n4 = nnflsw(5,i)
          nd = nnflsw(6,i)
          ib2 = nnflsw(7,i)
          ib3 = nnflsw(8,i)
C
          ib4 = ib3 +  in
          ib5 = ib4 +  in
          ib6 = ib5 +  in
          ib10= ib6 +  in
          ib11= ib10+  in
          ib12= ib11+  in
          pre = -(elbuf(ib2)+elbuf(ib2+1)+elbuf(ib2+2))*third
          ene = elbuf(ib3)
          rho = elbuf(ib4)
          vol = elbuf(ib6)
          tke = elbuf(ib10)
          tem = elbuf(ib11)
          tde = elbuf(ib12)
C
          x1 = x(1,n1)
          y1 = x(2,n1)
          z1 = x(3,n1)
C
          x2 = x(1,n2)
          y2 = x(2,n2)
          z2 = x(3,n2)
C
          x3 = x(1,n3)
          y3 = x(2,n3)
          z3 = x(3,n3)
C
          x4 = x(1,n4)
          y4 = x(2,n4)
          z4 = x(3,n4)
C
          vx1 = v(1,n1)
          vy1 = v(2,n1)
          vz1 = v(3,n1)
C
          vx2 = v(1,n2)
          vy2 = v(2,n2)
          vz2 = v(3,n2)
C
          vx3 = v(1,n3)
          vy3 = v(2,n3)
          vz3 = v(3,n3)
C
          vx4 = v(1,n4)
          vy4 = v(2,n4)
          vz4 = v(3,n4)
C
          vmx = (vx1+vx2+vx3+vx4)/nd
          vmy = (vy1+vy2+vy3+vy4)/nd
          vmz = (vz1+vz2+vz3+vz4)/nd
C
          omx = .25*(x1 +x2 +x3 +x4) - crx
          omy = .25*(y1 +y2 +y3 +y4) - cry
          omz = .25*(z1 +z2 +z3 +z4) - crz
C
          sfx = half*((y3-y1)*(z4-z2)-
     1                 (z3-z1)*(y4-y2))
          sfy =half*((z3-z1)*(x4-x2)-
     1                 (x3-x1)*(z4-z2))
                    sfz = half*((x3-x1)*(y4-y2)-
     1                 (y3-y1)*(x4-x2))
C
          vl2 = vmx*vmx+vmy*vmy+vmz*vmz
          vlm = sqrt(vl2)
          sfm = sqrt(sfx*sfx+sfy*sfy+sfz*sfz)
          flv = vmx*sfx + vmy*sfy + vmz*sfz
          sw1 = sfx * (omy*vmz - omz*vmy) +
     1             sfy * (omz*vmx - omx*vmz) +
     2             sfz * (omx*vmy - omy*vmx)
          sw2 = (omx*omx + omy*omy + omz*omz)*
     1             (cnx   *sfx + cny   *sfy + cnz   *sfz)-
     2             (omx*cnx    + omy*cny    + omz*cnz   )*
     3             (omx*sfx + omy*sfy + omz*sfz)
C
          srt        = srt        + sfm
          sv2        = sv2        + sfm*vl2
          flvm       = flvm       + flv
          flsw(1,is) = flsw(1,is) + flv*rho
          flsw(2,is) = flsw(2,is) + sw1
          flsw(3,is) = flsw(3,is) + sw2
          flsw(4,is) = flsw(4,is) + vlm*sfm
          flsw(5,is) = flsw(5,is) + rho*sfm
          flsw(6,is) = flsw(6,is) + pre*sfm
          flsw(7,is) = flsw(7,is) + ene*sfm
          flsw(8,is) = flsw(8,is) + tke*sfm
          flsw(9,is) = flsw(9,is) + tde*sfm
C
  400   CONTINUE
C
        flsw(2,is) = flsw(2,is)/flsw(3,is)
        flsw(4,is) = flsw(4,is)/srt
        flsw(5,is) = flsw(5,is)/srt
        flsw(6,is) = flsw(6,is)/srt
        flsw(7,is) = flsw(7,is)/srt
        flsw(8,is) = flsw(8,is)/srt
        flsw(9,is) = flsw(9,is)/srt
C
        flsw(3,is) = sqrt( sv2/srt - (flvm/srt)**2 )
        IF(flsw(4,is)==zero)THEN
         flsw(3,is)=zero         
        ELSE
         flsw(3,is)=flsw(3,is)/flsw(4,is)
        ENDIF
C
  500 CONTINUE
C
      WRITE (iout, 1000) tt
      DO 600 is = 1, nsflsw
        WRITE (iout,1100) is, (flsw(j,is),j=1,9)
  600 CONTINUE
C
 1000 FORMAT (3h t=,1pe9.3,11h  m. flux  ,11h  rotation ,
     .                     11h  std error,11h  velocity ,
     .                     11h  density  ,11h  pressure ,
     .                     11h  energie  ,11h  ener. t. ,
     .                     11h  diss. t. )
 1100 FORMAT (4h set,i5,3x,1p9e11.3)
      RETURN