# Flood fill for phase field
# This is a FORTRAN code to identify region outlined by a phase-field variable

module flood_fill
!
 implicit none
!
  contains
!
!-------------------------------------------------------------------------------------

! This is a 2D code, assume the system size is
! System_Size(1)xSystem_Size(2)xSystem_Size(3)
!

Subroutine flood_fill_PSD(yita, area)
  implicit none
  integer, parameter :: MAXPN=9999, System_Size(3)=(/1, 128, 128/)
  real :: System_Volume=real(System_Size(1)*System_Size(2)*System_Size(3))
  integer i,j,k, ka, kb, jr, jl, nprep, SSeedA(MAXPN,50,3), SSeedB(MAXPN,50,3), use_stack(MAXPN), Sstack(MAXPN, 2)
  integer :: js,jt, jtt, landing_P(3)
  logical :: safety_valve
  real :: boundvalue=0.5
  real :: area(MAXPN), yita(System_Size(1), System_Size(2), System_Size(3))
  integer :: Is_scanning_P(MAXPN), counted(System_Size(1), System_Size(2), System_Size(3))

  SSeedA=0; SSeedB=0; Sstack=1; use_stack=0
  area=0; counted=0
  nprep=1; landing_P=0
  i=1; j=1; k=1
  safety_valve=.TRUE.

  do while (k.le.System_Size(3).and.safety_valve)

     if(yita(i,j,k).ge.boundvalue.and.counted(i,j,k).eq.0)then
        jr=j
     else
        jr=j+1
        if(jr.eq.System_Size(2)+1)jr=1
        do while ((yita(i,jr,k).lt.boundvalue.or.counted(i,jr,k).eq.1).and.jr.ne.j)
           jr=jr+1
           if(jr.eq.System_Size(2)+1) jr=1
        enddo
     endif
   
     If(jr.eq.j.and.(yita(i,j,k).lt.boundvalue.or.counted(i,j,k).eq.1))then  ! this means no uncounted particle in this line
        k=k+1; j=1 ! if no particle in this line, move up one line and continue scanning
     else ! else there is a uncounted precipitate now 

!        write(12,*)"found a seed at ", jr, k
        Is_scanning_P(nprep)=1
        jl=jr;
        do while (yita(i,jl,k).ge.boundvalue.and.counted(i,jl,k).ne.1)
           jl=jl-1
           if(jl.eq.0)jl=System_Size(2)
        enddo
        jl=jl+1
        if(jl.eq.System_Size(2)+1)jl=1
        
        do while (yita(i,jr,k).ge.boundvalue.and.counted(i,jr,k).ne.1)
           jr=jr+1
           if(jr.eq.System_Size(2)+1)jr=1
        enddo
        jr=jr-1
        if(jr.eq.0) jr=System_Size(2)
                ! Now jl and jr mark the left and right bound of particle line

        if(jr.ge.jl)then
           counted(i, jl:jr,k)=1
           area(nprep)=area(nprep)+ jr-jl+1
        else
           counted(i, 1:jr,k)=1; counted(i, jl: System_Size(2),k)=1
           area(nprep)=area(nprep)+jr; area(nprep)=area(nprep)+System_Size(2)-jl+1
        endif
