!**********************************************************************!
!
! File: fire.f90 (21-Jun-2001)
!                (21-Jun-2005)
!
! Waldbrand-Simulation (Gittermodell fuer selbstorganisierte
! Kritizitaet), nach B. D. Malamud & D. L. Turcotte, Computing in
! Science and Engineering 2/3, pp. 42-51 (2000)
!
! nx=ny............Gitterdimension
! f_m..............Zuendfrequenz [nach je (1/f_m - 1) Baumpflanzungen
!                  wird ein Zuendholz weggeworfen]
! nsteps...........Anzahl der Zuendversuche pro Lauf
! "fire_in.dat"....Startkonfiguration fuer Lauf
! "fire_out.dat"...Endkonfiguration des Laufs
!
!**********************************************************************!

module fire_m

!**********************************************************************!

  implicit none

  public::init,finish,match,tree

! Globale Variablen

  integer,parameter,public::iin=1,iout=2

  integer,public::istep,nran,ns,nsteps,nx,ny
  integer,dimension(:),allocatable,public::ihist,iran
  integer,dimension(:,:),allocatable,public::lattice

contains

!**********************************************************************!

  subroutine init()     ! Lauf initialisieren

!**********************************************************************!

    integer::istat,mode

! Eingabeaufforderung

    do
      write(unit=*,fmt="("" Select mode=1,2,3 to"")")
      write(unit=*,fmt="(""   1...start new simulation"")")
      write(unit=*,fmt="(""   2...re-initialize simulation"")")
      write(unit=*,fmt="(""   3...continue simulation""/)")
      write(unit=*,fmt="(""   mode="")",advance="no")
      read(unit=*,fmt=*) mode
      if(1<=mode.and.mode<=3) then
        exit
      end if
    end do

    call random_seed(size=nran) ! Speicherplatz f. RNG seed
    allocate(iran(nran))

    if(mode==1) then ! Neue Simulation starten

      write(unit=*,fmt="(""  nx=ny="")",advance="no")
      read(unit=*,fmt=*) nx
      ny=nx
      do
        write(unit=*,fmt="(""  1/f_m="")",advance="no")
        read(unit=*,fmt=*) ns
        if(ns>=2) then
          exit
        end if
      end do
      write(unit=*,fmt="("" nsteps="")",advance="no")
      read(unit=*,fmt=*) nsteps

      allocate(lattice(0:nx+1,0:ny+1),ihist(1:nx*ny),stat=istat)
      if(istat/=0) then
        write(unit=*,fmt="("" fire: can't allocate"")")
        stop
      end if

      istep=0
      lattice=0
      ihist=0

    else if(mode==2) then ! Simulation re-initialisieren

      write(unit=*,fmt="("" nsteps="")",advance="no")
      read(unit=*,fmt=*) nsteps

      open(unit=iin,file="fire_in.dat",status="old",action="read", &
           form="unformatted")

      read(unit=iin) nx,istep,ns,iran
      ny=nx

      allocate(lattice(0:nx+1,0:ny+1),ihist(1:nx*ny),stat=istat)
      if(istat/=0) then
        write(unit=*,fmt="("" fire: can't allocate"")")
        stop
      end if

      read(unit=iin) lattice(1:nx,1:ny)

      close(unit=iin)

      call random_seed(put=iran) ! RNG mit seed von File initialisieren

      istep=0
      ihist=0

    else if(mode==3) then ! Simulation fortsetzen

      write(unit=*,fmt="("" nsteps="")",advance="no")
      read(unit=*,fmt=*) nsteps

      open(unit=iin,file="fire_in.dat",status="old",action="read", &
           form="unformatted")

      read(unit=iin) nx,istep,ns,iran
      ny=nx

      allocate(lattice(0:nx+1,0:ny+1),ihist(1:nx*ny),stat=istat)
      if(istat/=0) then
        write(unit=*,fmt="("" fire: can't allocate"")")
        stop
      end if

      read(unit=iin) lattice(1:nx,1:ny)
      read(unit=iin) ihist

      close(unit=iin)

      call random_seed(put=iran) ! RNG mit seed von File initialisieren

    end if

    return

  end subroutine init

!**********************************************************************!

  subroutine tree()     ! Baum an zufaelligem Gitterpunkt pflanzen

!**********************************************************************!

    integer::ix,iy
    real::ran
    
    call random_number(ran)
    ix=min(int(ran*nx+1.0),nx)
    call random_number(ran)
    iy=min(int(ran*ny+1.0),ny)

    lattice(ix,iy)=1

    return

  end subroutine tree

!**********************************************************************!

  subroutine match()     ! Zuendholz auf zufaelligen Gitterpunkt werfen

!**********************************************************************!

    integer::i,ix,iy,nburn,nn,nn_old
    integer,dimension(nx*ny)::nlistx,nlistx_old,nlisty,nlisty_old
    real::ran
    
    call random_number(ran)
    ix=min(int(ran*nx+1.0),nx)
    call random_number(ran)
    iy=min(int(ran*ny+1.0),ny)

    if(lattice(ix,iy)==1) then ! Baum auf Gitterplatz?

! Cluster von Baeumen bestimmen, die Feuer fangen

      nn=1             ! Brandherd
      nlistx(1)=ix
      nlisty(1)=iy
      lattice(ix,iy)=0

      nburn=0          ! Nachbarn u. Nachbarn der Nachbarn...
      do
        if(nn<1) then
          exit
        end if
        nburn=nburn+nn
        nn_old=nn
        nlistx_old=nlistx
        nlisty_old=nlisty
        nn=0
        do i=1,nn_old
          ix=nlistx_old(i)
          iy=nlisty_old(i)
          if(lattice(ix+1,iy)==1) then
            lattice(ix+1,iy)=0
            nn=nn+1
            nlistx(nn)=ix+1
            nlisty(nn)=iy
          end if
          if(lattice(ix-1,iy)==1) then
            lattice(ix-1,iy)=0
            nn=nn+1
            nlistx(nn)=ix-1
            nlisty(nn)=iy
          end if
          if(lattice(ix,iy+1)==1) then
            lattice(ix,iy+1)=0
            nn=nn+1
            nlistx(nn)=ix
            nlisty(nn)=iy+1
          end if
          if(lattice(ix,iy-1)==1) then
            lattice(ix,iy-1)=0
            nn=nn+1
            nlistx(nn)=ix
            nlisty(nn)=iy-1
          end if
        end do
      end do

      ihist(nburn)=ihist(nburn)+1     ! Histogramm

    end if

  end subroutine match

!**********************************************************************!

  subroutine finish()     ! Lauf abschliessen

!**********************************************************************!

    call random_seed(get=iran) ! RNG auslesen

    open(unit=iout,file="fire_out.dat",status="replace",action="write", &
         form="unformatted")

    write(unit=iout) nx,istep,ns,iran
    write(unit=iout) lattice(1:nx,1:ny)
    write(unit=iout) ihist

    close(unit=iout)

    return

  end subroutine finish

end module fire_m

!**********************************************************************!

program fire

!**********************************************************************!

  use fire_m

  implicit none

  integer::i,j

  call init()

  do i=1,nsteps
    istep=istep+1
    do j=1,ns-1
      call tree()
    end do
    call match()     ! Nach je (ns-1) Baumpflanzungen ein Zuendholz
  end do

  call finish()

end program fire

!**********************************************************************!