fortranでDLA

fortranでDLAを書きました。


![test.dat.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3360283/10f60290-a954-d837-d8ff-2d8303037a65.png)

```fortran
  implicit none
  integer,parameter::r=1000
  integer,parameter::n=100000
  real*8::c
  character*4::num
  integer::seedsize,k,x,y,a
  integer,allocatable::seed(:)
  integer::p(r,r),mark(r,r),t,i,j,ii,jj
  integer::kmax,check,count,anum

  
  call random_seed(size=seedsize) 
  allocate(seed(seedsize))        
  call random_seed(get=seed)

  p=0
  p(r/2,r/2)=2
  t=0
  timeloop:do
     t=t+1
     write(0,*)t
     do x=1,r
        do y=1,r
           if(x==1 .or. x==r .or. y==1 .or. y==r)then
              if(p(x,y)==2) then
                 write(0,*)"ok"
                 exit timeloop
              endif
           endif
        enddo
     enddo
     mark=0
     check=0
     do x=1,r
        do y=1,r
           if(p(x,y)==1)check=check+1
        enddo
     enddo
     kmax=n-check
     if(kmax==0)kmax=n+1

     check=0
     aloop:do k=1,kmax
        if(kmax==n+1) exit aloop
        do
           do
              call random_number(c)
              if(c==0d0)cycle
              a=nint(c*real(r))
              if(a==0) cycle
              exit
           enddo
           call random_number(c)
           if(0d0<c .and. c<=0.25d0) then
              x=1
              y=a
           endif
           if(0.25d0<c .and. c<=0.50d0) then
              x=r
              y=a
           endif
           if(0.50d0<c .and. c<=0.75d0) then
              x=a
              y=1
           endif
           if(0.75d0<c .and. c<=1d0)then
              x=a
              y=r
           endif
           if(p(x,y)/=0)then
              if(check==r*r) exit aloop
              check=check+1
              cycle
           endif
           p(x,y)=1
           mark(x,y)=1

           loops1:do ii=-1,1
              do jj=-1,1
                 if(x+ii==0 .or. x+ii==r+1 .or. y+jj==0 .or. y+jj==r+1)cycle
                 if(p(x+ii,y+jj)==2) p(x,y)=2
                 exit loops1
              enddo
           enddo loops1
           exit
        enddo
     enddo aloop

     check=0
     do y=1,r
        do x=1,r
           if(p(x,y)==1 .and. mark(x,y)==0)check=check+1
        enddo
     enddo
     anum=check
     count=0
     do
        do
           call random_number(c)
           x=nint(c*real(r))
           if(x==0) cycle
           exit
        enddo
        do
           call random_number(c)
           y=nint(c*real(r))
           if(y==0) cycle
           exit
        enddo
        if(p(x,y)==1 .and. mark(x,y)==0) then
           count=count+1
              check=0
              do ii=-1,1
                 loopj:do jj=-1,1
                    if(x+ii==0 .or. x+ii==r+1 .or. y+jj==0 .or. y+jj==r+1)then
                       check=check+1
                       cycle
                    endif
                    if(p(x+ii,y+jj)/=0)check=check+1
                 enddo loopj
              enddo
              if(check==9)then
           !      write(0,*)"full"
                 i=0
                 j=0
              else
                 do
                    call random_number(c)
                    if(c==0d0)cycle
                    if(0d0<c .and. c<=0.125d0              )then
                       i=-1
                       j=-1
                    endif
                    if(0.125d0<c .and. c<=0.250d0)then
                       i=0
                       j=-1
                    endif
                    if(0.250d0<c .and. c<=0.375d0)then
                       i=1
                       j=-1
                    endif
                    if(0.375d0<c .and. c<=0.500d0)then
                       i=-1
                       j=0
                    endif
                    if(0.500d0<c .and. c<=0.625d0)then
                       i=1
                       j=0
                    endif
                    if(0.625d0<c .and. c<=0.750d0)then
                       i=-1
                       j=1
                    endif
                    if(0.750d0<c .and. c<=0.875d0)then
                       i=0
                       j=1
                    endif
                    if(0.875d0<c .and. c<=1d0               )then
                       i=1
                       j=1
                    endif
                    if(x+i==0 .or. x+i==r+1 .or. y+j==0 .or. y+j==r+1)cycle
                    if(p(x+i,y+j)/=0) cycle
                    mark(x+i,y+j)=1
                    p(x,y)=0
                    p(x+i,y+j)=1
                    exit
                 enddo
              endif
              loopm:do ii=-1,1
                 do jj=-1,1
                    if(x+i+ii==0 .or. x+i+ii==r+1 .or. y+j+jj==0 .or. y+j+jj==r+1)cycle
                    if(p(x+i+ii,y+j+jj)==2) then
                       p(x+i,y+j)=2
                       exit loopm
                    endif
                 enddo
              enddo loopm
           end if
           if(count==anum)exit
        enddo
     

    if(t/100==real(t)/real(100))then
       !       write(num,'(i4.4)') t
       open(11,file="hoge.dat",status='replace')
       do x=1,r
          do y=1,r
             write(11,*)x,y,p(x,y)
           enddo
        enddo
        close(11)
     endif
  end do timeloop
  
  do x=1,r
     do y=1,r
        print*,x,y,p(x,y)
     enddo
  enddo
end program main
```