!**************************************************************************************************************************************
!**************************************************************************************************************************************
!***								            	 FRAGMENTED LANDSCAPES PROGRAM			  							            ***
!**************************************************************************************************************************************
!**************************************************************************************************************************************

!	Juliana Rosa Matias Ciccheto and Dr. Sabrina L. B. Araujo
!------------------------------------------------------------------------------------------------------------------------------------------------------------------------
	

! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨ Defining variables
Program MCmodel
	
	Implicit none
 
   integer:: Ni, Nj, B, Ttotal, Repetitions ! Model variabels - Fixed Numbers
   integer:: i,j ,ii, jj, x, y, t, xy, n, nrep, aux, aux2
   integer:: xtry, ytry, xcoord, ycoord, xcoordoffs, ycoordoffs, xcoorddad, ycoorddad
   integer:: Mo, Mt, Dif, m, G, Birth
   integer:: Rside, Neig, Nneig, Nneig1, Nneig2, Nsearch, Nsurv, Dispersion
   integer:: CountPatch, CountSite, Moffspring, FoundSite, FoundFather
   integer:: Nmaxfather, IDInd, IDMother, IDFatherTry, IDFather ! Auxiliary variabels
   integer:: mesp, ncontesp
   
   real:: Tmin, Tmax, K, R, PG, Tmut ! Model variabels - Fixed Numbers
   real:: sij, rand, Ti , h2, Ttry 
   
   integer, allocatable:: CoordDisp(:,:), P(:,:), Psurv(:,:), Poffspring(:,:), Sparents(:,:), Soffspring(:,:), Sfather(:,:,:)
   integer, allocatable:: SeqRandOffs(:), SeqRandDad (:), ListIDDisp(:), Vaux(:), NV(:)

	real, allocatable:: S(:,:), SK(:,:) 
	
	!Allocate(Seed(10))
	!Seed=1
	!call_random_seed(put = Seed)

!------------------------------------------------------------------------------------------------------------------------------------------------------------------------

! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨ 
! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨	Program	¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨
! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨

! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨ Parameters for individuals

	Open(Unit=20,File='inputP.txt',Status='old') 
		read(20,*) K, Ttotal 
		read(20,*)  Tmin, Tmax 
		read(20,*) Tmut 
	Close(20) 

!!!!!!!!!!!!!!!!!!!!!!
	
	OPEN(UNIT=22,FILE='inputRPGB.txt',STATUS='old')
	
	OPEN(UNIT=9,FILE="Temporal.txt",STATUS='new')
	write(9,*) "R ", "PG ", "B ", "Time ", "Individual ", "Latitude ", "Longitude ", "Site_temperature ", "Species ", "Repetitions "
	close(9)
	
	OPEN(UNIT=10,FILE="Species.txt",STATUS='new')
	write(10,*) "R ", "PG ", "B ", "Time ", "ncontesp ", "mesp ", "Repetitions "
	close(10)

! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨ Modelling environmental space (S) 
	
	Open(Unit=21,File='inputS.txt',Status='old') 
		read(21,*) Ni, Nj  
		allocate(S(Ni,Nj), SK(Ni, Nj), Sparents(Ni,Nj), Soffspring(Ni,Nj))
		S=0.0 
		SK=0
	CountPatch=0
	Do i=1, Ni 
		Do j=1, Nj 
		 read(21,*) ii, jj, sij 
		 S(ii,jj)=sij 
		if (sij .ge. Tmin .and. sij .le. Tmax) CountPatch=CountPatch +1
		EndDo	
	EndDo
	Moffspring=NInt(K*CountPatch)
	Mo=NInt(K*Ni*Nj) 
	Close(21)
	
	R=1
	Do while( R .gt. 0) 
		read (22,*) R, PG, B, Repetitions 
		if(R .lt. 1.0) exit 
		G= NInt((1.0-PG)*B) 
		OPEN(UNIT=9,FILE="Temporal.txt",STATUS='old', Access = 'append')
		OPEN(UNIT=10,FILE="Species.txt",STATUS='old', Access = 'append')
	
! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨ Modelling radius (R) and number of neighbors

	!! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨ R
	
	Nneig=0 
	Rside= Int(R)+1 
	Allocate(CoordDisp(4*Rside**2,2)) 
		CoordDisp=0 
	Do i= -Rside, Rside 
		Do j= -Rside, Rside 
			h2= real(i**2 + j**2) !
			if (h2 .le. R**2 .and. h2 .ne. 0.0) then 
				Nneig= Nneig +1 
				CoordDisp(Nneig,1)=i 
				CoordDisp(Nneig,2)=j 
			Endif
		EndDo
	EndDo

	!! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨ R + 1
	
	Nneig1=Nneig
	Do i= -Rside, Rside 
		Do j= -Rside, Rside 
			h2= real(i**2 + j**2) 
			if (h2 .gt. R**2 .and. h2 .le. (R+1)**2) then  
				Nneig= Nneig +1 
				CoordDisp(Nneig,1)=i 
				CoordDisp(Nneig,2)=j 
			Endif
		EndDo
	EndDo
	Nneig2=Nneig
	
! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨ Modelling SK parameter (SK=Capacidade suporte por sítio)

	Do i=1, Ni
		Do j=1, Nj
		CountSite=0
		if(S(i,j) .ge. Tmin .and. S(i,j) .le. Tmax) then
			CountSite=CountSite +1
			Do Neig=1, Nneig1 
				xcoord=i+CoordDisp(Neig,1)
				ycoord=j+CoordDisp(Neig,2)
				if (xcoord .gt. Ni) xcoord=xcoord-Ni 
				if (xcoord .lt. 1) xcoord=Ni+xcoord 
				if (ycoord .gt. Nj) ycoord=ycoord-Nj 
				if (ycoord .lt. 1) ycoord=Nj+ycoord  
				Ttry= S(xcoord,ycoord)
				if (Ttry .ge. Tmin .and. Ttry .le. Tmax) CountSite= CountSite + 1
			EndDo
		Endif
			SK(i,j)=CountSite*K
		EndDo
	EndDo
	
! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨ Modelling individuals 

	Nmaxfather=NInt(K*(Nneig1+1))
	Allocate(SeqRandOffs(Nneig1+1), Sfather(Ni, Nj, Nmaxfather), SeqRandDad(Nneig2+1), NV(Nneig2))
		
	nrep=0	
	Do while(nrep .lt. Repetitions)
	nrep=nrep +1
	
	Allocate(P(Mo, B+3))
	P=0 
	P(:,B+1)=1 
	ncontesp=0
	Do i=1, Mo 
		call random_number (rand) 
		x=1+Int(Ni*rand) 
		call random_number (rand) 
		y=1+Int(Nj*rand) 
		P(i,B+2)=x 
		P(i,B+3)=y 
	EndDo
	
!------------------------------------------------------------------------------------------------------------------------------------------------------------------------

! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨ Programming model dynamics

	Mt=Mo 
	Do t=1, Ttotal 
		Nsurv=0
		Allocate(Psurv(Mt,B+3), ListIDDisp(Mt))
		Psurv=0
		Sfather=0
		Sparents=0
		aux=0
		Do i=1, Mt	
			x=P(i,B+2) 
			y=P(i,B+3) 
			Ti= S(x, y) 
			if (Ti .le. Tmax .and. Ti .ge. Tmin) then 
				Sparents(x, y) = Sparents(x, y) +1 
				Nsurv= Nsurv+1
				Psurv(Nsurv,:)=P(i,:)
				Do n=1, Nmaxfather
					if (Sfather(x, y, n) .eq. 0) then
						Sfather(x, y, n) = Nsurv 
						Exit	
					Endif 
				EndDo
				Do Neig=1, Nneig1
					xcoord= x + CoordDisp(Neig,1)
					ycoord= y + CoordDisp(Neig,2)
					if (xcoord .gt. Ni) xcoord=xcoord-Ni 
					if (xcoord .lt. 1) xcoord=Ni+xcoord 
					if (ycoord .gt. Nj) ycoord=ycoord-Nj 
					if (ycoord .lt. 1) ycoord=Nj+ycoord 					
					Sparents(xcoord, ycoord)= Sparents(xcoord, ycoord) +1
				EndDo
			else
				aux=aux+1
				ListIDDisp(aux)= i
			Endif
		EndDo

	! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨ Dispersion ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨
	
		if(aux .gt. 0) then
			Allocate(Vaux(aux)) 
			call shuffle_1 (Vaux, aux) 
			Do i=1, aux
				IDInd=ListIDDisp(Vaux(i))
				x=P(IDInd,B+2) 
				y=P(IDInd,B+3) 
				Nsearch=0 
				Dispersion=0 
				call shuffle_Nneig(NV, Nneig1, Nneig2)
					Do while(Nsearch .lt. Nneig2 .and. Dispersion .eq. 0) 
						Nsearch=Nsearch+1
						xtry=x+CoordDisp(NV(Nsearch),1) 
						ytry=y+CoordDisp(NV(Nsearch),2) 
						if (xtry .gt. Ni) xtry=xtry-Ni 
						if (xtry .lt. 1) xtry=Ni+xtry 
						if (ytry .gt. Nj) ytry=ytry-Nj 
						if (ytry .lt. 1) ytry=Nj+ytry 
						Ttry= S(xtry,ytry)
						if (Ttry .ge. Tmin .and. Ttry .le. Tmax .and. Sparents(xtry,ytry) .lt. SK(xtry, ytry)) then
							Dispersion=1
							Nsurv= Nsurv+1
							Psurv(Nsurv,:)=P(i,:)
							Psurv(Nsurv, B+2) = xtry
							Psurv(Nsurv, B+3) = ytry
							Do Neig=1, Nneig1
								xcoord= xtry + CoordDisp(Nneig,1)
								ycoord= ytry + CoordDisp(Nneig,2)
								if (xcoord .gt. Ni) xcoord=xcoord-Ni 
								if (xcoord .lt. 1) xcoord=Ni+xcoord 
								if (ycoord .gt. Nj) ycoord=ycoord-Nj 
								if (ycoord .lt. 1) ycoord=Nj+ycoord 	
								Sparents(xcoord, ycoord)= Sparents(xcoord, ycoord) +1
							EndDo
						Endif
					EndDo
			EndDo
			Deallocate(Vaux)
		Endif
		Deallocate(ListIDDisp, P)
		Allocate(P(Nsurv,B+3))
		P=Psurv(1:Nsurv,:)
		Mt=Nsurv
		Deallocate (Psurv)

! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨ Reproduction ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨

		Allocate(Poffspring(Moffspring, B+3))
		Poffspring=0
		Soffspring=0
		Birth=0
		Do i=1, Moffspring 
			call random_number (rand) 
			IDMother=1+Int(Mt*rand)
			call shuffle_0 (SeqRandOffs,Nneig1 + 1, Nneig1 + 1)
			FoundSite=0
			aux=0
			Do while(FoundSite .eq. 0 .and. aux .le. Nneig1) 
				aux=aux +1 
				if(SeqRandOffs(aux) .eq. 0) then
					xcoordoffs=P(IDMother, B+2)
					ycoordoffs=P(IDMother, B+3)
				else
					xcoordoffs=P(IDMother, B+2)+CoordDisp(SeqRandOffs(aux),1)
					ycoordoffs=P(IDMother, B+3)+CoordDisp(SeqRandOffs(aux),2)
					if (xcoordoffs .gt. Ni) xcoordoffs=xcoordoffs-Ni 
					if (xcoordoffs .lt. 1) xcoordoffs=Ni+xcoordoffs 
					if (ycoordoffs .gt. Nj) ycoordoffs=ycoordoffs-Nj 
					if (ycoordoffs .lt. 1) ycoordoffs=Nj+ycoordoffs 
				Endif
				if(Soffspring(xcoordoffs, ycoordoffs) .lt. SK(xcoordoffs, ycoordoffs)) then
					FoundSite=1 

					
! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨ Find a father ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨				

					call shuffle_0 (SeqRandDad, Nneig1 + 1, Nneig2 + 1)
					FoundFather=0
					aux2=0
					Do while(FoundFather .eq. 0 .and. aux2 .le. Nneig2)
						aux2= aux2 +1	
						if(SeqRandDad(aux2) .eq. 0) then
							xcoorddad=P(IDMother, B+2) 
							ycoorddad=P(IDMother, B+3)
						else
							xcoorddad=P(IDMother, B+2)+CoordDisp(SeqRandDad(aux2),1)
							ycoorddad=P(IDMother, B+3)+CoordDisp(SeqRandDad(aux2),2)
							if (xcoorddad .gt. Ni) xcoorddad=xcoorddad-Ni 
							if (xcoorddad .lt. 1) xcoorddad=Ni+xcoorddad 
							if (ycoorddad .gt. Nj) ycoorddad=ycoorddad-Nj 
							if (ycoorddad .lt. 1) ycoorddad=Nj+ycoorddad 	
						Endif
						n=0
						Do while (n .lt. Nmaxfather .and. FoundFather .eq. 0) 
							n=n+1
							IDFatherTry= Sfather(xcoorddad, ycoorddad, n)
							if (IDFatherTry .ne. 0 .and. IDFatherTry .ne. IDMother) then
								Dif=0 
								Do m=1, B
									Dif=Dif+ abs(P(IDMother, m) - P(IDFatherTry, m))
								EndDo
								if ( Dif .le. G) then
									IDFather = IDFatherTry
									FoundFather=1
									Birth=Birth+1
									Do m=1, B
										call random_number (rand) 
										if (rand .lt. 0.5) then 
											Poffspring(Birth, m)=P(IDMother, m)
										else
											Poffspring(Birth, m)=P(IDFather, m)
										Endif
										call random_number (rand) 
										if (rand .le. Tmut) Poffspring(Birth, m) = 1-Poffspring(Birth, m)
									EndDo
									Poffspring(Birth,B+1)=P(IDMother, B+1)
									Poffspring(Birth,B+2)=xcoordoffs
									Poffspring(Birth,B+3)=ycoordoffs
									
! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨ Update Soffspring¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨	

									Soffspring(xcoordoffs, ycoordoffs)= Soffspring(xcoordoffs, ycoordoffs) +1
									Do Neig=1, Nneig1
										xcoord= xcoordoffs + CoordDisp(Neig,1)
										ycoord= ycoordoffs + CoordDisp(Neig,2)
										if (xcoord .gt. Ni) xcoord=xcoord-Ni 
										if (xcoord .lt. 1) xcoord=Ni+xcoord 
										if (ycoord .gt. Nj) ycoord=ycoord-Nj 
										if (ycoord .lt. 1) ycoord=Nj+ycoord 	
										Soffspring(xcoord, ycoord)= Soffspring(xcoord, ycoord) +1
									EndDo
								Endif
							Endif 							
						EndDo
					EndDo
				Endif
			EndDo	
		EndDo
		
		if(Birth .eq. 0) exit
		Mt=Birth
		Deallocate (P)
		Allocate(P(Birth, B+3))
		P=Poffspring(1:Birth,:)
		Deallocate(Poffspring)
		if(t .eq. 100 .or. t .eq. 500 .or. t .eq. 1000 .or. t .eq. 2000 .or.  &
		t .eq. 1000) then
			call species(P, Mt, B+3, G, ncontesp, mesp) 
			Do i=1, Mt
				x= P(i, B+2)
				y= P(i, B+3)
				Ti= S(x, y) 
				write(9,100) R, PG, B, t, i, x, y, Ti, P(i, B+1), nrep 
			EndDo
			write(10,101) R, PG, B, t, ncontesp, mesp, nrep
		Endif
		
	EndDo
	
	Deallocate(P)
	EndDo ! End Do Repetitions
	
	close(9)
	close(10)
	
	Deallocate(CoordDisp, SeqRandDad, SeqRandOffs, NV, Sfather)
	EndDo
	close(22)
	
	100 format(f4.1,1X,f5.3,1X,5(I5,1X), f5.2,1X, I5,1X,I5)
	101 format(f4.1,1X,f5.3,1X,4(I5,1X),1X,I5)
