      SUBROUTINE INITAL
C...
C...  H2 - HEAT CONDUCTION IN ONE DIMENSIONAL SOLIDS
C...
C...  DOUBLE PRECISION CODING IS USED
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  COMMON AREA WITH THE COORDINATE SYSTEM
      COMMON/I/       NC
C...
C...  SELECT COORDINATE SYSTEM
C...
C...     NC = 0, CARTESIAN COODINATES
C...
C...     NC = 1, CYLINDRICAL COORDINATES
C...
C...     NC = 2, SPHERICAL COORDINATES
C...
      NC=0
C...
C...  INITIAL CONDITIONS
         IF(NC.EQ.0)CALL INIT10
         IF(NC.EQ.1)CALL INIT11
         IF(NC.EQ.2)CALL INIT12
      RETURN
      END
      SUBROUTINE DERV
      COMMON/I/       NC
C...
C...  CALCULATION OF TEMPORAL DERIVATIVES
         IF(NC.EQ.0)CALL DERV10
         IF(NC.EQ.1)CALL DERV11
         IF(NC.EQ.2)CALL DERV12
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/I/       NC
C...
C...  WRITING AND PLOTTING OF NUMERICAL SOLUTIONS
         IF(NC.EQ.0)CALL PRT10(NI,NO)
         IF(NC.EQ.1)CALL PRT11(NI,NO)
         IF(NC.EQ.2)CALL PRT12(NI,NO)
      RETURN
      END
      SUBROUTINE INIT10
C...
C...  INITIAL CONDITIONS FOR CARTESIAN COORDINATES
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21)
      COMMON/T/     TIME,     NSTOP,     NORUN
     +      /Y/    T(NX)
     +      /F/   TT(NX)
     +      /S/   TX(NX),   TXX(NX)
     +      /P/     HX0K
     +      /I/       NC
C...
C...  SELECT BIOT NUMBER, H*X0/K
      IF(NORUN.EQ.1)HX0K=0.0D0
      IF(NORUN.EQ.2)HX0K=1.0D0/5.0D0
      IF(NORUN.EQ.3)HX0K=1.0D0/2.0D0
      IF(NORUN.EQ.4)HX0K=1.0D0/1.0D0
      IF(NORUN.EQ.5)HX0K=1.0D0/0.5D0
      IF(NORUN.EQ.6)HX0K=1.0D0/0.05D0
C...
C...  INITIAL CONDITION
      DO 1 I=1,NX
         T(I)=1.0D0
1     CONTINUE
      RETURN
      END
      SUBROUTINE DERV10
C...
C...  TEMPORAL DERIVATIVES FOR CARTESIAN COORDINATES
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21)
      COMMON/T/     TIME,     NSTOP,     NORUN
     +      /Y/    T(NX)
     +      /F/   TT(NX)
     +      /S/   TX(NX),   TXX(NX)
     +      /P/     HX0K
     +      /I/       NC
C...
C...  BOUNDARY CONDITION AT X = 0
      NL=2
      TX(1)=0.0D0
C...
C...  BOUNDARY CONDITION AT X = X0
      IF(NORUN.NE.6)THEN
         NU=2
         TX(NX)=-HX0K*T(NX)
      ELSE IF(NORUN.EQ.6)THEN
         NU=1
         T(NX)=0.0D0
      END IF
C...
C...  TXX
      CALL DSS044(0.0D0,1.0D0,NX,T,TX,TXX,NL,NU)
C...
C...  FOURIER'S SECOND LAW
      DO 1 I=1,NX
         TT(I)=TXX(I)
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRT10(NI,NO)
C...
C...  WRITING AND PLOTTING SOLUTION FOR CARTESIAN COORDINATES
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21)
      COMMON/T/     TIME,     NSTOP,     NORUN
     +      /Y/    T(NX)
     +      /F/   TT(NX)
     +      /S/   TX(NX),   TXX(NX)
     +      /P/     HX0K
     +      /I/       NC
C...
C...  MONITOR OUTPUT
      IF(TIME.LT.0.001D0)WRITE(*,2)HX0K
