      SUBROUTINE INITAL
C...
C...  H3 - HEAT CONDUCTION IN A CIRCULAR FIN
C...
C...  DOUBLE PRECISION CODING IS USED
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NR=21)
      COMMON/T/      T,  NSTOP,  NORUN
     +      /Y/  U(NR)
     +      /F/ UT(NR)
     +      /P/  ALPHA,   BETA,     UA,     U0,     DR,    DRS,
     +           R(NR)
     +      /I/     IP
C...
C...  PROBLEM PARAMETERS
C...
C...     FIN INSIDE RADIUS
         R0=0.5D0
C...
C...     FIN OUTSIDE RADIUS
         R1=2.0D0
C...
C...     FIN THICKNESS
         B=0.2D0
C...
C...     HEAT TRANSFER COEFFICIENT
         H=1.0D0
C...
C...     RHO*CP
         RHOCP=2.0D0
C...
C...     THERMAL DIFFUSIVITY
         ALPHA=0.01D0
C...
C...     INITIAL TEMPERATURE
         UI=25.0D0
C...
C...     AMBIENT TEMPERATURE
         UA=25.0D0
C...
C...     INSIDE TEMPERATURE
         U0=150.0D0
C...
C...  RECIPROCAL TIME FOR HEAT TRANSFER
      BETA=(2.0D0*H)/(B*RHOCP)
C...
C...  INCREMENT FOR RADIAL GRID, SQUARE OF INCREMENT
      DR=(R1-R0)/DFLOAT(NR-1)
      DRS=DR**2
C...
C...  RADIAL GRID AND INITIAL CONDITION
      DO 1 I=1,NR
         R(I)=R0+DFLOAT(I-1)*DR
         U(I)=UI
1     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NR=21)
      COMMON/T/      T,  NSTOP,  NORUN
     +      /Y/  U(NR)
     +      /F/ UT(NR)
     +      /P/  ALPHA,   BETA,     UA,     U0,     DR,    DRS,
     +           R(NR)
     +      /I/     IP
C...
C...  BOUNDARY CONDITION AT R = R0
       U(1)=U0
      UT(1)=0.0D0
C...
C...  BOUNDARY CONDITION AT R = R1
C...
C...     RADIAL GROUP NOT INCLUDED
         IF(NORUN.EQ.1)THEN
         UT(NR)=-BETA*(U(NR)-UA)
C...
C...     RADIAL GROUP INCLUDED
         ELSE IF(NORUN.EQ.2)THEN
         UT(NR)=ALPHA*2.0D0*(U(NR-1)-U(NR))/DRS-BETA*(U(NR)-UA)
         END IF
C...
C...  INTERIOR POINTS
      DO 1 I=2,NR-1
         UT(I)=ALPHA*(U(I+1)-2.0D0*U(I)+U(I-1))/DRS
     +        +ALPHA*(1.0D0/R(I))*(U(I+1)-U(I-1))/(2.0D0*DR)
     +        -BETA*(U(NR)-UA)
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NR=21)
      COMMON/T/      T,  NSTOP,  NORUN
     +      /Y/  U(NR)
     +      /F/ UT(NR)
     +      /P/  ALPHA,   BETA,     UA,     U0,     DR,    DRS,
     +           R(NR)
     +      /I/     IP
C...
C...  MONITOR OUTPUT
      WRITE(*,*)' NORUN = ',NORUN,'  T = ',T
C...
C...  IF THE SOLUTION HAS A PHYSICALLY UNREALISTIC VALUE, TERMINATE
C...  EXECUTION OF CURRENT RUN
      DO 2 I=1,NR
         IF((U(I).GT.U0).OR.(U(I).LT.UA))THEN
         WRITE(NO,1)NORUN,T
1        FORMAT(//,' RUN NO = ',I2,' AT T = ',F7.2,' TERMINATED',//)
         NSTOP=1
         END IF
2     CONTINUE
C...
C...  PRINT THE SOLUTION
      WRITE(NO,3)T,(R(I),I=1,NR,4),(U(I),I=1,NR,4)
3     FORMAT(' t =   ',F6.2,/,'      r ',6F9.2,/,' u(r,t) ',6F9.3,/)
      RETURN
      END
