*DECK HEAT16
      SUBROUTINE INITAL
C...
C...  THERMAL BONDING OF SOLIDS
C...
C...  TWO SHEETS OF A SOLID MATERIAL ARE TO BE BONDED TOGETHER BY
C...  USING A THIN LAYER OF THERMOPLASTIC MATERIAL WHICH FUSES AND
C...  FORMS A BOND AT TEMPERATURE UG.  THE TWO PLATES OF SOLID
C...  MATERIAL ARE INSERTED IN A PRESS, WITH BOTH PLATTENS OF THE
C...  PRESS HEATED AT A CONSTANT RATE, Q.  THE THERMOPLASTIC MATER-
C...  IAL BETWEEN THE PLATES THEREFORE IS HEATED, AND THE FOLLOWING
C...  ANALYSIS IS INTENDED TO DETERMINE THE TEMPERATURE OF THE
C...  MATERIAL, U(O,T), AS A FUNCTION OF TIME, T.  THE TIME REQUIRED
C...  TO REACH THE CRITICAL TEMPERATURE, UG, CAN THEREFORE BE DETER-
C...  MINED SO THAT THAT TIME REQUIRED IN THE PRESS TO FORM THE BOND
C...  CAN BE ESTIMATED.
C...
C...  THE FOLLOWING MATHEMATICAL MODEL IS PROGRAMMED
C...
C...     U  = A*U
C...      T      XX
C...
C...     U (0,T) = 0
C...      X
C...
C...     K*U (XL,T) = -Q
C...        X
C...
C...     U(X,0) = U0
C...
C...  WHERE
C...
C...     U       TEMPERATURE IN EITHER PLATE (C)
C...
C...     T       TIME (SEC)
C...
C...     X       DISTANCE WITHIN EITHER PLATE (X = 0 CORRESPONDS TO
C...             THE SURFACE FOR THE THERMOPLASTIC BONDING MATERIAL)
C...
C...     A       THERMAL DIFFUSIVITY OF THE PLATES (CM**2/SEC)
C...
C...     K       THERMAL CONDUCTIVITY OF THE PLATES (J/SEC-CM**2-C)
C...
C...     Q       HEAT FLUX APPLIED BY THE PRESS TO THE PLATES
C...             (J/SEC-CM**2)
C...
C...     XL      THICKNESS OF THE PLATES (CM)
C...
C...     U0      INITIAL TEMPERATURE OF THE PLATES (C)
C...
C...  THE FIRST EQUATION IS JUST FOURIER*S SECOND LAW FOR HEAT CON-
C...  DUCTION IN THE PLATES WITH THERMAL DIFFUSIVITY A.  THE SECOND
C...  EQUATION IS A BOUNDARY CONDITION REFLECTING SYMMETRY WHERE THE
C...  TWO PLATES COME TOGETHER (AND WHERE THE THERMOPLASTIC MATERIAL
C...  IS LOCATED).  THE THIRD EQUATION IS A BOUNDARY CONDITION AT
C...  X = XL SPECIFYING THE HEAT FLUX Q AT THE INTERFACE BETWEEN THE
C...  PRESS PLATTENS AND THE PLATES (THIS IS FOURIER*S FIRST LAW).
C...  THE FOURTH EQUATION IS AN INITIAL CONDITION INDICATING THAT THE
C...  PLATES START AT TEMPERATURE U0.
C...
C...  THIS MODEL IS PROGRAMMED BELOW FOR A SPATIAL GRID OF 11 AND
C...  21 POINTS, USING THREE AND FIVE-POINT APPROXIMATIONS OF THE
C...  DERIVATIVES IN X (SUBROUTINES DSS002 AND DSS004).  THUS, FOUR
C...  RUNS ARE PROGRAMMED.  ALSO, IN THE FOLLOWING PROGRAMMING, B =
C...  Q/K.
C...
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     U(21)
     2      /F/    UT(21)
     3      /S/    UX(21),   UXX(21)
     4      /C/         A,         B,        XL,         N,        IP
C...
C...  SET THE LENGTH OF THE SOLID, NUMBER OF GRID POINTS
      XL=2.
      IF((NORUN.EQ.1).OR.(NORUN.EQ.3))N=11
      IF((NORUN.EQ.2).OR.(NORUN.EQ.4))N=21
C...
C...  SET THE MODEL CONSTANTS
      A=0.009
      B=-1000.
C...
C...  SET THE MODEL INITIAL CONDITION
      DO 1 I=1,N
      U(I)=20.
1     CONTINUE
C...
C...  INITIALIZE THE CALCULATIONS IN SUBROUTINE DERV
      CALL DERV
C...
C...  INITIALIZE THE COUNTER FOR THE PLOTTED SOLUTION
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     U(21)
     2      /F/    UT(21)
     3      /S/    UX(21),   UXX(21)
     4      /C/         A,         B,        XL,         N,        IP
C...
C...  COMPUTE THE DERIVATIVE UX
      IF((NORUN.EQ.1).OR.(NORUN.EQ.2))CALL DSS002(0.,XL,N,U,UX)
      IF((NORUN.EQ.3).OR.(NORUN.EQ.4))CALL DSS004(0.,XL,N,U,UX)
C...
C...  APPLY THE BOUNDARY CONDITIONS AS CONSTRAINTS
      UX(1)=0.
      UX(N)=-B
C...
C...  COMPUTE THE DERIVATIVE UXX
      IF((NORUN.EQ.1).OR.(NORUN.EQ.2))CALL DSS002(0.,XL,N,UX,UXX)
      IF((NORUN.EQ.3).OR.(NORUN.EQ.4))CALL DSS004(0.,XL,N,UX,UXX)
C...
C...  ASSEMBLE THE PDE
      DO 1 I=1,N
      UT(I)=A*UXX(I)
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     U(21)
     2      /F/    UT(21)
     3      /S/    UX(21),   UXX(21)
     4      /C/         A,         B,        XL,         N,        IP
C...
C...  DIMENSION THE ARRAYS TO STORE THE PLOTTED SOLUTION
      DIMENSION   TP(101),   UP(101)
C...
C...  PRINT A HEADING FOR THE NUMERICAL SOLUTION
      IF(IP.EQ.0)WRITE(NO,2)
2     FORMAT(1H ,//,9X,1HT,4X,6HU(0,T),3X,7HU(XL,T),
     1                     3X,7HUX(0,T),2X,8HUX(XL,T))
C...
C...  PRINT THE NUMERICAL SOLUTION EVERY FIFTH CALL TO PRINT
      IP=IP+1
      IF(((IP-1)/5*5).EQ.(IP-1))WRITE(NO,1)T,U(1),U(N),UX(1),UX(N)
1     FORMAT(F10.1,4F10.2)
C...
C...  STORE THE NUMERICAL SOLUTION FOR PLOTTING
      TP(IP)=T
      UP(IP)=U(1)
C...
C...  PLOT THE SOLUTION AT THE END OF THE RUN
      IF(IP.LT.101)RETURN
      CALL TPLOTS(1,IP,TP,UP)
      WRITE(NO,3)
3     FORMAT(1H ,//,22H U(0,T) (C) VS T (SEC))
      RETURN
      END