2     FORMAT(//,' H*X0/K = ',F8.3,'  FROM PRT10')
      WRITE(*,1)TIME,T(1),T(NX)
C...
C...  OPEN A FILE FOR MATLAB PLOTTING
      IF((NORUN.EQ.1).AND.(TIME.LT.0.001D0))THEN
         OPEN(8,FILE='h2.out')
      END IF
C...
C...  HEADING FOR THE NUMERICAL SOLUTION
      IF(TIME.LT.0.001D0)WRITE(NO,3)
3     FORMAT(8X,'t ',3X,' x = 0 ',3X,'x = 0.2',3X,'x = 0.4',
     +               3X,'x = 0.6',3X,'x = 0.8',3X,' x = 1 ')
C...
C...  WRITE NUMERICAL SOLUTION FOR 21 GRID POINTS
         NXI=4
C...
C...     WRITE SOLUTION ON FILE OUTPUT
         WRITE(NO,1)TIME,(T(I),I=1,NX,NXI)
C...
C...     WRITE SOLUTION ON FILE H2.OUT
         WRITE( 8,1)TIME,(T(I),I=1,NX,NXI)
1        FORMAT(F10.2,6F10.5)
      RETURN
      END
      SUBROUTINE INIT11
C...
C...  INITIAL CONDITIONS FOR CYLINDRICAL COORDINATES
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NR=21)
      COMMON/T/     TIME,     NSTOP,     NORUN
     +      /Y/    T(NR)
     +      /F/   TT(NR)
     +      /S/   TR(NR),   TRR(NR)
     +      /P/     HR0K
     +      /I/       NC
C...
C...  SELECT BIOT NUMBER, H*R0/K
      IF(NORUN.EQ.1)HR0K=0.0D0
      IF(NORUN.EQ.2)HR0K=1.0D0/5.0D0
      IF(NORUN.EQ.3)HR0K=1.0D0/2.0D0
      IF(NORUN.EQ.4)HR0K=1.0D0/1.0D0
      IF(NORUN.EQ.5)HR0K=1.0D0/0.5D0
      IF(NORUN.EQ.6)HR0K=1.0D0/0.05D0
C...
C...  INITIAL CONDITION
      DO 1 I=1,NR
         T(I)=1.0D0
1     CONTINUE
      RETURN
      END
      SUBROUTINE DERV11
C...
C...  TEMPORAL DERIVATIVES FOR CYLINDRICAL COORDINATES
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NR=21)
      COMMON/T/     TIME,     NSTOP,     NORUN
     +      /Y/    T(NR)
     +      /F/   TT(NR)
     +      /S/   TR(NR),   TRR(NR)
     +      /P/     HR0K
     +      /I/       NC
C...
C...  DIRICHLET BOUNDARY CONDITION AT R = R0
      IF(NORUN.EQ.6)THEN
         NU=1
         T(NR)=0.0D0
      END IF
C...
C...  TR
      CALL DSS004(0.0,1.0D0,NR,T,TR)
C...
C...  BOUNDARY CONDITION AT R = 0
      NL=2
      TR(1)=0.0D0
C...
C...  NEUMANN BOUNDARY CONDITION AT R = R0
      IF(NORUN.NE.6)THEN
         NU=2
         TR(NR)=-HR0K*T(NR)
      END IF
C...
C...  TRR
      CALL DSS044(0.0D0,1.0D0,NR,T,TR,TRR,NL,NU)
C...
C...  FOURIER'S SECOND LAW
      DO 1 I=1,NR
         IF(I.EQ.1)THEN
            TT(I)=2.0D0*TRR(I)
         ELSE
            R=DFLOAT(I-1)/DFLOAT(NR-1)
            TT(I)=TRR(I)+(1.0D0/R)*TR(I)
         END IF
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRT11(NI,NO)
C...
C...  WRITING AND PLOTTING SOLUTION FOR CYLINDRICAL COORDINATES
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NR=21)
      COMMON/T/     TIME,     NSTOP,     NORUN
     +      /Y/    T(NR)
     +      /F/   TT(NR)
     +      /S/   TR(NR),   TRR(NR)
     +      /P/     HR0K
     +      /I/       NC
C...
C...  MONITOR OUTPUT
      IF(TIME.LT.0.001D0)WRITE(*,2)HR0K
2     FORMAT(//,' H*R0/K = ',F8.3,'  FROM PRT11')
      WRITE(*,1)TIME,T(1),T(NR)
C...
C...  OPEN A FILE FOR MATLAB PLOTTING
      IF((NORUN.EQ.1).AND.(TIME.LT.0.001D0))THEN
         OPEN(8,FILE='h2.out')
      END IF
C...
C...  HEADING FOR THE NUMERICAL SOLUTION
      IF(TIME.LT.0.001D0)WRITE(NO,3)
3     FORMAT(8X,'t ',3X,' r = 0 ',3X,'r = 0.2',3X,'r = 0.4',
     +               3X,'r = 0.6',3X,'r = 0.8',3X,' r = 1 ')
C...
C...  WRITE NUMERICAL SOLUTION FOR 21 GRID POINTS
         NRI=4
C...
C...     WRITE SOLUTION ON FILE OUTPUT
         WRITE(NO,1)TIME,(T(I),I=1,NR,NRI)
C...
C...     WRITE SOLUTION ON FILE H2.OUT
         WRITE( 8,1)TIME,(T(I),I=1,NR,NRI)
1        FORMAT(F10.2,6F10.5)
      RETURN
      END
      SUBROUTINE INIT12
C...
C...  INITIAL CONDITIONS FOR SPHERICAL COORDINATES
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NR=21)
      COMMON/T/     TIME,     NSTOP,     NORUN
     +      /Y/    T(NR)
     +      /F/   TT(NR)
     +      /S/   TR(NR),   TRR(NR)
     +      /P/     HR0K
     +      /I/       NC
