cc---6400, cm=3.0, tstop=1.0, eps=0.95 or 0.9
      parameter(NN=1004)
      DIMENSION DENS(NN), XMOM(NN), ENER(NN), 
     *     FM(NN), FP(NN), FE(NN),
     *     DD(NN),DX(NN),DE(NN)
      COMMON/PAR/ CK, PRAN

      PRAN=1.0
      PI=4.0*atan(1.0)
      GAM = 5.0/3.0
      CK = 2.0
      RGAMH = SQRT(GAM/2.)
      PROT1  = 1.0E-10
      EPS=0.4
      AA=1.0
      EP=1.0E-20

c      write(*,*) 'IIMAX,CM,TSTOP,EPS'
c      read(*,*) IIMAX,CM,TSTOP,EPS

c      write(*,*) 'IIMAX,TSTOP,EPS'
c      read(*,*) IIMAX,TSTOP,EPS

      IIMAX=200
      TSTOP=0.1
      EPS=0.5
     

      UNU=0.0005
C      EPS=0.5
      DXX=1.0/IIMAX     
      
      IMAX=200

      IMAX=IMAX+4

C      TSTOP=0.25


      TIME   = 0.0
      TAUMAX = 10.0
      UBD    = 0.0
      IT     = 1
      OPEN(UNIT=1,FILE='maa.out',STATUS='unknown')
c      OPEN(UNIT=2,FILE='old.out',STATUS='unknown')
C      OPEN(UNIT=3,FILE='new.out',STATUS='unknown')
      
      IMAXM  = IMAX - 1
      IMAXM2 = IMAX - 2

C     INITIAL CONDITIONS:
      

      D1=1.0
      U1=1.0
      AP1=1/(GAM*CM*CM)
      TT1=UNU/AP1

      
      D2=(GAM+1.0)*CM*CM*D1/(2.0+(GAM-1.0)*CM*CM)
      U2=((GAM-1.0)/(GAM+1.0)+2.0/((GAM+1.0)*CM*CM))*U1
      AP2=(2.0*GAM*CM*CM/(GAM+1.0)-(GAM-1.0)/(GAM+1.0))*AP1

      AE00=D1/(2.0*AP1)
      TT2=UNU/AP2

      CU=1.0

c      DO I=3,IMAXM2
c         READ(2,*)  DENS(I), XMOM(I), ENER(I)
c      enddo



      DO I = 2, IMAXM
         DENS(I)=1.0
         XMOM(I)=-2.0
         ENER(I)=3.0

         IF(FLOAT(I-2).GE.100.0) THEN
            DENS(I)=1.0
           XMOM(I)=2.0
            ENER(I)=3.0
         ENDIF


c         DENS(I)=1.0
c         XMOM(I)=0.0
c         ENER(I)=2.5

c         IF(FLOAT(I-2).GE.100.0) THEN
c            DENS(I)=0.125
c           XMOM(I)=0.0
c            ENER(I)=0.25
c         ENDIF

      ENDDO


      XMOM(1) =  XMOM(4)
      XMOM(2) =  XMOM(3)
      XMOM(IMAXM) =  XMOM(IMAXM2)
      XMOM(IMAX) =  XMOM(IMAXM2-1)

      DENS(1) =  DENS(4)
      ENER(1) =  ENER(4)
      DENS(2) =  DENS(3)
      ENER(2) =  ENER(3)
      DENS(IMAXM) =  DENS(IMAXM2)
      ENER(IMAXM) =  ENER(IMAXM2)
      DENS(IMAX) =  DENS(IMAXM2-1)
      ENER(IMAX) =  ENER(IMAXM2-1)

