program fire
      include "fire.inc"
      include 'mpif.h'
      integer i_rank,i_size,ierr
      double precision f_avbn(n_prob)
      double precision f_prob(n_prob)
      double precision f_stor(n_prob)
      integer i_stat(MPI_STATUS_SIZE)

      call MPI_Init(ierr)
      call MPI_Comm_rank(MPI_COMM_WORLD,i_rank,ierr)
      call MPI_Comm_size(MPI_COMM_WORLD,i_size,ierr)
      
      call init()


c     ------------------------- loop over probabilities
      f_minp=0.0
      f_maxp=1.0
      f_delp = (f_maxp-f_minp)/(n_prob-1)
      do i_prob=1,n_prob
          f_prob(i_prob)=f_minp+(i_prob-1)*f_delp;
          f_avbn(i_prob) = 0.0;
c         --------------------- loop over trials
c         --------------------- split loop over nodes
          do i_try=1,n_try/i_size
              call burn(f_prob(i_prob))
              f_avbn(i_prob) = f_avbn(i_prob) + f_prbn()
          enddo
          f_avbn(i_prob)=f_avbn(i_prob)/n_try
      enddo

      if (i_rank.eq.0) then
          do i=1,i_size-1
              call MPI_Recv(f_stor,n_prob,MPI_DOUBLE_PRECISION,i,999,
     ~            MPI_COMM_WORLD,i_stat,ierr)
              do j=1,n_prob
                  f_avbn(j)=f_avbn(j)+f_stor(j)
              enddo
          enddo
      else
          call MPI_Send(f_avbn,n_prob,MPI_DOUBLE_PRECISION,0,999,
     ~        MPI_COMM_WORLD,ierr)
      endif
      
c     ------------------------- output results
      if (i_rank.eq.0) then
          do i_prob=1,n_prob
              write(*,*) f_prob(i_prob),f_avbn(i_prob)
          enddo
      endif
      
      call MPI_Finalize(ierr)
      end

c---------------------------------------
c     burn the forest til fire is out
c---------------------------------------       
      subroutine burn(f_prob)
      double precision f_prob
      include "fire.inc"
      
      call init_f()
 10   if (i_fout().eq.0) then
          call f_burn(f_prob)
       goto 10
      endif
     
      end
      
c---------------------------------------
c     initialize program
c---------------------------------------      
      subroutine init
      include "fire.inc"
      integer i_time(3)
      
      call itime(i_time)
      i_seed = i_time(1)+i_time(2)+i_time(3)
      call srand(i_seed)

c     set up definitions of forest state 
      i_unbt=0;
      i_smld=1;
      i_brng=2;
      i_brnt=3;
      
      end
      
c---------------------------------------
c     set up the forest
c---------------------------------------      
      subroutine init_f
      include "fire.inc"
      
c     initialize forest      
      do i=1,n_frst
          do j=1,n_frst
              i_frst(i,j)=i_unbt
          enddo
      enddo
      
c     start middle tree smoldering
      i_frst(n_frst/2,n_frst/2)=i_smld
      
      end


c---------------------------------------
c     is the fire out?
c---------------------------------------
      integer function i_fout()
      include "fire.inc"
      
      i_fout=1;
      do i=1,n_frst
          do j=1,n_frst
              if(i_frst(i,j).eq.i_smld.or.i_frst(i,j).eq.i_brng) then
                  i_fout=0;
              endif
          enddo
      enddo
      
      end

c---------------------------------------
c     calculate percentage burned
c--------------------------------------- 
      double precision function f_prbn()
      include "fire.inc"
      
      f_sum=0.0;
      f_size=n_frst*n_frst;
      do i=1,n_frst
          do j=1,n_frst
              if (i_frst(i,j).eq.i_brnt) f_sum=f_sum+1.0;
          enddo
      enddo
      f_prbn = f_sum/f_size;
      
      end

c---------------------------------------
c     burn one step through the fire
c---------------------------------------       
      subroutine f_burn(f_prob)
      double precision f_prob
      include "fire.inc"
      
      do i=1,n_frst
          do j=1,n_frst
              if (i_frst(i,j).eq.i_brng) i_frst(i,j)=i_brnt;
              if (i_frst(i,j).eq.i_smld) i_frst(i,j)=i_brng;
          enddo
      enddo
      
      do i=1,n_frst
          do j=1,n_frst
              if (i_frst(i,j).eq.i_brng) then
                  if (i.gt.1) then
                      if(rand().lt.f_prob.and.
     ~                    i_frst(i-1,j).eq.i_unbt) then
                          i_frst(i-1,j)=i_smld
                      endif
                  endif
                  if (i.lt.n_frst) then
                      if(rand().lt.f_prob.and.
     ~                    i_frst(i+1,j).eq.i_unbt) then
                          i_frst(i+1,j)=i_smld
                      endif
                  endif
                  if (j.gt.1) then
                      if(rand().lt.f_prob.and.
     ~                    i_frst(i,j-1).eq.i_unbt) then
                          i_frst(i,j-1)=i_smld
                      endif
                  endif
                  if (j.lt.n_frst) then
                      if(rand().lt.f_prob.and.
     ~                    i_frst(i,j+1).eq.i_unbt) then
                          i_frst(i,j+1)=i_smld
                      endif
                  endif
              endif
          enddo
      enddo
      
      end