      PROGRAM cool
*     numerical solution of Newton's law of cooling
      IMPLICIT NONE
      REAL*8 T_coffee, T_room, delta_t, r, t, tmax
      INTEGER counter, nshow
      CALL initial(t,T_coffee,T_room,r,delta_t,tmax,nshow)
      counter = 0                           ! initialize counter for clarity
      DO WHILE(t .LE. tmax)
         CALL Euler(t,T_coffee,T_room,r,delta_t)
         counter = counter + 1              ! number of iterations
         IF(mod(counter,nshow) .EQ. 0) THEN
            CALL output(t,T_coffee,T_room)
         END IF
      END DO
      END
C
      SUBROUTINE initial(t,T_init,T_room,r,delta_t,tmax,nshow)
*     time is keyword in True BASIC so cannot be used as a variable
      IMPLICIT NONE
      REAL*8 tshow, t, T_init, T_room, r, delta_t, tmax
      INTEGER nshow
      t = 0                                 ! time
      T_init = 82.3D0                       ! initial coffee temperature (C)
      T_room = 17                           ! room temperature (C)
      WRITE(*,'(A,$)') 'cooling constant r = '
      READ(*,*) r
      WRITE(*,'(A,$)') 'time step dt = '
      READ(*,*) delta_t
      WRITE(*,'(A,$)') 'duration of run = '
      READ(*,*) tmax                        ! minutes
      WRITE(*,'(A,$)') 'time between output of data = '
      READ(*,*) tshow
      nshow = nint(tshow/delta_t)
      CALL output(t,T_init,T_room)
      END
C
      SUBROUTINE Euler(t,T_coffee,T_room,r,delta_t)
      IMPLICIT NONE
      REAL*8 change, t, T_coffee, T_room, r, delta_t
      change = -r*(T_coffee - T_room)*delta_t
      T_coffee = T_coffee + change
      t = t + delta_t
      END
C
      SUBROUTINE output(t,T_coffee,T_room)
      IMPLICIT NONE
      REAL*8 t, T_coffee, T_room
      IF(t .EQ. 0) THEN
         WRITE(*,*)
         WRITE(*,*) '  time (min)    T_coffee  T_coffee - T_room'
         WRITE(*,*)
      END IF
      WRITE(*,'(3F13.6)') t, T_coffee, T_coffee - T_room 
                                            ! show results on screen
      END
