      SUBROUTINE INITAL
C...
C...  TEMPERATURE DISTRIBUTION IN A NUCLEAR FUEL ROD ASSEMBLY
C...
C...  THE FOLLOWING MODEL HAS BEEN PROPOSED FOR THE TEMPERATURE
C...  DISTRIBUTION IN A NUCLEAR FUEL ROD ASSEMBLY
C...
C...  TF  = (KF/(RHOF*CPF))*(TF   + (1/R)*TF ) - (RD*HR/(RHOF*CPF)) (1)
C...    T                      RR           R
C...
C...  TF(R,0) = TF0, TF (0,T) = 0, TF(RF,T) = TA(RF,T)        (2)(3)(4)
C...                   R
C...
C...  TA  = (KA/(RHOA*CPA))*(TA   + (1/R)*TA )                      (5)
C...    T                      RR           R
C...
C...  TA(R,0) = TA0, TA(RA,T) = TC, -KA*TA (RF,T) = -KF*TF (RF,T)
C...                                      R               R
C...                                                          (6)(7)(8)
C...  WHERE
C...
C...  TF(R,T)  FUEL TEMPERATURE, 0 LE R LE RF
C...
C...  TA(R,T)  ALUMINUM CLADDING TEMPERATURE, RF LE R LE RA
C...
C...  R        RADIAL POSITION IN THE ASSEMBLY
C...
C...  T        TIME
C...
C...  TF0      INITIAL FUEL TEMPERATURE
C...
C...  TA0      INITIAL ALUMINUM TEMPERATURE
C...
C...  TC       COOLANT TEMPERATURE
C...
C...  RHOF,    DENSITY, HEAT CAPACITY AND THERMAL CONDUCTIVITY OF
C...  CPF,KF   THE FUEL
C...
C...  RHOA,    DENSITY, HEAT CAPACITY AND THERMAL CONDUCTIVITY OF
C...  CPA,KA   THE ALUMINUM CLADDING
C...
C...  RD,HR    VOLUMETRIC DECAY RATE AND HEAT OF REACTION OF THE
C...           NUCLEAR FUEL
C...
C...  DEVELOP A COMPUTER SIMULATION OF THIS SYSTEM WHICH PRINTS AND
C...  AND PLOTS TF(0,T) AND TA(RF,T) VS T SO THAT THE HEATING OF THE
C...  FUEL AND CLADDING CAN BE DETERMINED.  THE FOLLOWING PARAMETERS
C...  APPLY
C...
C...  RHOF = 0.5, CPF = 0.8, KF = 1.0, RHOA = 1.1, CPA = 0.7, KA = 10,
C...
C...  RD = 0.1, HR = -10000, RF = 2.5, RA = 5.0, TC = 85, TF0 = 25,
C...
C...  TA0 = 25
C...
      COMMON/T/         T,      NFIN,     NORUN
     1      /Y/    TF(11),    TA(11)
     2      /F/   TFT(11),   TAT(11)
     3      /S/   TFR(11),  TFRR(11),   TAR(11),  TARR(11)
     4      /C/      RFOF,       CPF,        KF,      RHOA,       CPA,
     5                 KA,        RD,        HR,        RF,        RA,
     6                 TC,       TF0,       TA0,    R1(11),    R2(11),
     7                CF1,       CF2,       CA1,        NF,        NA,
     8                 IP
      REAL             KF,        KA
C...
C...  SET THE MODEL PARAMETERS
      RHOF=0.5
      CPF=0.8
      KF=1.0
      RHOA=1.1
      CPA=0.7
      KA=10.
      RD=0.1
      HR=-10000.
      RF=2.5
      RA=5.0
      TC=85.
      TF0=25.
      TA0=25.
C...
C...  PRECOMPUTE SOME CONSTANTS FOR USE IN SUBROUTINE DERV
      CF1=KF/(RHOF*CPF)
      CF2=RD*HR/(RHOF*CPF)
      CA1=KA/(RHOA*CPA)
C...
C...  INITIAL CONDITIONS AND RADIAL GRID FOR THE FUEL
      NF=11
      DO 1 I=1,NF
      TF(I)=TF0
      R1(I)=RF*FLOAT(I-1)/FLOAT(NF-1)
1     CONTINUE
C...
C...  INITIAL CONDITIONS AND RADIAL GRID FOR THE CLADDING
      NA=11
      DO 2 I=1,NA
      TA(I)=TA0
      R2(I)=RF+(RA-RF)*FLOAT(I-1)/FLOAT(NA-1)
2     CONTINUE
C...
C...  COMPUTE THE INITIAL DERIVATIVES
      CALL DERV
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/         T,      NFIN,     NORUN
     1      /Y/    TF(11),    TA(11)
     2      /F/   TFT(11),   TAT(11)
     3      /S/   TFR(11),  TFRR(11),   TAR(11),  TARR(11)
     4      /C/      RFOF,       CPF,        KF,      RHOA,       CPA,
     5                 KA,        RD,        HR,        RF,        RA,
     6                 TC,       TF0,       TA0,    R1(11),    R2(11),
     7                CF1,       CF2,       CA1,        NF,        NA,
     8                 IP
      REAL             KF,        KA
