      PROGRAM SLOPES
C     
C     THIS IS MY PROGRAM TO STUDY THE SLOPES SCHEME.
C     RUSSELL AND LERNER, JAM, 1981
C     
C
C     MY NOTATION IS AN ANALOGUE TO THAT IN RUSSELL AND LERNER.
C
      PARAMETER J$=40
      REAL*4 RN(J$), RN1(J$), V(J$), X(J$), P(J$), E(J$)
     1,      SN(J$), SN1(J$), F(J$), A(J$), M(J$), RHO(J$) 
     2,      MRN1(J$), MSN1(J$) 
      LOGICAL LDELT
C
C     SET UP THE INITIAL CONDITION
C
C
      J$1=J$-1
      T=0.
      ISTEP=0
      DO 100 J=1,J$
      RN(J)=0.
      X(J)=FLOAT(J)
      RHO(J)=1.
      V(J)=1.
  100 CONTINUE
      DO 125 J=5,15
      RN(J)=.5
      RN(J)=.5-.4/5.*ABS(10.-J)
  125 CONTINUE
C
C     SET EPS=U DELTT/ DELTX <.5
C
      PRINT 25
      READ(5,20) LDELT
      IF (LDELT) THEN
      READ (5,26) DELTT
      ELSE
      DELTT=.5
      END IF
      DELTX=1.
      IOSTEP=5/DELTT
      ISTOP=50/DELTT
      EPS= DELTT/DELTX
      WRITE(6,30) EPS
C
C     CALCULATE THE INITIAL SLOPES SUCH THAT 
C     RN(J)+.5*SN(J)=.5*(R(J)+R(J+1)
C                   = AVERAGE CONC. ON THE RHS OF THE BOX.
C
      DO J=1,J$1
      SN(J)=RN(J+1)-RN(J)
      END DO
      J=J$
      SN(J)=RN(1)-RN(J)
C
C      EXPERIMENT WITH THE SLOPE
C
C       DO J=2,J$
C       SN(J)=RN(J-1)-RN(J)
C       END DO
C       J=1
C       SN(J)=RN(J$)-RN(J)
 1000 CONTINUE
C
C     CALCULATE THE A AND M FUNCTIONS WHICH ARE SORT OF MASS FLUXES
C
      DO J=1,J$
      A(J)=RHO(J)*V(J)*DELTT
      M(J)=RHO(J)*DELTX
      END DO
C
C     CALCULATE THE F FUNCTION
C
      DO J=1,J$1
      IF (A(J).GE.0.) THEN
      F(J)=A(J)*(RN(J)+.5*(1.-A(J)/M(J))*SN(J))
      ELSE
      F(J)=A(J)*(RN(J+1)-.5*(1.+A(J)/M(J+1))*SN(J+1))
      END IF
      END DO
      J=J$
      IF (A(J).GE.0.) THEN
      F(J)=A(J)*(RN(J)+.5*(1.-A(J)/M(J))*SN(J))
      ELSE
      F(J)=A(J)*(RN(1)-.5*(1.+A(J)/M(1))*SN(1))
      END IF
C
C     CALCULATE THE P FUNCTION
C
      DO J=1,J$1
      IF (A(J).GE.0.) THEN
      P(J)=A(J)*(A(J)*A(J)*SN(J)/M(J)-6.*F(J))
      ELSE 
      P(J)=A(J)*(A(J)*A(J)*SN(J+1)/M(J+1)-6*F(J))
      END IF
      END DO
      J=J$    
      IF (A(J).GE.0.) THEN
      P(J)=A(J)*(A(J)*A(J)*SN(J)/M(J)-6.*F(J))
      ELSE 
      P(J)=A(J)*(A(J)*A(J)*SN(1)/M(1)-6*F(J))
      END IF
C
C     CALCULATE THE E FUNCTION
C
      DO J=1,J$1
      IF (A(J).GE.0.) THEN
      E(J)=A(J)*SN(J)
      ELSE
      E(J)=A(J)*SN(J+1)
      END IF
      END DO
      J=J$
      IF (A(J).GE.0.) THEN
      E(J)=A(J)*SN(J)
      ELSE
      E(J)=A(J)*SN(J+1)
      END IF
C
C     NOW DO THE ADVECTION
C
      DO J=2,J$
      MRN1(J)=M(J)*RN(J)-F(J)+F(J-1)
      MSN1(J)=M(J)*SN(J)-(
     1             M(J)*(6.*(F(J-1)+F(J))-(A(J)-A(J-1))*SN(J))
     2            -6.*(A(J-1)+A(J))*(M(J)*RN(J)-F(J)+F(J-1))
     3            +P(J)-P(J-1))/(M(J)-A(J)+A(J-1))
     4            -E(J)+E(J-1)
      END DO
      J=1
      MRN1(J)=M(J)*RN(J)-F(J)+F(J$)
      MSN1(J)=M(J)*SN(J)-(
     1             M(J)*(6.*(F(J$)+F(J))-(A(J)-A(J$))*SN(J))
     2            -6.*(A(J$)+A(J))*(M(J)*RN(J)-F(J)+F(J$))
     3            +P(J)-P(J$))/(M(J)-A(J)+A(J$))
     4            -E(J)+E(J$)
C
C     FACTOR OUT M
C
      DO J=1,J$
      RN1(J)=MRN1(J)/M(J)
      SN1(J)=MSN1(J)/M(J)      
      RN(J)=RN1(J)
      SN(J)=SN1(J)
      END DO
      T =T+DELTT
      ISTEP=ISTEP+1
      TOTMASS=0.
      DO J=1,J$
      TOTMASS=TOTMASS+RN1(J)
      END DO
      WRITE(6,11) T,TOTMASS
      IF (JMOD(ISTEP,IOSTEP).EQ.0) THEN
      CALL NOGRID
      CALL USXUSY(1.,40.,.FALSE.,-.5,1.,.FALSE.
     1 ,'SLOPES SCHEME TRANSPORT$')
      CALL IPLOT(X,RN1,J$,'*',.FALSE.)
      IF (ISTEP.EQ.ISTOP) STOP 50
      GO TO 1000
      ELSE
      IF (ISTEP.EQ.ISTOP) STOP 50
      GO TO 1000
      END IF
 1500 STOP 1500
   10 FORMAT(' TIME=',F10.1,'   GRID POINT=',I4,'  DENSITY=',F10.3)
   11 FORMAT(' TIME=',F10.1,'  TOTAL MASS=',F10.3)
   20 FORMAT(L5)
   25 FORMAT(' CHANGE THE TIME STEP TRUE/FALSE (F6.4)')
   26 FORMAT(F6.4)
   30 FORMAT(' DELTT/DELTX =  ',F15.5)
      END
