implicit double precision (a-h,o-z)
c
      parameter(npo=4)
      parameter(mlax=51,nlax=mlax+1)
      parameter(mlay=51,nlay=mlay+1)
      dimension p(npo,0:nlax,0:nlay),pp(npo,0:nlax,0:nlay)
      dimension aij(npo,npo),wij(npo,npo),tij(npo,npo),eij(npo,npo)
      dimension mfij(npo,mlax,mlay),nfij(mlax,mlay)
      dimension ls(npo),ms(npo)
c
      dx    = 0.5e-6
      dy    = dx
      gamma = 1.
      delta = 7.*dx
      pree  = 6.2e-6
      activ = 2.08e-19
      tempe = 800
      boltz = 1.38e-23
      amobi = pree*exp(-activ/(boltz*tempe))
      eee   = 1.e+6
      pi    = 3.141592654
c
      aaa   = 2./pi*sqrt(2.*delta*gamma)
	www   = 4.*gamma/delta
      pmobi = amobi*pi*pi/(8.*delta)
c
      dt    = dx*dx/(5.*pmobi*aaa*aaa)
      write(*,*) 'Time increment =',dt,'[s]'
      nend  = 500
c
      do i=1,npo
      do j=1,npo
         wij(i,j)=www
         aij(i,j)=aaa
         tij(i,j)=pmobi
         eij(i,j)=0.
         if(i.eq.j) then
            wij(i,j)=0.
            aij(i,j)=0.
            tij(i,j)=0.
         endif
         if(i.eq.npo.or.j.eq.npo) eij(i,j)=eee
         if(i.gt.j) eij(i,j)=-eij(i,j)
      enddo
      enddo
c
c
      ls(1)=10
      ms(1)=10
      ls(2)=40
      ms(2)=20
      ls(3)=25
      ms(3)=40
c
      r0=3.*dx
      do m=1,mlay
      do l=1,mlax
         phi=0.
         do i=1,npo-1
            r=sqrt((real(l-ls(i))*dx)**2+(real(m-ms(i))*dy)**2)
            r=r-r0
            a=sqrt(2.*www)/aaa*r
            p(i,l,m)=0.5*(1.-sin(a))
            if(a.ge. pi/2.) p(i,l,m)=0.
            if(a.le.-pi/2.) p(i,l,m)=1.
            phi=phi+p(i,l,m)
         enddo
         p(npo,l,m)=1.-phi
      enddo
      enddo
c
      do m=1,mlay
      do l=1,mlax
         do i=1,npo
            pp(i,l,m)=0.
         enddo
      enddo
      enddo
c
      do m=0,nlay
      do i=1,npo
         p(i,0,m)   =p(i,mlax,m)
         p(i,nlax,m)=p(i,1,m)
      enddo
      enddo
      do l=0,nlax
      do i=1,npo
         p(i,l,0)   =p(i,l,mlay)
         p(i,l,nlay)=p(i,l,1)
      enddo
      enddo
c
      nstep=0
      do nstep=1,nend
      write(*,*) 'Step No. =',nstep
c
      do l=1,mlax
      do m=1,mlay
         n=0
         do i=1,npo
            if(p(i,l  ,m  ).gt.0..or.
     &         (p(i,l  ,m  ).eq.0..and.p(i,l+1,m  ).gt.0..or.
     &                                 p(i,l-1,m  ).gt.0..or.
     &                                 p(i,l  ,m+1).gt.0..or.
     &                                 p(i,l  ,m-1).gt.0.)   ) then
            n=n+1
            mfij(n,l,m)=i
            endif
         enddo
         nfij(l,m)=n
      enddo
      enddo
c
      do m=1,mlay
      do l=1,mlax
         do n1=1,nfij(l,m)
            i=mfij(n1,l,m)
            dpi=0.
         do n2=1,nfij(l,m)
            j=mfij(n2,l,m)
            ppp=0.
         do n3=1,nfij(l,m)
            k=mfij(n3,l,m)
            ppp=ppp+
     &        (aij(i,k)*aij(i,k)-aij(j,k)*aij(j,k))/2.
     &        *((p(k,l+1,m)-2.*p(k,l,m)+p(k,l-1,m))/(dx*dx)
     &         +(p(k,l,m+1)-2.*p(k,l,m)+p(k,l,m-1))/(dy*dy))
     &         +(wij(i,k)-wij(j,k))*p(k,l,m)
         enddo
            pee=p(i,l,m)*p(j,l,m)
            dpi=dpi-2.*tij(i,j)/real(nfij(l,m))
     &              *(ppp-8./pi*sqrt(pee)*eij(i,j))
         enddo
            pp(i,l,m)=p(i,l,m)+dpi*dt
         enddo
      enddo
      enddo
c
      do m=1,mlay
      do l=1,mlax
         phi=0.
         do i=1,npo
            if(pp(i,l,m).le.0.) pp(i,l,m)=0.
            if(pp(i,l,m).ge.1.) pp(i,l,m)=1.
            phi=phi+pp(i,l,m)
         enddo
         a=1.
         if(phi.ne.1.) a=1./phi
         do i=1,npo
            p(i,l,m) =a*pp(i,l,m)
            pp(i,l,m)=0.
         enddo
      enddo
      enddo
c
      do m=0,nlay
      do i=1,npo
         p(i,0,m)   =p(i,mlax,m)
         p(i,nlax,m)=p(i,1,m)
      enddo
      enddo
      do l=0,nlax
      do i=1,npo
         p(i,l,0)   =p(i,l,mlay)
         p(i,l,nlay)=p(i,l,1)
      enddo
      enddo
c
      enddo
c
c
      stop
	end