C     COURANT CONDITIONS:

 601  CONTINUE
      DO 201 I = 2, IMAXM
         UU=XMOM(I)/DENS(I)
         ELAM = 0.25 * (CK + 1.0) * DENS(I)/(ENER(I) 
     *        - 0.5 *(XMOM(I)**2)/DENS(I))
         IF(ELAM.LE.0.0) THEN 
            WRITE(*,*) 'ELAM NEGATIVE',D
     *           ENS(I),XMOM(I),ENER(I), I
            stop
         ENDIF
         SIG = 1.0/(SQRT(ELAM))
         UABS = ABS(UU)
         S    = UABS + RGAMH * SIG
         UBD  = MAX(UBD,S)
         TBD  = 1.0/(UBD + PROT1)
 201  CONTINUE

c--- TAU IS THE TIME STEP

      TAU =DXX*(1.0-EPS)*MIN(TAUMAX,TBD)
      UBD = 0.0

      WRITE(*,*) ' IT, TIME, TAU: ',IT,TIME,TAU


C------ OUTPUT (AP: PRESSURE)




      TIME=TIME+TAU
      IT=IT+1

c Initial Reconstrcution for Conservative variables (van Leer limiter is used)
c--- (other limiters can be also used here)

      DO I= 2,IMAXM2      
         I0=I-1
         I1=I+1
         
         DFD1=DENS(I)-DENS(I0)
         DFX1=XMOM(I)-XMOM(I0)
         DFE1=ENER(I)-ENER(I0)
                  
         DFD2=DENS(I1)-DENS(I)
         DFX2=XMOM(I1)-XMOM(I)
         DFE2=ENER(I1)-ENER(I)
         

         DD(I)=(SIGN(AA,DFD1)+SIGN(AA,DFD2))
     *        *((ABS(DFD1)*ABS(DFD2))/(ABS(DFD1)+ABS(DFD2)+EP))/DXX
         DX(I)=(SIGN(AA,DFX1)+SIGN(AA,DFX2))
     *        *((ABS(DFX1)*ABS(DFX2))/(ABS(DFX1)+ABS(DFX2)+EP))/DXX
         DE(I)=(SIGN(AA,DFE1)+SIGN(AA,DFE2))
     *        *((ABS(DFE1)*ABS(DFE2))/(ABS(DFE1)+ABS(DFE2)+EP))/DXX
         
      ENDDO


      IF(TIME.GT.TSTOP) THEN
         
         DO I=4,IMAXM2-1
            XI = FLOAT(I-2)*DXX
            AP=(GAM-1.0)*(ENER(I)
     *           -0.5*XMOM(I)*XMOM(I)/DENS(I))
            AE0=DENS(I)/(2.0*AP)
            CC=SQRT(GAM*AP/DENS(I))
            T=1.0/AE0

            UX=(DX(I)-XMOM(I)*DD(I)/DENS(I))/DENS(I)
           
            UN=UNU*(AE00/AE0)**0.8            
            TNN=4.0*UN*UX/(6.0*AP)
            
            TX=(4.0/(CK+1.0))*(DE(I)+
     *           (XMOM(I)*XMOM(I)/DENS(I)-ENER(I))*DD(I)/DENS(I)
     *           -XMOM(I)*DX(I)/DENS(I))/DENS(I)

            CHNN=1.25*UN*TX/(PRAN*AP*CC)

            WRITE(1,88)  XI, DENS(I), XMOM(I), XMOM(I)/DENS(I)
 88         FORMAT(4E13.5)
         enddo

         DO I=3,IMAXM2           
            WRITE(3,*) DENS(I), XMOM(I), ENER(I)
         enddo

         write(*,*) 'IMAX=', IMAX 
         GOTO 965
         
      ENDIF


c---time-dependent numerical fluxes evaluation at cell interface i+1/2

      TMIN=10000.0
      TMAX=0.0

      DO 999 I= 2,IMAXM2

         I0=I-1
         I1=I+1
         I2=I+2

