      PROGRAM VANLEER1
C
C     PROGRAM TO TEST THE VAN LEER SCHEME AS DEVELOPED THROUGH 
C     VAN LEER # 2
C     THIS SCHEME IS SIMILAR IN PHILOSOPHY TO PDM AND HARTEN.
C     SEE ALSO TAKACS AND FROMM
C
      PARAMETER MAG=5
      PARAMETER J$=40*MAG
      REAL*4 RHON(J$), RHON1(J$), S(J$), DIFF(J$)
     1, V(J$), X(J$)
      LOGICAL LAX
C
C     SET UP THE INITIAL CONDITION
C
C
C     J$=40 IS THE STANDARD TEST.  TO DOUBLE THE RESOLUTION I HAVE TO 
C     RESOLVE THE SAME DISTRIBUTION WITH TWICE AS MANY POINTS.
C     THIS MEANS SET UP THE INITIAL CONDITON APPROPRIATELY.  IT ALSO MEANS
C     DELTA T MUST BE CHANGED AND HENCE TOTAL INTEGRATION.
C
      aj=float(j$)
C
      T=0.
      TOUT=0.
      DO 100 J=1,J$
      RHON(J)=0.
      X(J)=FLOAT(J)
      V(J)=1.
  100 CONTINUE
      DO 125 J=MAG*15,MAG*25
      RHON(J)=.5
      RHON(J)=.5-.5/(MAG*5.)*ABS(MAG*20.-J)
  125 CONTINUE
      TOTMASS=0.
      DO 150 J=1,J$
      TOTMASS=TOTMASS+RHON(J)/FLOAT(MAG)
  150 CONTINUE
      WRITE(6,11) T,TOTMASS
      write (25) t,rhon
C
C     SET EPS=U DELTT/ DELTX <.5
C
      EPS=.5
      DELTX=1./FLOAT(MAG)
      DELTT=DELTX*EPS
      J$1=J$-1
 1000 CONTINUE
C
C     FORM THE DIFFERENCES
C
      DO J=1,J$1
      DIFF(J)=RHON(J+1)-RHON(J)
      END DO
      J=J$
      DIFF(J)=RHON(1)-RHON(J)
C
C     CALCULATE THE SMOOTHNESS PARAMETER
C
      J=1
      CHK=ABS(DIFF(J))+ABS(DIFF(J$))
      IF (CHK.EQ.0.) THEN
      S(J)=0.
      ELSE
      S(J)=(ABS(DIFF(J))-ABS(DIFF(J$)))/CHK
      END IF
      DO J=2,J$
      CHK=ABS(DIFF(J))+ABS(DIFF(J-1))
      IF (CHK.EQ.0.) THEN
      S(J)=0.
      ELSE
      S(J)=(ABS(DIFF(J))-ABS(DIFF(J-1)))/CHK
      END IF
      END DO
C
C     NOW DO THE ADVECTION
C
      J=1
      RHON1(J)=RHON(J)-V(J)*EPS*DIFF(J$)
     1                -V(J)*EPS/4.*(1.-V(J)*EPS)*(DIFF(J)-DIFF(J$1))
     2                +V(J)*EPS/4.*(1.-V(J)*EPS)*(
     3                 S(J)*(DIFF(J)-DIFF(J$))-
     4                 S(J$)*(DIFF(J$)-DIFF(J$1)))
      J=2
      RHON1(J)=RHON(J)-V(J)*EPS*DIFF(J-1)
     1                -V(J)*EPS/4.*(1.-V(J)*EPS)*(DIFF(J)-DIFF(J$))
     2                +V(J)*EPS/4.*(1.-V(J)*EPS)*(
     3                 S(J)*(DIFF(J)-DIFF(J-1))-
     4                 S(J-1)*(DIFF(J-1)-DIFF(J$)))
      DO J=3,J$
      RHON1(J)=RHON(J)-V(J)*EPS*DIFF(J-1)
     1                -V(J)*EPS/4.*(1.-V(J)*EPS)*(DIFF(J)-DIFF(J-2))
     2                +V(J)*EPS/4.*(1.-V(J)*EPS)*(
     3                 S(J)*(DIFF(J)-DIFF(J-1))-
     4                 S(J-1)*(DIFF(J-1)-DIFF(J-2)))
      END DO
C
C     SET UP NEXT TIME STEP
C
      DO 425 J=1,J$
      RHON(J)=RHON1(J)
  425 CONTINUE
      T=T+DELTT
      TOUT=TOUT+DELTT
      TOTMASS=0.
      DO J=1,J$
      TOTMASS=TOTMASS+RHON1(J)/FLOAT(MAG)
      END DO
      WRITE(6,11) T,TOTMASS
      IF (IFIX(TOUT).EQ.400) THEN
      TOUT=0.
      write (25) t,rhon
c      CALL RSCALE(60,20)
c      CALL NOGRID
c      CALL USXUSY(1.,aj,.FALSE.,-.5,1.,.FALSE.
c     1 ,'VAN LEER (THROUGH II)$')
c      CALL IPLOT(X,RHON,J$,'*',.FALSE.)
      IF(ifix(T).eq.400) STOP50
      GO TO 1000
      ELSE
      IF(T.GE.500) STOP50
      GO TO 1000
      END IF
   10 FORMAT(' TIME=',F10.1,'   GRID POINT=',I4,'  DENSITY=',F10.3)
   11 FORMAT(' TIME=',F10.1,'  TOTAL MASS=',F10.3)
   20 FORMAT(L5)
      END