C...
C...  COMPUTE THE SPATIAL DERIVATIVES IN TF
      TF(NF)=TA(1)
      CALL DSS004(0.,RF,NF,TF ,TFR )
      TFR(1)=0.
      CALL DSS004(0.,RF,NF,TFR,TFRR)
C...
C...  COMPUTE THE SPATIAL DERIVATIVES IN TA
      TA(NA)=TC
      CALL DSS004(RF,RA,NA,TA ,TAR )
      TAR(1)=(KF/KA)*TFR(NF)
      CALL DSS004(RF,RA,NA,TAR,TARR)
C...
C...  PDE FOR TF
      TFT(1)=CF1*2.*TFRR(1)-CF2
      DO 1 I=2,NF
      TFT(I)=CF1*(TFRR(I)+(1./R1(I))*TFR(I))-CF2
1     CONTINUE
C...
C...  PDE FOR TA
      DO 2 I=1,NA
      TAT(I)=CA1*(TARR(I)+(1./R2(I))*TAR(I))
2     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/         T,      NFIN,     NORUN
     1      /Y/    TF(11),    TA(11)
     2      /F/   TFT(11),   TAT(11)
     3      /S/   TFR(11),  TFRR(11),   TAR(11),  TARR(11)
     4      /C/      RFOF,       CPF,        KF,      RHOA,       CPA,
     5                 KA,        RD,        HR,        RF,        RA,
     6                 TC,       TF0,       TA0,    R1(11),    R2(11),
     7                CF1,       CF2,       CA1,        NF,        NA,
     8                 IP
      REAL             KF,        KA
C...
C...  DIMENSION THE ARRAYS FOR THE PLOTTED SOLUTION
      DIMENSION     TP(31),   TFP(31),   TAP(31),
     1                     TFPR(7,31),TAPR(7,31)
C...
C...  INITIALIZE COUNTERS FOR THE PLOTTED SOLUTION
      DATA IP,IR/0,0/
C...
C...  PRINT THE NUMERICAL SOLUTION
      IF((IP/5*5).NE.IP)GO TO 10
      NS=2
      WRITE(NO,1)T,
     1 (TF(I),I=1,NF,NS),(TA(I),I=1,NA,NS)
1     FORMAT(1H ,//,5H  T  ,F9.1,/,
     1              5H TF  ,6F10.1,/,
     2              5H TA  ,6F10.1)
C...
C...  PRINT THE TIME DERIVATIVES
      IF(IP.GT.0)GO TO 10
      WRITE(NO,2)(TFT(I),I=1,NF,NS),(TAT(I),I=1,NA,NS)
2     FORMAT(1H ,/, 5H TFT ,6F10.1,/,
     1              5H TAT ,6F10.1)
C...
C...  PRINT THE SPATIAL DERIVATIVES
      WRITE(NO,3)( TFR(I),I=1,NF,NS),( TAR(I),I=1,NA,NS),
     1           (TFRR(I),I=1,NF,NS),(TARR(I),I=1,NA,NS)
3     FORMAT(1H ,/, 5H TFR ,6F10.1,/,
     1              5H TAR ,6F10.1,/,
     2              5H TFRR,6F10.1,/,
     3              5H TARR,6F10.1)
C...
C...  STORE TF(0,T) AND TA(RF,T) FOR PLOTTING
10    IP=IP+1
      TP(IP)=T
      TFP(IP)=TF(1)
      TAP(IP)=TA(1)
      IF((IP-1)/5*5.NE.(IP-1))RETURN
C...
C...  STORE TF(R,T) AND TA(R,T) FOR PLOTTING OF THE TEMPERATURE
C...  PROFILES
      IR=IR+1
      DO 4 I=1,NF
      TFPR(IR,I)=TF(I)
4     CONTINUE
      DO 5 I=1,NA
      TAPR(IR,I)=TA(I)
5     CONTINUE
C...
C...  PLOT THE NUMERICAL SOLUTION
      IF(IP.LT.31)RETURN
      CALL TPLOTS(1,IP,TP,TFP)
      WRITE(NO,11)
11    FORMAT(1H ,//,41H TF(0,T) VS T                            )
      CALL TPLOTS(1,IP,TP,TAP)
      WRITE(NO,12)
12    FORMAT(1H ,//,41H TA(RF,T) VS T                           )
      CALL TPLOTS(IR,NF,R1,TFPR)
      WRITE(NO,13)
13    FORMAT(1H ,//,41H TF(R,T) VS R, T = 0, 0.5, 1, ..., 3     )
      CALL TPLOTS(IR,NA,R2,TAPR)
      WRITE(NO,14)
14    FORMAT(1H ,//,41H TA(R,T) VS R, T = 0, 0.5, 1, ..., 3     )
      IP=0
      IR=0
      RETURN
      END

