!*****************************************************
!*  dr. dale a. hudson         office: 104 pond
!*  assist. prof: penn state     phone: 863-1945
!*    fortran 90/95, c++       dhudson@cse.psu.edu
!*    engineering mechanics:dynamics
!*****************************************************
!* comments and formatting added by bob harris
!*****************************************************

!     last change:  dah   5 dec 2003   10:35 am
program hotmetal
implicit none

integer, parameter :: maxx = 100
integer, parameter :: maxy = maxx
real (kind=selected_real_kind(14,300)), dimension(maxx, maxy) :: f, df
real :: sum_sq_error, pi=3.14157265
logical :: not_at_edge(maxx, maxy)
integer :: i, j, iteration

!---------
! open files for results
!---------

open (10, file="hotmetal.dat",   status="unknown", form="formatted")
open (12, file="hotmetal.stats", status="unknown", form="formatted")

!---------
! create initial temperature pattern on the plate
!
! note that this program treats the 2D indices as <column,row> instead of
! <row,column>.  <column,row> matches the internal ordering in which fortran
! stores the cells in memory.  It is similar to cartesian <x,y> coordinates
! except that y increases from top to bottom.
!---------

f = 0.

! left column, x = 1 boundary

f(1,:) = (/ (400.+(970.-400.)*(real(j)/maxy),j=1,maxy) /)

! right column, x = maxx boundary

f(maxx,:) = (/ (400.+(1090.-400.)*(real(j)/maxy),j=1,100) /)

! top row, y = 1 boundary

f(1:50,1) = (/ (400.+(700.-400.)*(real(i)/50),i=1,50) /)
do i=51,100
    f(i,1) = 700.+(400.-700.)*(real(i-50)/50)
enddo

! bottom row, y = maxy boundary

f(:,maxy) = (/ (800*sin(pi*real(i+1)/100)+1000,i=1,100) /)

!---------
! create a mask array that's true in the interior and false around the
! boundary
!---------

not_at_edge             = .false.
not_at_edge(2:99, 2:99) = .true.

!---------
! iterate until temperature pattern reaches steady state
!---------

iteration = 0

do
    iteration = iteration + 1

    ! compute the amount of change at each location where the mask is true,
    ! and apply the change to update the temperature pattern on the plate

    df = 0.
    where (not_at_edge)
        df = 0.25*(cshift(f,+1,1) + cshift(f,-1,1) &
                 + cshift(f,+1,2) + cshift(f,-1,2)) - f
        f = f + df
    endwhere

    ! figure out our measure of how much things are changing

    sum_sq_error = sum(df*df)

    ! every 100th iteration, report how we're doing 

    if (mod(iteration, 100) .eq. 0) then
        write(12,*) iteration, sum_sq_error
        print *, iteration, sum_sq_error
    endif

    ! if things aren't changing much, we're done

    if (sum_sq_error .lt. 1.e-12) exit
end do

!---------
! write results to files, so we can plot 'em in matlab
!---------

write(12,*) iteration, sum_sq_error

write(10,100) ((f(i,j),j=1,maxy),i=1,maxx)
100 format( 100e15.5 )

end program hotmetal