C--- MASS, MOMENTUM AND ENERGY DENSITIES AT LEFT AND RIGHT SIDES OF A CELL 
C--- INTERFACE

         ADE1=DENS(I)+0.5*DD(I)*DXX
         AXM1=XMOM(I)+0.5*DX(I)*DXX
         AEN1=ENER(I)+0.5*DE(I)*DXX
         
         ADE2=DENS(I1)-0.5*DD(I1)*DXX
         AXM2=XMOM(I1)-0.5*DX(I1)*DXX
         AEN2=ENER(I1)-0.5*DE(I1)*DXX
         
c--- \lambda and U on both sides (this is needed for initial f_0)
      
         AE1=0.25*(CK+1)*ADE1/(AEN1
     *        -0.50*AXM1*AXM1/ADE1)
         AU1=AXM1/ADE1
         
         AE2=0.25*(CK+1)*ADE2/(AEN2
     *        -0.50*AXM2*AXM2/ADE2)
         AU2=AXM2/ADE2


c--- if the initial reconstruction gives negative temperature, ignore the
c--- reconstrcution)

         
         IF(AE1.LE.0.0) THEN
            ADE1=DENS(I)
            AXM1=XMOM(I)
            AEN1=ENER(I)
            DD(I)=0.0
            DX(I)=0.0
            DE(I)=0.0
            AE1=0.25*(CK+1)*ADE1/(AEN1
     *           -0.50*AXM1*AXM1/ADE1)
            AU1=AXM1/ADE1
         ENDIF

         IF(AE2.LE.0.0) THEN
            ADE2=DENS(I1)
            AXM2=XMOM(I1)
            AEN2=ENER(I1)
            DD(I1)=0.0
            DX(I1)=0.0
            DE(I1)=0.0
            AE2=0.25*(CK+1)*ADE2/(AEN2
     *           -0.50*AXM2*AXM2/ADE2)
            AU2=AXM2/ADE2
         ENDIF
         
C--- From slopes in the macroscopic variables to obtain the slopes in
c--- the gas-distribution function of f_0.         

         CALL TRAN(AU1,AE1,DD(I),DX(I),DE(I),X1,Y1,Z1)
         CALL TRAN(AU2,AE2,DD(I1),DX(I1),DE(I1),X2,Y2,Z2)


c-- define the moments for g_L and g_R IN f_0

         RAE1=1.0/AE1
         RAE2=1.0/AE2
         
         TEU0=0.5*DERFC(-AU1*SQRT(AE1))
         TEU1=AU1*TEU0+0.50*EXP(-AE1*AU1*AU1)/SQRT(AE1*PI)
         TEU2=AU1*TEU1+0.5*TEU0*RAE1
         TEU3=AU1*TEU2+1.0*TEU1*RAE1
         TEU4=AU1*TEU3+1.5*TEU2*RAE1
         TEU5=AU1*TEU4+2.0*TEU3*RAE1
         TEU6=AU1*TEU5+2.5*TEU4*RAE1
          
         TGU0=0.5*(DERFC(AU2*SQRT(AE2)))
         TGU1=AU2*TGU0-0.50*EXP(-AE2*AU2*AU2)/SQRT(AE2*PI)
         TGU2=AU2*TGU1+0.5*TGU0*RAE2
         TGU3=AU2*TGU2+1.0*TGU1*RAE2
         TGU4=AU2*TGU3+1.5*TGU2*RAE2
         TGU5=AU2*TGU4+2.0*TGU3*RAE2
         TGU6=AU2*TGU5+2.5*TGU4*RAE2

