SUBROUTINE INITAL C... C... M3 - NONLINEAR, FRONT-SHARPENING, CONVECTIVE FLOW C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/T/ T, NSTOP, NORUN + /I/ N, NCASE, IP, NT C... C... SELECT CASE NCASE=3 C... C... SELECT INITIALIZATION ROUTINE IF(NCASE.EQ.1)CALL INIT1 IF(NCASE.EQ.2)CALL INIT2 IF(NCASE.EQ.3)CALL INIT3 RETURN END SUBROUTINE INIT1 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NX=101) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NX) + /F/ UT(NX) + /S/ UX(NX) + /X/ X(NX) + /C/ A, B, C, XL, DX + /I/ N, NCASE, IP, NT C... C... SELECT THE SOLUTION PARAMETERS C... C... FRONT FLATTENING IF(NORUN.EQ.1)THEN A=1.0D0 B=1.0D0 C=0.0D0 NT=11 C... C... FRONT SHARPENING ELSE + IF(NORUN.EQ.2)THEN A= 1.01D0 B=-1.0D0 C= 1.0D0 NT=11 END IF C... C... NUMBER OF GRID POINTS N=11 C... C... TOTAL LENGTH IN X, GRID SPACING XL=1.0D0 DX=XL/DFLOAT(N-1) C... C... INITIAL CONDITION DO 1 I=1,N X(I)=XL*DFLOAT(I-1)/DFLOAT(N-1) U(I)=ANAL(X(I),0.0D0) 1 CONTINUE C... C... INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE INIT2 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NX=101) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NX) + /F/ UT(NX) + /S/ UX(NX) + /X/ X(NX) + /C/ A, B, C, XL, DX + /I/ N, NCASE, IP, NT C... C... SELECT THE SOLUTION PARAMETERS C... C... FRONT FLATTENING IF(NORUN.EQ.1)THEN A=1.0D0 B=1.0D0 C=0.0D0 NT=11 C... C... FRONT SHARPENING ELSE + IF(NORUN.EQ.2)THEN A= 1.0D0 B=-1.0D0 C= 1.0D0 NT=13 END IF C... C... NUMBER OF GRID POINTS N=101 C... C... TOTAL LENGTH IN X, GRID SPACING XL=1.0D0 DX=XL/DFLOAT(N-1) C... C... INITIAL CONDITION DO 1 I=1,N X(I)=XL*DFLOAT(I-1)/DFLOAT(N-1) U(I)=ANAL(X(I),0.0D0) 1 CONTINUE C... C... INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE INIT3 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NX=101) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NX) + /F/ UT(NX) + /S/ UX(NX) + /X/ X(NX) + /C/ A, B, C, XL, DX + /I/ N, NCASE, IP, NT C... C... SELECT THE SOLUTION PARAMETERS C... C... FRONT FLATTENING IF(NORUN.EQ.1)THEN A=1.0D0 B=1.0D0 C=0.0D0 NT=11 C... C... FRONT SHARPENING ELSE + IF(NORUN.EQ.2)THEN A= 1.0D0 B=-1.0D0 C= 1.0D0 NT=11 END IF C... C... NUMBER OF GRID POINTS N=101 C... C... TOTAL LENGTH IN X, GRID SPACING XL=1.0D0 DX=XL/DFLOAT(N-1) C... C... INITIAL CONDITION DO 1 I=1,N X(I)=XL*DFLOAT(I-1)/DFLOAT(N-1) U(I)=ANAL(X(I),0.0D0) 1 CONTINUE C... C... INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE DERV IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/T/ T, NSTOP, NORUN + /I/ N, NCASE, IP, NT C... C... SELECT DERIVATIVE SUROUTINE IF(NCASE.EQ.1)CALL DERV1 IF(NCASE.EQ.2)CALL DERV2 IF(NCASE.EQ.3)CALL DERV3 RETURN END SUBROUTINE DERV1 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NX=101) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NX) + /F/ UT(NX) + /S/ UX(NX) + /X/ X(NX) + /C/ A, B, C, XL, DX + /I/ N, NCASE, IP, NT C... C... BOUNDARY CONDITION AT X = 0 U(1)=ANAL(0.0D0,T) UT(1)=0.0D0 C... C... PDE DO 1 I=2,N UT(I)=-U(I)*(U(I)-U(I-1))/DX 1 CONTINUE RETURN END SUBROUTINE DERV2 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NX=101) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NX) + /F/ UT(NX) + /S/ UX(NX) + /X/ X(NX) + /C/ A, B, C, XL, DX + /I/ N, NCASE, IP, NT C... C... BOUNDARY CONDITION AT X = 0 U(1)=ANAL(0.0D0,0.0D0) UT(1)=0.0D0 C... C... PDE DO 1 I=2,N UT(I)=-U(I)*(U(I)-U(I-1))/DX 1 CONTINUE RETURN END SUBROUTINE DERV3 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NX=101) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NX) + /F/ UT(NX) + /S/ UX(NX) + /X/ X(NX) + /C/ A, B, C, XL, DX + /I/ N, NCASE, IP, NT C... C... BOUNDARY CONDITION AT X = 0 U(1)=ANAL(0.0D0,0.0D0) UT(1)=0.0D0 C... C... PDE C... X NE 0 OR XL DO 1 I=2,N-1 UT(I)=-U(I)*(U(I+1)-U(I-1))/(2.0D0*DX) 1 CONTINUE C... C... X = XL UT(N)=-U(N)*(U(N)-U(N-1))/DX RETURN END SUBROUTINE PRINT(NI,NO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/T/ T, NSTOP, NORUN + /I/ N, NCASE, IP, NT C... C... SELECT PRINT ROUTINE IF(NCASE.EQ.1)CALL PRINT1(NI,NO) IF(NCASE.EQ.2)CALL PRINT2(NI,NO) IF(NCASE.EQ.3)CALL PRINT3(NI,NO) RETURN END SUBROUTINE PRINT1(NI,NO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NX=101) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NX) + /F/ UT(NX) + /S/ UX(NX) + /X/ X(NX) + /C/ A, B, C, XL, DX + /I/ N, NCASE, IP, NT C... C... MONITOR OUTPUT ON SCREEN WRITE(*,*)NORUN,T C... C... PRINT NUMERICAL AND ANALYTICAL SOLUTIONS WRITE(NO,3)T 3 FORMAT(/,' t = ',F5.2,/,9X,'x',4X,'u(x,t)',3X,'ua(x,t)') DO 1 I=1,N UA=ANAL(X(I),T) WRITE(NO,2)X(I),U(I),UA 2 FORMAT(F10.2,2F10.5) 1 CONTINUE C... C... PLOT THE SOLUTION IP=IP+1 CALL PLOTM RETURN END SUBROUTINE PRINT2(NI,NO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NX=101) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NX) + /F/ UT(NX) + /S/ UX(NX) + /X/ X(NX) + /C/ A, B, C, XL, DX + /I/ N, NCASE, IP, NT C... C... MONITOR OUTPUT ON SCREEN WRITE(*,*)NORUN,T C... C... PRINT NUMERICAL SOLUTION WRITE(NO,3)T 3 FORMAT(/,' t = ',F5.2,/,9X,'x',4X,'u(x,t)') DO 1 I=1,N WRITE(NO,2)X(I),U(I) 2 FORMAT(F10.2,F10.5) 1 CONTINUE C... C... PLOT THE SOLUTION IP=IP+1 CALL PLOTM RETURN END SUBROUTINE PRINT3(NI,NO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NX=101) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NX) + /F/ UT(NX) + /S/ UX(NX) + /X/ X(NX) + /C/ A, B, C, XL, DX + /I/ N, NCASE, IP, NT C... C... MONITOR OUTPUT ON SCREEN WRITE(*,*)NORUN,T C... C... PRINT NUMERICAL SOLUTION WRITE(NO,3)T 3 FORMAT(/,' t = ',F5.2,/,9X,'x',4X,'u(x,t)') DO 1 I=1,N WRITE(NO,2)X(I),U(I) 2 FORMAT(F10.2,F10.5) 1 CONTINUE C... C... PLOT THE SOLUTION IP=IP+1 CALL PLOTM RETURN END DOUBLE PRECISION FUNCTION ANAL(X,T) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/C/ A, B, C, XL, DX ANAL=(C+B*X)/(A+B*T) RETURN END SUBROUTINE PLOTM IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(NX=101) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NX) + /F/ UT(NX) + /S/ UX(NX) + /X/ X(NX) + /C/ A, B, C, XL, DX + /I/ N, NCASE, IP, NT C... C... OPEN A FILE FOR MATLAB PLOTTING IF((IP.EQ.1).AND.(NORUN.EQ.1))THEN OPEN(1,FILE='m3.out') END IF C... C... WRITE THE NUMERICAL SOLUTION FOR SUBSEQUENT MATLAB PLOTTING DO 2 I=1,N WRITE(1,3)X(I),U(I) 3 FORMAT(F10.3,F10.5) 2 CONTINUE RETURN END