module common
   public :: initial, Euler, output
   integer, parameter, public :: double = 2
   real (kind = double), public :: r,t

contains

subroutine initial(T_init,T_room,delta_t,tmax,nshow)
   real (kind = double), intent (out) :: T_init,T_room,delta_t,tmax
   integer, intent (out) :: nshow
   real (kind = double) :: tshow
   t = 0.0                   ! time
   T_init = 82.3             ! initial coffee temperature (C)
   T_room = 17.0             ! room temperature (C)
   print *, "cooling constant r ="
   read *, r
   print *, "time step dt ="
   read *, delta_t
   print *, "duration of run ="
   read *, tmax              ! minutes
   print *, "time between output of data ="
   read *, tshow
   nshow = int(tshow/delta_t)
   call output(T_init,T_room)
end subroutine initial

subroutine Euler(T_coffee,T_room,delta_t)
   real (kind = double), intent (in) :: T_room,delta_t
   real (kind = double), intent (in out) :: T_coffee
   real (kind = double) :: change
   change = -r*(T_coffee - T_room)*delta_t
   T_coffee = T_coffee + change
   t = t + delta_t
end subroutine Euler

subroutine output(T_coffee,T_room)
   real (kind = double), intent (in) :: T_coffee,T_room
   if (t == 0) then
      print *, ""
      print "(t7,a,t16,a,t28,a)", "time","T_coffee","T_coffee - T_room"
      print *, ""
   end if
   print "(f10.2,2f13.4)",t,T_coffee,T_coffee - T_room
end subroutine output

end module common

program cool
! numerical solution of Newton's law of cooling
   use common
   real (kind = double) :: T_coffee,T_room,delta_t,tmax
   integer :: nshow,counter

   call initial(T_coffee,T_room,delta_t,tmax,nshow)
   counter = 0
   do
      if (t >= tmax) then
         exit
      end if
      call Euler(T_coffee,T_room,delta_t)
      counter = counter + 1      ! number of iterations
      if (modulo(counter,nshow) == 0) then
         call output(T_coffee,T_room)
      end if
   end do
end program cool