C-----------
         TLU0=1.0
         TLU1=AU1
         TLU2=AU1*TLU1+0.5*TLU0*RAE1
         TLU3=AU1*TLU2+1.0*TLU1*RAE1
         TLU4=AU1*TLU3+1.5*TLU2*RAE1
         TLU5=AU1*TLU4+2.0*TLU3*RAE1
         TLU6=AU1*TLU5+2.5*TLU4*RAE1

         TRU0=1.0
         TRU1=AU2
         TRU2=AU2*TRU1+0.5*TRU0*RAE2
         TRU3=AU2*TRU2+1.0*TRU1*RAE2
         TRU4=AU2*TRU3+1.5*TRU2*RAE2
         TRU5=AU2*TRU4+2.0*TRU3*RAE2
         TRU6=AU2*TRU5+2.5*TRU4*RAE2
          
         TL02=0.5*CK*RAE1
         TL04=0.25*CK*(CK+2.0)*RAE1*RAE1

         TR02=0.5*CK*RAE2
         TR04=0.25*CK*(CK+2.0)*RAE2*RAE2


         ADE=ADE1*TEU0+ADE2*TGU0
         AXM=ADE1*TEU1+ADE2*TGU1
         AEN=0.5*(ADE1*(TEU2+TEU0*TL02)
     *        +ADE2*(TGU2+TGU0*TR02))

cc-- \lambda_0 and U_0 in g_0

         AE0=0.25*(CK+1)*ADE/(AEN-0.50*AXM*AXM/ADE)
         AU0=AXM/ADE

c--   collision time (normalized by time-step)         
         APP=0.5*ADE/AE0
         APP1=0.5*ADE1/AE1
         APP2=0.5*ADE2/AE2
         TE=0.05*TAU+
     *        +TAU*(ABS(APP1-APP2)/(APP1+APP2))         

c         TE=UNU*(AE00/AE0)**0.8/APP
c     *        +TAU*(ABS(APP1-APP2)/(APP1+APP2))         
c         TE=UNU/APP+TAU*(ABS(APP1-APP2)/(APP1+APP2))         


c--   some time integrals (gammas)

         TMIN=MIN(TMIN,TE/TAU)
         TMAX=MAX(TMAX,TE/TAU)
         

         SP=-TAU/TE
         SW=EXP(SP)
         SE1=TE*(1.0-SW)
         SET1=-TE*TAU*SW+TE*SE1

c-- the moments for f0(-ut)         

         ETRI1U0= X1*TEU0+Y1*TEU1+Z1*(TEU2+TEU0*TL02)
         ETRI1U1= X1*TEU1+Y1*TEU2+Z1*(TEU3+TEU1*TL02)
         ETRI1U2= X1*TEU2+Y1*TEU3+Z1*(TEU4+TEU2*TL02)
         ETRI1U3= X1*TEU3+Y1*TEU4+Z1*(TEU5+TEU3*TL02)
         ETRI1EU0=0.5*(X1*(TEU2+TEU0*TL02)+Y1*(TEU3+TEU1*TL02)+
     *        Z1*(TEU4+2.0*TEU2*TL02+TEU0*TL04))         
         ETRI1EU1=0.5*(X1*(TEU3+TEU1*TL02)+Y1*(TEU4+TEU2*TL02)+
     *        Z1*(TEU5+2.0*TEU3*TL02+TEU1*TL04))
         ETRI1EU2=0.5*(X1*(TEU4+TEU2*TL02)+Y1*(TEU5+TEU3*TL02)+
     *        Z1*(TEU6+2.0*TEU4*TL02+TEU2*TL04))

         GTRI1U0= X2*TGU0+Y2*TGU1+Z2*(TGU2+TGU0*TR02)
         GTRI1U1= X2*TGU1+Y2*TGU2+Z2*(TGU3+TGU1*TR02)
         GTRI1U2= X2*TGU2+Y2*TGU3+Z2*(TGU4+TGU2*TR02)
         GTRI1U3= X2*TGU3+Y2*TGU4+Z2*(TGU5+TGU3*TR02)
         GTRI1EU0=0.5*(X2*(TGU2+TGU0*TR02)+Y2*(TGU3+TGU1*TR02)+
     *        Z2*(TGU4+2.0*TGU2*TR02+TGU0*TR04))         
         GTRI1EU1=0.5*(X2*(TGU3+TGU1*TR02)+Y2*(TGU4+TGU2*TR02)+
     *        Z2*(TGU5+2.0*TGU3*TR02+TGU1*TR04))
         GTRI1EU2=0.5*(X2*(TGU4+TGU2*TR02)+Y2*(TGU5+TGU3*TR02)+
     *        Z2*(TGU6+2.0*TGU4*TR02+TGU2*TR04))


         FMI1=ADE1*SE1*TEU0-SET1*ETRI1U1
         FPI1=ADE1*SE1*TEU1-SET1*ETRI1U2
         FEI1=ADE1*0.5*SE1*(TEU2+TEU0*TL02)-SET1*ETRI1EU1

         FMI2=ADE2*SE1*TGU0-SET1*GTRI1U1
         FPI2=ADE2*SE1*TGU1-SET1*GTRI1U2
         FEI2=ADE2*0.5*SE1*(TGU2+TGU0*TR02)-SET1*GTRI1EU1