End program MCmodel


!------------------------------------------------------------------------------------------------------------------------------------------------------------------------

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨ Subroutines 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨ Dispersion ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨

subroutine modelling_dispersion(Nsearch, Nneig, Dispersion, x, y, xy, xtry, ytry, Ni, Nj, Ttry, Tmin, Tmax, CoordDisp, S, Rside)
	
	Implicit none
	
	real::Tmax, Tmin, Ttry, rand, S
	integer:: Nsearch, Nneig, Dispersion, x, y, xy, xtry, ytry, Ni, Nj, Rside, CoordDisp
	dimension::S(Ni,Nj), CoordDisp(4*Rside**2,2)

	Do while(Nsearch .le. Nneig .and. Dispersion .eq. 0) 
		Nsearch=Nsearch+1
		call random_number (rand) 
		xy= 1+Int(Nneig*rand) 
		xtry=x+CoordDisp(xy,1) 
		ytry=y+CoordDisp(xy,2) 
		if (xtry .gt. Ni) xtry=xtry-Ni 
		if (xtry .lt. 1) xtry=Ni+xtry
		if (ytry .gt. Nj) ytry=ytry-Nj  
		if (ytry .lt. 1) ytry=Nj+ytry 
		Ttry= S(xtry,ytry)
		if (Ttry .ge. Tmin .and. Ttry .le. Tmax) Dispersion=1
	EndDo
	
endsubroutine modelling_dispersion

!------------------------------------------------------------------------------------------------------------------------------------------------------------------------

! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨ Shuffle: without 0 value ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨

subroutine shuffle_1 (NV,Ni)

	Implicit real*8 (A-H,O-Z)
	dimension::NV(Ni)
  
	Do n=1,Ni
		NV(n)=n
	EndDo
	Do n=1, Ni
		call random_number(h1)
		call random_number(h2)
		i1=1+int(h1*Ni)
		i2=1+int(h2*Ni)
		vh1=NV(i1)
		NV(i1)=NV(i2)
		NV(i2)=vh1
	EndDo	  	
	
endsubroutine shuffle_1

! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨ Shuffle: with 0 value ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨

subroutine shuffle_0(NV,Nneig1, Nneig2)

	Implicit real*8 (A-H,O-Z)
	dimension::NV(Nneig2)
  
	Do n=1,Nneig2
		NV(n)=n-1
	EndDo
	Do n=1, Nneig1
		call random_number(h1)
		call random_number(h2)
		i1=1+int(h1*Nneig1)
		i2=1+int(h2*Nneig1)
		vh1=NV(i1)
		NV(i1)=NV(i2)
		NV(i2)=vh1
	EndDo
	Do n=(Nneig1+1), Nneig2
		call random_number(h1)
		call random_number(h2)
		i1=Nneig1+1+int(h1*(Nneig2-Nneig1))
		i2=Nneig1+1+int(h2*(Nneig2-Nneig1))
		vh1=NV(i1)
		NV(i1)=NV(i2)
		NV(i2)=vh1
	EndDo	  	
	
endsubroutine shuffle_0

! ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨ Shuffle: radius for dispersion and mating ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨

subroutine shuffle_Nneig(NV,Nneig1, Nneig2)

	Implicit real*8 (A-H,O-Z)
	dimension::NV(Nneig2)
  
	Do n=1,Nneig2
		NV(n)=n
	EndDo
	Do n=1, Nneig1
		call random_number(h1)
		call random_number(h2)
		i1=1+int(h1*Nneig1)
		i2=1+int(h2*Nneig1)
		vh1=NV(i1)
		NV(i1)=NV(i2)
		NV(i2)=vh1
	EndDo
	Do n=(Nneig1+1), Nneig2
		call random_number(h1)
		call random_number(h2)
		i1=Nneig1+1+int(h1*(Nneig2-Nneig1))
		i2=Nneig1+1+int(h2*(Nneig2-Nneig1))
		vh1=NV(i1)
		NV(i1)=NV(i2)
		NV(i2)=vh1
	EndDo	  	
	