!        write(12,*)"integrate line", k," jl:jr ", jl, jr, "area:", area(nprep)

        if(use_stack(nprep).eq.0)then
            landing_P(2)=jr+1; landing_P(3)=k
            if(jr+1.eq.System_Size(2)+1)then
                landing_P(2)=1; landing_P(3)=k+1
            endif
        endif

        ka=k+1
        if(ka.eq.System_Size(3)+1) ka=1
        js=jl
        do while(js.ne.jr.and.(yita(i,js,ka).lt.boundvalue.or.counted(i,js,ka).eq.1))
           js=js+1
           if(js.eq.System_Size(2)+1)js=1
        enddo
        if((js.ne.jr.and.yita(i,js,ka).ge.boundvalue).and.counted(i,js,ka).eq.0) then
           SSeedA(nprep,Sstack(nprep,1),2)= js; SSeedA(nprep,Sstack(nprep,1),3)= ka
           Sstack(nprep,1)=Sstack(nprep,1)+1
           jt=js+1 
           if(jt.eq.System_Size(2)+1)jt=1
           jtt=jt+1
           if(jtt.eq.System_Size(2)+1)jtt=1

           if(js.ne.jr.and.jt.ne.jr)then
             do while (jtt.ne.jr)
                if(yita(i,jt,ka).lt.boundvalue.and.yita(i,jtt,ka).ge.boundvalue.and.counted(i,jtt,ka).eq.0)then
                    SSeedA(nprep,Sstack(nprep,1),2)= jtt; SSeedA(nprep,Sstack(nprep,1),3)= ka
                    Sstack(nprep,1)=Sstack(nprep,1)+1
                endif
                jt=jt+1
                if(jt.eq.System_Size(2)+1)jt=1
                jtt=jt+1
                if(jtt.eq.System_Size(2)+1)jtt=1
              enddo
           endif
        else
           if(js.eq.jr.and.js.eq.jl.and.yita(i,js,ka).ge.boundvalue.and.counted(i,js,ka).eq.0)then
                SSeedA(nprep,Sstack(nprep,1),2)= js; SSeedA(nprep,Sstack(nprep,1),3)= ka
                Sstack(nprep,1)=Sstack(nprep,1)+1
           endif
        endif

        kb=k-1
        if(kb.eq.0) kb=System_Size(3)
        js=jl
        do while(js.ne.jr.and.(yita(i,js,kb).lt.boundvalue.or.counted(i,js,kb).eq.1))
           js=js+1
           if(js.gt.System_Size(2))js=1
        enddo
        if((js.ne.jr.and.yita(i,js,kb).ge.boundvalue).and.counted(i,js,kb).eq.0) then
           SSeedB(nprep,Sstack(nprep,2),2)= js; SSeedB(nprep,Sstack(nprep,2),3)= kb
           Sstack(nprep,2)=Sstack(nprep,2)+1

           jt=js+1
           if(jt.eq.System_Size(2)+1)jt=1
           jtt=jt+1
           if(jtt.eq.System_Size(2)+1)jtt=1
           if(js.ne.jr.and.jt.ne.jr)then
             do while (jtt.ne.jr)
                if(yita(i,jt,kb).lt.boundvalue.and.yita(i,jtt,kb).ge.boundvalue.and.counted(i,jtt,kb).eq.0)then
                    SSeedB(nprep,Sstack(nprep,2),2)= jtt; SSeedB(nprep,Sstack(nprep,2),3)= kb
                    Sstack(nprep,2)=Sstack(nprep,2)+1
                endif
                jt=jt+1
                if(jt.eq.System_Size(2)+1)jt=1
                jtt=jt+1
                if(jtt.eq.System_Size(2)+1)jtt=1
             enddo
           endif
        else
           if(js.eq.jr.and.js.eq.jl.and.yita(i,js,kb).ge.boundvalue.and.counted(i,js,kb).eq.0)then
                SSeedB(nprep,Sstack(nprep,2),2)=js; SSeedB(nprep,Sstack(nprep,2),3)= kb
                Sstack(nprep,2)=Sstack(nprep,2)+1
           endif
        endif

        if(Sstack(nprep,1).gt.1)then  ! There is particle line above
           use_stack(nprep)=1
           j=SSeedA(nprep,Sstack(nprep,1)-1,2); k=SSeedA(nprep,Sstack(nprep,1)-1,3)
           Sstack(nprep,1)=Sstack(nprep,1)-1
        else 
             if(Sstack(nprep,2).gt.1)then  ! There is particle line below
                use_stack(nprep)=1 
                j=SSeedB(nprep,Sstack(nprep,2)-1,2); k=SSeedB(nprep,Sstack(nprep,2)-1,3)
                Sstack(nprep,2)=Sstack(nprep,2)-1
             else
                Is_scanning_P(nprep)=0
                j=landing_P(2); k=landing_P(3)
                nprep=nprep+1
             endif
        endif

     Endif

     if(area(nprep).gt.System_volume/2.or.nprep.ge.MAXPN-2.or.Sstack(nprep,1).gt.40.or.Sstack(nprep,2).gt.40)then
        safety_valve=.False.
        write(12,*)'Safety valve is off, NOT ready for PSD', 'n & area(n) & Stack12:', nprep, &
        area(nprep), Sstack(nprep,1:2)
     endif
     
  Enddo
    
End Subroutine flood_fill_PSD

end module flood_fill