c--- fluxes from f_0 (the non-equilibrium state )

         FM1=ADE1*SE1*TEU1-(SET1+SE1*TE)*ETRI1U2
         FP1=ADE1*SE1*TEU2-(SET1+SE1*TE)*ETRI1U3
         FE1=ADE1*0.5*SE1*(TEU3+TEU1*TL02)-(SET1+SE1*TE)*ETRI1EU2
         FM2=ADE2*SE1*TGU1-(SET1+SE1*TE)*GTRI1U2
         FP2=ADE2*SE1*TGU2-(SET1+SE1*TE)*GTRI1U3
         FE2=ADE2*0.5*SE1*(TGU3+TGU1*TR02)-(SET1+SE1*TE)*GTRI1EU2


         A1=-(X1*TLU1+Y1*TLU2+Z1*
     *        (TLU3+TLU1*TL02))
         A2=-(X1*TLU2+Y1*TLU3+Z1*
     *        (TLU4+TLU2*TL02))
         A3=-0.5*(X1*(TLU3+TLU1*TL02)+
     *        Y1*(TLU4+TLU2*TL02)+
     *        Z1*(TLU5+2.0*TLU3*TL02+TLU1*TL04))

         B1=-(X2*TRU1+Y2*TRU2+Z2*
     *        (TRU3+TRU1*TR02))
         B2=-(X2*TRU2+Y2*TRU3+Z2*
     *        (TRU4+TRU2*TR02))
         B3=-0.5*(X2*(TRU3+TRU1*TR02)+
     *        Y2*(TRU4+TRU2*TR02)+
     *        Z2*(TRU5+2.0*TRU3*TR02+TRU1*TR04))

         CALL TRAN(AU1,AE1,A1,A2,A3,A11,B11,C11)
         CALL TRAN(AU2,AE2,B1,B2,B3,A22,B22,C22)

         AETRI1U1=A11*TEU1+B11*TEU2+C11*(TEU3+TEU1*TL02)
         AETRI1U2=A11*TEU2+B11*TEU3+C11*(TEU4+TEU2*TL02)         
         AETRI1U3=A11*TEU3+B11*TEU4+C11*(TEU5+TEU3*TL02)         
         AETRI1EU1=0.5*(A11*(TEU3+TEU1*TL02)+B11*(TEU4+
     *        TEU2*TL02)+C11*(TEU5+2.0*TEU3*TL02+TEU1*TL04))

         AGTRI1U1=A22*TGU1+B22*TGU2+C22*(TGU3+TGU1*TR02)
         AGTRI1U2=A22*TGU2+B22*TGU3+C22*(TGU4+TGU2*TR02)         
         AGTRI1U3=A22*TGU3+B22*TGU4+C22*(TGU5+TGU3*TR02)         
         AGTRI1EU1=0.5*(A22*(TGU3+TGU1*TR02)+B22*(TGU4+
     *        TGU2*TR02)+C22*(TGU5+2.0*TGU3*TR02+TGU1*TR04))



         AMI=FMI1+FMI2
         API=FPI1+FPI2
         AEI=FEI1+FEI2
          

         FMLL=(-SE1*TE)*AETRI1U1
         FPLL=(-SE1*TE)*AETRI1U2
         FELL=(-SE1*TE)*AETRI1EU1

         FMRR=(-SE1*TE)*AGTRI1U1
         FPRR=(-SE1*TE)*AGTRI1U2
         FERR=(-SE1*TE)*AGTRI1EU1