endsubroutine shuffle_Nneig

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Speciation

	subroutine species(MaV,np,nl,nG,ncontesp,mesp) !(P,num,nl,intq,ncontesp,mesp)
	IMPLICIT REAL*8 (A-H,O-Z)
	Dimension::MaV(np,nl)
	Integer, ALLOCATABLE :: NES(:),Nes2(:),Nfreq(:,:),Nlo(:),Nv(:)
	Allocate(Nlo(np), Nv(np), NES(mesp))
	nfim=0
	nGes=NG
	ngrupo=0
	Nv=0
	Nlo=0
	Do while(nfim.eq.0)
		nfim=1
		Do l=1,np 
			if(Nv(l).eq.0)then 
				nfim=0
				ngrupo=ngrupo+1
				Nv(l)=ngrupo 
				go to 22	
			endif
		Enddo
22  	continue
		nbus=0
		Do while (nbus.eq.0) 
			nbus=1
			Do i=1,np
			if(Nv(i).eq.ngrupo.and.Nlo(i).eq.0) then	
					Nlo(i)=1 
					Do ii=1,np
						if(Nv(ii).eq.0) then 
							ndiff=0
							Do j=1,nl-3 
							 n1=(MaV(i,j))
							 n2=(MaV(ii,j))
							 ndiff=ndiff+abs(n1-n2)
							Enddo
							if(ndiff.le.NG)then 
								nbus=0
								Nv(ii)=ngrupo
							else
								if(ndiff.le.NGes) then
									if(MaV(i,nl-2).eq.MaV(ii,nl-2)) then !!!! antes = nl-1
										nbus=0
										Nv(ii)=ngrupo
									endif
								endif
							endif
						endif
					enddo
				endif
			Enddo
		Enddo
	Enddo
	

	ALLOCATE(Nfreq(ngrupo,2))
	Nes=0
	nmax=mesp
	ncontesp=0 
	Do n=1,nmax
		Nfreq=0
		Do i=1,np
			mva=(MaV(i,nl-2))
			if(mva.eq.n)then 
				Nfreq(nv(i),2)=Nfreq(nv(i),2)+1	
				Nfreq(nv(i),1)=nv(i)	
			endif
		Enddo
		call sortc(Nfreq,ngrupo)
		ngesco=Nfreq(1,1) 
		if(Nes(n).eq.0)then
			Do i=1,np
				IF(Nv(i).eq.ngesco)then
					MaV(i,nl-2)=n 
					Nes(n)=1 
				endif
			Enddo
		endif
		if(Nes(n).gt.0)	ncontesp=ncontesp+1 
		if(ngrupo.gt.1) then
			do nn=2,ngrupo
					if(Nfreq(nn,2).ge.1)then 
						mesp=mesp+1
						ncontesp=ncontesp+1 
						Do i=1,np
							if(Nv(i).eq.Nfreq(nn,1))then
								MaV(i,nl-2)=mesp 
							endif
						enddo
					Endif
			Enddo
		endif
	enddo
endsubroutine species


subroutine sortc(a,n)
    implicit none
    INTEGER :: n, i, j, x, x2, a
    Dimension :: a(n,2)
    DO i = 2, n
        x = a(i,2)
		x2 = a(i,1)
        j = i - 1
        DO WHILE (j.ge.1)
            if (a(j,2).ge.x) exit
            a(j + 1, 2) = a(j,2)
			a(j + 1, 1) = a(j,1)
			j = j - 1
        end do
        a(j + 1, 2) = x
		a(j + 1, 1) = x2
	ENDDO	
	end subroutine sortc