C...
C...  SELECT BIOT NUMBER, H*R0/K
      IF(NORUN.EQ.1)HR0K=0.0D0
      IF(NORUN.EQ.2)HR0K=1.0D0/5.0D0
      IF(NORUN.EQ.3)HR0K=1.0D0/2.0D0
      IF(NORUN.EQ.4)HR0K=1.0D0/1.0D0
      IF(NORUN.EQ.5)HR0K=1.0D0/0.5D0
      IF(NORUN.EQ.6)HR0K=1.0D0/0.05D0
C...
C...  INITIAL CONDITION
      DO 1 I=1,NR
         T(I)=1.0D0
1     CONTINUE
      RETURN
      END
      SUBROUTINE DERV12
C...
C...  TEMPORAL DERIVATIVES FOR SPHERICAL COORDINATES
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NR=21)
      COMMON/T/     TIME,     NSTOP,     NORUN
     +      /Y/    T(NR)
     +      /F/   TT(NR)
     +      /S/   TR(NR),   TRR(NR)
     +      /P/     HR0K
     +      /I/       NC
C...
C...  DIRICHLET BOUNDARY CONDITION AT R = R0
      IF(NORUN.EQ.6)THEN
         NU=1
         T(NR)=0.0D0
      END IF
C...
C...  TR
      CALL DSS004(0.0,1.0D0,NR,T,TR)
C...
C...  BOUNDARY CONDITION AT R = 0
      NL=2
      TR(1)=0.0D0
C...
C...  NEUMANN BOUNDARY CONDITION AT R = R0
      IF(NORUN.NE.6)THEN
         NU=2
         TR(NR)=-HR0K*T(NR)
      END IF
C...
C...  TRR
      CALL DSS044(0.0D0,1.0D0,NR,T,TR,TRR,NL,NU)
C...
C...  FOURIER'S SECOND LAW
      DO 1 I=1,NR
         IF(I.EQ.1)THEN
            TT(I)=3.0D0*TRR(I)
         ELSE
            R=DFLOAT(I-1)/DFLOAT(NR-1)
            TT(I)=TRR(I)+(2.0D0/R)*TR(I)
         END IF
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRT12(NI,NO)
C...
C...  WRITING AND PLOTTING SOLUTION FOR SPHERICAL COORDINATES
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NR=21)
      COMMON/T/     TIME,     NSTOP,     NORUN
     +      /Y/    T(NR)
     +      /F/   TT(NR)
     +      /S/   TR(NR),   TRR(NR)
     +      /P/     HR0K
     +      /I/       NC
C...
C...  MONITOR OUTPUT
      IF(TIME.LT.0.001D0)WRITE(*,2)HR0K
2     FORMAT(//,' H*R0/K = ',F8.3,'  FROM PRT12')
      WRITE(*,1)TIME,T(1),T(NR)
C...
C...  OPEN A FILE FOR MATLAB PLOTTING
      IF((NORUN.EQ.1).AND.(TIME.LT.0.001D0))THEN
         OPEN(8,FILE='h2.out')
      END IF
C...
C...  HEADING FOR THE NUMERICAL SOLUTION
      IF(TIME.LT.0.001D0)WRITE(NO,3)
3     FORMAT(8X,'t ',3X,' r = 0 ',3X,'r = 0.2',3X,'r = 0.4',
     +               3X,'r = 0.6',3X,'r = 0.8',3X,' r = 1 ')
C...
C...  WRITE NUMERICAL SOLUTION FOR 21 GRID POINTS
         NRI=4
C...
C...     WRITE SOLUTION ON FILE OUTPUT
         WRITE(NO,1)TIME,(T(I),I=1,NR,NRI)
C...
C...     WRITE SOLUTION ON FILE H2.OUT
         WRITE( 8,1)TIME,(T(I),I=1,NR,NRI)
1        FORMAT(F10.2,6F10.5)
      RETURN
      END