c--- the slopes for the macroscopic variables across cell interface

         AW1=2.0*(ADE-DENS(I))/DXX
         BW1=2.0*(AXM-XMOM(I))/DXX
         CW1=2.0*(AEN-ENER(I))/DXX
         
         AW2=-2.0*(ADE-DENS(I1))/DXX
         BW2=-2.0*(AXM-XMOM(I1))/DXX
         CW2=-2.0*(AEN-ENER(I1))/DXX
         
c--- the slopes in g_0 at both sides of the cell interface

         CALL TRAN(AU0,AE0,AW1,BW1,CW1,XL,YL,ZL)
         CALL TRAN(AU0,AE0,AW2,BW2,CW2,XR,YR,ZR)
         
         RAE0=1.0/AE0
         

c-- define some moments for g0

         T0=1.0
         T1=AU0
         T2=AU0*T1+0.5*T0*RAE0
         T3=AU0*T2+1.0*T1*RAE0
         T4=AU0*T3+1.5*T2*RAE0
         T5=AU0*T4+2.0*T3*RAE0
         T6=AU0*T5+2.5*T4*RAE0

         TLU0=0.5*DERFC(-AU0*SQRT(AE0))
         TLU1=AU0*TLU0+0.50*EXP(-AE0*AU0*AU0)/SQRT(AE0*PI)
         TLU2=AU0*TLU1+0.5*TLU0*RAE0
         TLU3=AU0*TLU2+1.0*TLU1*RAE0
         TLU4=AU0*TLU3+1.5*TLU2*RAE0
         TLU5=AU0*TLU4+2.0*TLU3*RAE0
         TLU6=AU0*TLU5+2.5*TLU4*RAE0
          
         TRU0=1.0-TLU0
         TRU1=T1-TLU1
         TRU2=T2-TLU2
         TRU3=T3-TLU3
         TRU4=T4-TLU4
         TRU5=T5-TLU5
         TRU6=T6-TLU6

         T02=0.5*CK*RAE0
         T04=0.25*CK*(CK+2.0)*RAE0*RAE0

C--- moments from uag_0
          
         TRI1U1=(XL*TLU1+YL*TLU2+ZL*(TLU3+TLU1*T02))+
     *        (XR*TRU1+YR*TRU2+ZR*(TRU3+TRU1*T02))
         TRI1U2=(XL*TLU2+YL*TLU3+ZL*(TLU4+TLU2*T02))+
     *        (XR*TRU2+YR*TRU3+ZR*(TRU4+TRU2*T02))
         TRI1U3=(XL*TLU3+YL*TLU4+ZL*(TLU5+TLU3*T02))+
     *        (XR*TRU3+YR*TRU4+ZR*(TRU5+TRU3*T02))
         TRI1EU1=0.5*(XL*(TLU3+TLU1*T02)+YL*(TLU4+TLU2*T02)+
     *        ZL*(TLU5+2.0*TLU3*T02+TLU1*T04))+
     *        0.5*(XR*(TRU3+TRU1*T02)+YR*(TRU4+TRU2*T02)+
     *        ZR*(TRU5+2.0*TRU3*T02+TRU1*T04))
         TRI1EU2=0.5*(XL*(TLU4+TLU2*T02)+YL*(TLU5+TLU3*T02)+
     *        ZL*(TLU6+2.0*TLU4*T02+TLU2*T04)) +
     *        0.5*(XR*(TRU4+TRU2*T02)+YR*(TRU5+TRU3*T02)+
     *        ZR*(TRU6+2.0*TRU4*T02+TRU2*T04))


C-- time variations of macroscopic variables at cell interface          

         AS=(-TRI1U1*(TE*(-TAU+SE1)+SET1)
     *        +(SE1*ADE-AMI))/(TE*(-TAU+SE1))
         BS=(-TRI1U2*(TE*(-TAU+SE1)+SET1)
     *        +(SE1*AXM-API))/(TE*(-TAU+SE1))
         CS=(-TRI1EU1*(TE*(-TAU+SE1)+SET1)
     *        +(SE1*AEN-AEI))/(TE*(-TAU+SE1))
         
c--- time variations in the distribution function (A_1, A_2, A_3)

         CALL TRAN(AU0,AE0,AS,BS,CS,A1,B1,C1)

c--- moments for {A g_0}
         
         ATRI1U1=A1*AU0+B1*T2+C1*(T3+T1*T02)
         ATRI1U2=A1*T2+B1*T3+C1*(T4+T2*T02)
         ATRI1EU1=0.5*(A1*(T3+T1*T02)+B1*(T4+
     *        T2*T02)+C1*(T5+2.0*T3*T02+T1*T04))
         

          A=ADE*(TAU-SE1)
          B=TE*(-TAU+SE1)+SET1
          D=0.5*TAU*TAU-TE*TAU+TE*SE1

          FM(I)=A*T1+B*TRI1U2+
     *         D*ATRI1U1
     *         +(FM1+FM2+FMLL+FMRR)
          FP(I)=A*T2+B*TRI1U3+
     *         D*ATRI1U2
     *         +(FP1+FP2+FPLL+FPRR)
          FE(I)=0.5*A*(T3+T1*T02)+
     *         B*TRI1EU2+
     *         D*ATRI1EU1
     *         +(FE1+FE2+FELL+FERR)
     *         +TAU*(1.0/PRAN-1.0)*(-TE)*
     *         (TRI1EU2+ATRI1EU1-T1*(TRI1U3+ATRI1U2))



 999   CONTINUE




       DO I = 1, IMAX
          DENS(I) = DENS(I)+(FM(I-1)-FM(I))/DXX
          XMOM(I) = XMOM(I)+(FP(I-1)-FP(I))/DXX
          ENER(I) = ENER(I)+(FE(I-1)-FE(I))/DXX
       ENDDO


c--- boundary condition

       XMOM(1) =  XMOM(4)
       XMOM(2) =  XMOM(3)
       XMOM(IMAXM) =  XMOM(IMAXM2)
       XMOM(IMAX) =  XMOM(IMAXM2-1)


       DENS(1) =  DENS(4)
       ENER(1) =  ENER(4)
       DENS(2) =  DENS(3)
       ENER(2) =  ENER(3)
       DENS(IMAXM) =  DENS(IMAXM2)
       ENER(IMAXM) =  ENER(IMAXM2)
       DENS(IMAX) =  DENS(IMAXM2-1)
       ENER(IMAX) =  ENER(IMAXM2-1)


       WRITE(*,*)'TE/TAU:', TMIN,'MAX', TMAX
       GOTO 601
       
 965   STOP
       END
       

      
      SUBROUTINE TRAN(AU,AL,A,B,C,XQ,YQ,ZQ)
      COMMON/PAR/ CK, PRAN
      ZQ=2.0*AL*AL*(2.0*C-2.0*AU*B+(AU*AU-0.5*(CK+1.0)/AL)*A)
     *     /(CK+1.0)
      YQ=2.0*AL*(B-AU*A)-2.0*AU*ZQ
      XQ=A-YQ*AU-ZQ*(AU*AU+0.5*(CK+1.0)/AL)
      RETURN
      END

