program sim_ind_micro
!Author: Toni Pujol
!-----------------------------------------------------------------------
!****f* sim_ind_micro
! NAME 
!    sim_ind_micro
! DESCRIPTION
!
! Main program: this is the main program.
!				It opens the input data file and calls the two main subroutines:
!				lpn_viewfa => calculates the sea-travel trajectory
!				lpn_solite => solves the model equations
!
! VERY IMPORTANT: IF YOU HAVE OPENED comparison.f90, PLEASE CLOSE THAT PROJECT BEFORE COMPILING AND EXECUTING THIS PROGRAM (FILE-->CLOSE WORKSPLACE IF YOU USE FORTRAN POWERSTATION 4.0)
!-----------------------------------------------------------------------
! IMPORTANT: for description of input file parameters: see below
!-----------------------------------------------------------------------
! IMPORTANT: prior to compile this file you should compile the 
!			 def_sim_lpn.f90 file
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
! IMPORTANT: def_sim_lpn.f90
!			 sim_ind_micro.f90
!								have been included into a single project
!								in Microsoft Fortran program and
!								run in Windows XP and Vista environments
!-----------------------------------------------------------------------
!				
! USES
!    lpn_viewfa
!	 lpn_solite
!
! USED BY
!
!-----------------------------------------------------------------------

use def_sim_lpn

implicit none 

integer(4) :: ielem2,npoin,mpoin,jpoin,iy,iboundfront,ictnpepvc
real(8) :: var,xlatmax,xlonmax

! input file

open(12,file='sim_ind_micro.inp')	  

! iboundfront => integer for determining the influence of spatial boundaries: 
!	if: iboundfront = 0 no spatial boundaries; if iboundfront = 1 then we accept mountain boundaries for lat > xlatmax and lon > xlonmax
!
! ictnpepvc => integer for determining the variation of the persistence pe as a function of latitude
!	if ictnpepvc = 0 then pe is constant and equal to cntpep; if ictnpepvc = 1 then pe varies in latitude

read(12,*) iboundfront,ictnpepvc 

! in case iboundfront = 1, mountain boundaris for latitutdes > xlatmax and longitudes > xlonmax

read(12,*) xlatmax,xlonmax

! ifile => integer for determining the input grid: (= 0 homogeneous; = 1 grid with seas and mountains above 1750 m, = 2 grid with seas and mountains above 1000 m)
! icondini => integer for determining the initial conditions: (= 0 Jericó as a single origin; = 1 Hemar as a single origin; = 2 initial condition at 16 PPNB/C)
! itask => integer value for sea travel type (= 1 by means of straight lines, = 0 through the coast line) 

read(12,*) ifile,icondini,itask

! distmax = maximum sea travel distance allowed (x 50 km); ctnpep = persistence 

read(12,*) distmax,ctnpep		

! net reproduction rate (farmers => rtermn; HG => rtermp)
 
read(12,*) rtermn,rtermp

! gamma parameter for vertical transmission
! coeff f coefficient horizontal transmission
! coefpnvisit => term in equation: pn'' = pn' + coefpnvisit pn visit + coefpp pp'
! coefpp => term in equation: pn'' = pn' + coefpnvisit pn visit + coefpp pp'
! coefalfa => frequency term

read(12,*) gamma,coeff,coefpnvisit,coefpp,coefalfa

! if imetodehoritzontalweight = 1 with weights ; if imetodehoritzontalweight = 0 without weights
! if icoeffp = 1, the last term coefpp pp included in equation coefpp, if icoeffp = 0, last term not included
! igenmax = number maximum of generations (= iterations)
														
read(12,*) imetodehoritzontalweight,icoeffp,igenmax

! distranshor = horizontal transmission distance (only for HG) (x 50 km)
! distranshormax = chosen as the maximum horizontal distance available (x 50 km)

read(12,*) distranshor,distranshormax

! outfilen = output file farmers	FORMAT OUTPUT FILE=> x,y,gen (x = x-grid position, y = y-grid position, gen = generations to reach that grid position as indicated in the lpn_output subroutine)

read(12,*) outfilen,outfilep

close (12)

! end input file

!stop


nneig = 8		! number of neighbors in a square lattice


! open grid file	
!
! The grid file format is:
!
! n								= integer that indicates the number of elements in the grid (nelem) (see the 1st read(11,*) command below)
! i   j							= where i is the i-th element (from 1 to n) and j  is the element type (0 = sea, 1 = land, 2 = mountain, 3 = coast line) (see the 2nd read(11,*) command below)
! i	  i1 i2 i3 i4 i5 i6 i7 i8	= where i1 to i8 are the eight neighbors of element i (see the 3rd read(11,*) below)
! i	  xi yi						= where xi and yi are grid coordinates of element i (see the 4th read(11,*) below)
!

if(ifile.eq.0) then

 open(11,file='raig_tot_uns.dat')					! homogeneous grid (no seas, no mountains)

else if(ifile.eq.1) then

 open(11,file='raig_Hemar_mun_revisat.dat')			! grid with seas and mountains above 1750 m

else if(ifile.eq.2) then

 open(11,file='raig_Hemar_mun_1000_revisat.dat')	! grid with seas and mountains above 1000 m

else
	write(*,*) "No input data file"
	stop

end if
		
		read(11,*) nelem							! number of elements in the grid
		
		do ielem = 1, nelem

			read(11,*) lnele(ielem),ltype(ielem)	! lnele = number that identifies the element; ltype = element type: 0 = sea, 1 = land, 2 = mountain, 3 = coast line

		end do

		do ielem = 1, nelem

			read(11,*) ix,(lneig(ineig,ielem),ineig = 1, nneig)		! lneig = neighbors of element ielem
			
		end do 

		do ielem = 1, nelem

			read(11,*) ix, coord(1,ielem), coord(2,ielem)			! grid coordinates of element ielem

		end do


close(11)


! close grid file


distmax = distmax + 0.1			! maximum sea travel distance allowed (x 50 km), add 0.1 for including points that are exactly at distance distmax as defined in the input file

iboun = 0

do ielem = 1, nelem

	if(ltype(ielem).eq.3) then		

		iboun = iboun + 1

		lboun(iboun) = ielem	! coast line points

	end if

end do

nboun = iboun

! Initialize

do ineig = 1, 4

	do ielem = 1, nelem

		lcount(ineig,ielem) = 0

		do ielem2 = 1, 100

			lelecont(ineig,ielem,ielem2) = 0

		end do

	end do

end do

! Initialize

do ielem = 1, nelem

	popdens(ielem) = 0.

	varpp1(ielem) = 0.
	varpp2(ielem) = 0.

	varpn1(ielem) = 0.
	varpn2(ielem) = 0.

end do

100 format(2i15)

101 format(9i15)

102 format(i15,2f15.7)

!! Inicialize neolithic values

densnmx = 1.28*50.*50.  
denspmx	= 0.064*50.*50. 

open(11,file='XYZ_Hemar_munt_revisat_ci.dat')		! initial conditions at 16 PPNB/C points

read(11,*) npoin,mpoin

if (icondini .eq. 0) then

	ielem = 0
	do jpoin = 1, npoin
		do ipoin = 1, mpoin
			ielem = ielem + 1
			if(ipoin.eq.112 .AND. jpoin.eq.23) varpn1(ielem)=densnmx   ! initial point at Jericó (31.8667 Lat, 35.45 Lon) for an Hemar based grid
		end do
	end do

else if(icondini .eq. 1) then

	ielem = 0
	do jpoin = 1, npoin
		do ipoin = 1, mpoin
			ielem = ielem + 1
			if(ipoin.eq.111 .AND. jpoin.eq.22) varpn1(ielem)=densnmx   ! initial point at Hemar (31.17 Lat, 35.18 Lon) for an Hemar based grid
		end do
	end do

else if(icondini .eq. 2) then

	ielem = 0

	do jpoin = 1, npoin

		do ipoin = 1, mpoin

			ielem = ielem + 1

			read(11,*) ix,iy,var
		
			varpn1(ielem)=var*densnmx
		
		end do

	end do


else

	write(*,*) 'No initial conditions data'
	stop

end if

close(11)

!! end initialize neolithic values

!! Initialize mesolithic values

do ielem = 1, nelem

	if(ltype(ielem).eq.1 .OR. ltype(ielem).eq. 3) then

		if(varpn1(ielem).gt.0.) then
		
			varpp1(ielem) = 0.
			
		else
		
			varpp1(ielem) = denspmx
			
		end if
		
	end if
	
end do

!! end initialize mesolitchi values

! Albers projection
!

! 1 and 2 reference latitudes

varpi = 3.141592654

reflat1 = 30.*varpi/180.
reflat2 = 65.*varpi/180.

! Hemar latitude and longitude (grid origin)

vorilat = 31.17*varpi/180.
vorilon = 35.18*varpi/180.

! Coefficients of Albers projections

coefn = .5*(sin(reflat1) + sin(reflat2))
coefc = cos(reflat1)*cos(reflat1)+2.*coefn*sin(reflat1)
coefr = sqrt(coefc-2.*coefn*sin(vorilat))/coefn

ielem = 0

do jpoin = 1, npoin

	termy = (real(jpoin) - 1.) * .00784806 -	0.16480926

	do ipoin = 1, mpoin

		ielem = ielem + 1

		termx = (real(ipoin) - 1.) * .00784806 - 0.8632866

		theta = atan(termx/(coefr - termy))

		xlon = vorilon + theta/coefn
		
		xlat = asin( (coefc - (termx / sin(theta) * coefn )**2.)/ (2.* coefn) )

		if(ictnpepvc.eq.1) then

			ctnpepv(ielem) = 0.44 * (xlat * 180. / varpi - 6.) / (67.- 6.)		! variation of pe as a function of latitude (change if necessary)
			ctnpenv(ielem) = 0.44 * (xlat * 180. / varpi - 6.) / (67.- 6.) 		! variation of pe as a function of latitude (change if necessary)

		else

! ctnpepv;ctnpenv:  minimum = 0.64, maximum = 0.12, mean = 0.38

			ctnpepv(ielem) = ctnpep
			ctnpenv(ielem) = ctnpep

		end if

! Spatial boundaries

		if(iboundfront.eq.1) then

			if(xlat*180./varpi .gt. xlatmax .AND. xlon*180./varpi.gt.xlonmax) then
			
				if(ltype(ielem).eq.1) ltype(ielem)=2								! all land => mountains

			end if

		end if

	end do

end do


!!! initialize variables

do ielem = 1, nelem

	weight(ielem) = 0.

	varpp2(ielem) = 0.

end do


!! End initial conditions


call lpn_viewfa

call lpn_solite

end 

!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lpn_viewfa
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-----------------------------------------------------------------------
!****f* lpn_viewfa
! NAME 
!    lpn_viewfa
!
! DESCRIPTION
!	 Search for the grid points available for sea travel    
!
! USES
!    lpn_elecon
!
! USED BY
!	sim_ind_micro
!	 
!***
!-----------------------------------------------------------------------

use def_sim_lpn

implicit none 

integer*4 :: ielemo

!!

do ielem = 1, nelem

	write(*,*) 'viewfa ', ielem

	ielemo = ielem

	call lpn_elecon(ielemo)
	
end do	


end subroutine lpn_viewfa


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lpn_elecon(ielemo)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-----------------------------------------------------------------------
!****f* lpn_elecon
! NAME 
!    lpn_elecon
!
! DESCRIPTION
!	 Obtains the points that can contribute to ielemo (= origin element) 
!
!		    x
!			!
!			v	
!	 x -->	O <-- x		   O = origin element (ielemo)
!			^			   x = ielemneighbor = lelecont[ineigo,ielemo,lcont(ineigo,ielemo)]
!			!
!			x			   
!						   
!
!	  if ltype(ielemo) = 1 (land) 
!			if ltype[lneig(ineigo,ielemo)] = 1, 3 (land or coast line) => lcount(ineigo,ielemo) = 1 
!			if ltype[lneig(ineigo,ielemo)] = 2 (mountain)  => lcount(ineigo,ielemo) = 0
!			if ltype[lneig(ineigo,ielemo)] = 0 (sea)  => Error (between sea and land there shoud be coast line (ltype = 3)
!
!	  if ltype(ielemo) = 3 (coast line) => if lcount(ineigo,ielemo) = 0 (sea) => call lpn_mardis
!			if ltype[lneig(ineigo,ielemo)] = 1, 3 (land or coast line) => lcount(ineigo,ielemo) = 1 
!			if ltype[lneig(ineigo,ielemo)] = 2 (mountain)  => lcount(ineigo,ielemo) = 0
!			if ltype[lneig(ineigo,ielemo)] = 0 call lpn_mardis (=> sea travel!!)
!
!
!
! USES
!     lpn_mardis
!
! USED BY
!	  lpn_viewfa
!
!***
!-----------------------------------------------------------------------

use def_sim_lpn

implicit none 

integer*4 ineigo,ielemo


if(ltype(ielemo).eq.0 .OR. ltype(ielemo).eq.2) then  !! sea or mountain

	popdens(ielemo) = 0.			! weight for population


else if(ltype(ielemo).eq.1) then		!! ground

	popdens(ielemo) = 1.

	do ineigo = 1, 4
		
		if(lneig(ineigo,ielemo).gt.0) then

			if(ltype(lneig(ineigo,ielemo)).eq.0) then

				write(*,*) 'go from 1 to 0 without 3 in the middle ',ielemo, ineigo,lneig(ineigo,ielemo)
				
				ltype(lneig(ineigo,ielemo)) = 3
				
				stop 

			else if(ltype(lneig(ineigo,ielemo)).eq.1) then

				lcount(ineigo,ielemo) = 1

				lelecont(ineigo,ielemo,lcount(ineigo,ielemo)) = lneig(ineigo,ielemo)

			else if(ltype(lneig(ineigo,ielemo)).eq.2) then

				lcount(ineigo,ielemo) = 0

			else if(ltype(lneig(ineigo,ielemo)).eq.3) then

				lcount(ineigo,ielemo) = 1

				lelecont(ineigo,ielemo,lcount(ineigo,ielemo)) = lneig(ineigo,ielemo)

			end if

		end if

	end do


else if(ltype(ielemo).eq.3) then		!! coast line

	popdens(ielemo) = 1.

	do ineigo = 1, 4

		if(lneig(ineigo,ielemo).gt.0) then

			if(ltype(lneig(ineigo,ielemo)).eq.0) then

				call lpn_mardis(ineigo,ielemo)
							
			else if(ltype(lneig(ineigo,ielemo)).eq.1) then

				lcount(ineigo,ielemo) = 1

				lelecont(ineigo,ielemo,lcount(ineigo,ielemo)) = lneig(ineigo,ielemo)

			else if(ltype(lneig(ineigo,ielemo)).eq.2) then

				lcount(ineigo,ielemo) = 0

			else if(ltype(lneig(ineigo,ielemo)).eq.3) then

				lcount(ineigo,ielemo) = 1

				lelecont(ineigo,ielemo,lcount(ineigo,ielemo)) = lneig(ineigo,ielemo)

			end if

		end if

	end do

end if

end subroutine lpn_elecon

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lpn_mardis(ineigo,ielemo)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-----------------------------------------------------------------------
!****f* lpn_mardis
! NAME 
!    lpn_mardis
! DESCRIPTION
!
!	 if itask = 0 => sea travel through the coast line
!	 if itask = 1 => sea travel through straight lines
!

! USES
!    lpn_marcos
!	 lpn_marcom

! USED BY
!	 lpn_elecon
!***
!-----------------------------------------------------------------------

use def_sim_lpn

implicit none 

integer*4 :: ineigo,ielemo,ivora
	   
ivora = 0

!! in case that the point corresponds to the edges

do ineig = 1, nneig

	if(lneig(ineig,ielemo).eq.0) ivora = 1

end do

if(ivora.eq.0) then

	SELECT CASE (itask)

	CASE(0)			  !! itask = 0 => travel through the coast line

		call lpn_marcos(ineigo,ielemo)

		write(*,*) ' call marcos = ', ielemo

	CASE(1)			  !! itask = 1 => travel trhough the sea (like straight lines)

		call lpn_marcom(ineigo,ielemo)
		
		write(*,*) ' call marcom = ', ielemo

	CASE(2)

!!


	END SELECT

end if

end subroutine lpn_mardis

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lpn_marcos(ineigo,ielemo)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!-----------------------------------------------------------------------
!****f* lpn_marcos
! NAME 
!    lpn_marcos
! DESCRIPTION
!		
!	 sea travel through the coast line
!
!
! USES
!   lpn_eledis

! USED BY
!   lpn_mardis
!***
!-----------------------------------------------------------------------

use def_sim_lpn

implicit none 

integer(4) :: ielemo,ineigo,idummy,icosta,icostamax
integer(4) :: ielem1,ielem2,ineigi3(8),ineig3,ielemneig,itestneig(8)

real*8 :: distelem,dist2elem,distini(8)

!!! Initialize

do ineig = 1, nneig

	itestneig(ineig) = 0.

	ineigi3(ineig) = 0.

	distini(ineig) = 0.

end do

!!! Coast distribution aroung point "ielemo"

icosta = 0

do ineig = 1, 4

	ielemneig = lneig(ineig,ielemo)

	if(ielemneig.gt.0) then

		if(ltype(ielemneig).eq.3) then

			icosta = icosta + 1
			
			ineigi3(icosta) = ineig

			itestneig(ineig) = 1

			if(ineigo.eq.1) then
			
				if(ineig.eq.3) then

					distini(icosta) = 3.
			  
			  	else if(ineig.eq.2 .OR. ineig.eq.4) then
			
					distini(icosta) = sqrt(2.)

			  	else if(ineig.eq.1) then
					
					distini(icosta) = 1.
			
				end if 
			
			else if(ineigo.eq.2) then

				if(ineig.eq.4) then

					distini(icosta) = 3.

			  	else if(ineig.eq.1 .OR. ineig.eq.3) then
			
					distini(icosta) = sqrt(2.)

			  	else if(ineig.eq.2) then

					distini(icosta) = 1.

				end if

			else if(ineigo.eq.3) then

				if(ineig.eq.1) then

					distini(icosta) = 3.

			  	else if(ineig.eq.2 .OR. ineig.eq.4) then
			
					distini(icosta) = sqrt(2.)

			  	else if(ineig.eq.3) then

					distini(icosta) = 1.

				end if

			else if(ineigo.eq.4) then

				if(ineig.eq. 2) then

					distini(icosta) = 3.

			  	else if(ineig.eq.1 .OR. ineig.eq.3) then
			
					distini(icosta) = sqrt(2.)

			  	else if(ineig.eq.4) then

					distini(icosta) = 1.

				end if

			end if

		end if

	end if

end do

do ineig = 5, nneig - 1 

	ielemneig = lneig(ineig,ielemo)

	if(ielemneig.gt.0) then

		if(ltype(ielemneig).eq.3 .AND. (itestneig(ineig - 4).eq.0 .OR. itestneig(ineig - 3).eq.0)) then

			icosta = icosta + 1
			
			ineigi3(icosta) = ineig

			itestneig(ineig) = 1

			if(ineigo.eq.1) then
			
				if(ineig.eq.5) then

					distini(icosta) = sqrt(2.)
			  
			  	else
			
					distini(icosta) = 2. + sqrt(2.)
			
				end if 
			
			else if(ineigo.eq.2) then

				if(ineig.eq.5 .OR. ineig.eq.6) then

					distini(icosta) = sqrt(2.)

				else

					distini(icosta) = 2. + sqrt(2.)

				end if

			else if(ineigo.eq.3) then

				if(ineig.eq.6 .OR. ineig.eq.7) then

					distini(icosta) = sqrt(2.)

				else

					distini(icosta) = 2. + sqrt(2.)

				end if

			else if(ineigo.eq.4) then

				if(ineig.eq. 7) then

					distini(icosta) = sqrt(2.)

				else

					distini(icosta) = 2. + sqrt(2.)

				end if

			end if


		end if

	end if

end do

ielemneig = lneig(nneig,ielemo)

if(ielemneig.gt.0) then

	if(ltype(ielemneig).eq.3 .AND. (itestneig(nneig - 4).eq.0 .OR. itestneig(1).eq.0)) then

		icosta = icosta + 1
			
		ineigi3(icosta) = nneig

		itestneig(nneig) = 1

		if(ineigo.eq.1 .OR. ineigo.eq.4) then

			distini(icosta) = sqrt(2.)

		else

			distini(icosta) = 2. + sqrt(2.)

		end if

	end if

end if

icostamax = icosta


!! Following the coastline


do icosta = 1, icostamax

	ineig3 = ineigi3(icosta)

	ielem1 = ielemo

	ielem2 = lneig(ineig3,ielem1)

	distelem = distini(icosta)

	if(distelem.lt.distmax) then

		lcount(ineigo,ielemo) = 1 + lcount(ineigo,ielemo)

		lelecont(ineigo,ielemo,lcount(ineigo,ielemo)) = ielem2
			
	end if

	idummy = 0

!! Find coast point for neighbor; following the coastline

	do while(distelem.lt.distmax .AND. idummy.lt.nboun)

		idummy = idummy + 1

		ineig3 = 0

		do ineig = 1, nneig
			
			ielemneig = lneig(ineig,ielem2)

			if(ielemneig.gt.0) then
				
					if(ltype(lneig(ineig,ielem2)).eq.3 .AND. lneig(ineig,ielem2).ne.ielem1) ineig3 = ineig

			end if

		end do

		if(ineig3.gt.0) then

			ielem1 = ielem2

			ielem2 = lneig(ineig3,ielem2)

			call lpn_eledis(ielem1,ielem2,dist2elem)

			distelem = distelem + dist2elem


			!! save
			
			if(distelem .lt. distmax) then

				lcount(ineigo,ielemo) = 1 + lcount(ineigo,ielemo)

				lelecont(ineigo,ielemo,lcount(ineigo,ielemo)) = ielem2
			
			end if

		end if

	end do

end do


end subroutine lpn_marcos

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lpn_marcom(ineigo,ielemo)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-----------------------------------------------------------------------
!****f* lpn_marcom
! NAME 
!    lpn_marcom
! DESCRIPTION
!	
!	 sea travel by straight lines
!
! USES
!    lpn_eledis
!	 lpn_marslo
!	 lpn_marcon
!	 lpn_marpat
!	 lpn_marwat
!

! USED BY
!	lpn_elecon
!***
!-----------------------------------------------------------------------

use def_sim_lpn

implicit none 

integer(4) :: ipoinmin1,ipoinmin2,ielem1,ielem2,ielem3,ielem4,ielem5
integer(4) :: ibouno,ibound,iwatpath,ideltax,ideltay,ielemo,ineigo

real*8 :: distelem

!! Find "ibouno" 

do iboun = 1, nboun

	if(ielemo.eq.lboun(iboun)) ibouno = iboun

end do

!! Rays for all boundaries (number 3)


do ibound = 1, ibouno - 1

	iter = 0

	ielem1 = lboun(ibouno)
	ielem2 = lboun(ibound)

	call lpn_eledis(ielem1,ielem2,distelem)
		
	if(distelem.lt.distmax) then

		call lpn_marslo(ibouno,ibound)

		ideltax = int(coord(1,lboun(ibouno))-coord(1,lboun(ibound)))
		ideltay = int(coord(2,lboun(ibouno))-coord(2,lboun(ibound)))

		itmax = max(abs(ideltax),abs(ideltay))

		iwatpath = 1

		do while(iter.lt.itmax .AND. iwatpath.ne.0)

			iter = iter + 1

			call lpn_marpat(ibouno)
			
			if(iter.eq.1) then
			
				ielem1 = lboun(ibouno)
	
				call lpn_marcon(ielem1,ipoinmin1,ipoinmin2)

				ielem2 = lneig(ipoinmin1,ielem1)

				ielem3 = lneig(ipoinmin2,ielem1)
				
				if(itmax.ne.1) then

					if(ltype(ielem2).ne.0 .AND.ltype(ielem3).ne.0) then

						iwatpath = 0

					end if

				end if

			else 
				 
				call lpn_marcon(ielem2,ipoinmin1,ipoinmin2)

				ielem4 = lneig(ipoinmin1,ielem2)
				 
				ielem5 = lneig(ipoinmin2,ielem2)
					
				call lpn_marwat(ielem2,ielem3,ielem4,ielem5,iwatpath)

				if(iter.eq.itmax) iwatpath = 1
				
				ielem2 = ielem4

				ielem3 = ielem5

			end if
		
		end do

		if (iwatpath.ne.0) then

			lcount(ineigo,ielemo) = 1 + lcount(ineigo,ielemo)

			lelecont(ineigo,ielemo,lcount(ineigo,ielemo)) = lboun(ibound)

		end if		

	end if

end do

!!! Repeat for an ibound greater than ibouno

do ibound = ibouno + 1, nboun

	iter = 0

	ielem1 = lboun(ibouno)
	ielem2 = lboun(ibound)

	call lpn_eledis(ielem1,ielem2,distelem)
		
	if(distelem.lt.distmax) then

		call lpn_marslo(ibouno,ibound)

		ideltax = int(coord(1,lboun(ibouno))-coord(1,lboun(ibound)))
		ideltay = int(coord(2,lboun(ibouno))-coord(2,lboun(ibound)))

		itmax = max(abs(ideltax),abs(ideltay))

		iwatpath = 1

		do while(iter.lt.itmax .AND. iwatpath.ne.0)

			iter = iter + 1

			call lpn_marpat(ibouno)
			
			if(iter.eq.1) then
			
				ielem1 = lboun(ibouno)
	
				call lpn_marcon(ielem1,ipoinmin1,ipoinmin2)

				ielem2 = lneig(ipoinmin1,ielem1)

				ielem3 = lneig(ipoinmin2,ielem1)
				
				if(itmax.ne.1) then

					if(ltype(ielem2).ne.0 .AND. ltype(ielem3).ne.0) then

						iwatpath = 0

					end if

				end if

			else 
				 
				call lpn_marcon(ielem2,ipoinmin1,ipoinmin2)

				ielem4 = lneig(ipoinmin1,ielem2)
				 
				ielem5 = lneig(ipoinmin2,ielem2)
					
				call lpn_marwat(ielem2,ielem3,ielem4,ielem5,iwatpath)
				
				if(iter.eq.itmax) iwatpath = 1

				ielem2 = ielem4

				ielem3 = ielem5

			end if
		
		end do

		if (iwatpath.ne.0) then

			lcount(ineigo,ielemo) = 1 + lcount(ineigo,ielemo)

			lelecont(ineigo,ielemo,lcount(ineigo,ielemo)) = lboun(ibound)

		end if		

	end if
	
end do

end subroutine lpn_marcom

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lpn_marslo(ibouno,ibound)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!-----------------------------------------------------------------------
!****f* lpn_marslo
! NAME 
!    lpn_marslo
! DESCRIPTION
!
!	 find the slope and the quadrant of the ray trajectory
!
!
! USES
!

! USED BY
!    lpn_marcom
!***
!-----------------------------------------------------------------------

use def_sim_lpn

implicit none 

integer(4) :: ibouno,ibound

xorigen = coord(1,lboun(ibouno))
yorigen = coord(2,lboun(ibouno))

xdestin = coord(1,lboun(ibound))
ydestin	= coord(2,lboun(ibound))

if(xdestin.eq.xorigen) then

	if(ydestin.gt.yorigen) then

		slope = 0.
		islop = 1	
		iquad = 0

	else

		slope =	 0.
		islop = -1
		iquad =  0

	end if

else    
	
	slope = (ydestin - yorigen) / (xdestin - xorigen)
	islop = 0

	if(xdestin.gt.xorigen .AND. ydestin.ge.yorigen) then

		iquad = 1

	else if(xdestin.lt.xorigen .AND. ydestin.ge.yorigen) then

		iquad = 2

	else if(xdestin.lt.xorigen .AND. ydestin.le.yorigen) then

		iquad = 3

	else if(xdestin.gt.xorigen .AND. ydestin.le.yorigen) then

		iquad = 4

	else

		write(*,*) 'marslo - quadrants'
		stop

	end if

end if


end subroutine lpn_marslo


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lpn_marpat(ibouno)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-----------------------------------------------------------------------
!****f* lpn_marpat
! NAME 
!    lpn_marpat
! DESCRIPTION
!
!	 find the coordinates for the ray trajectory at the corresponding iteration


! USES
!    

! USED BY
!  lpn_marcom
!***
!-----------------------------------------------------------------------

use def_sim_lpn

implicit none 

integer(4) :: ibouno



if(islop.eq.0) then

	if(slope.ge.-1. .AND. slope.lt.1.) then

	   if(iquad.eq.1) then

			coordraig(1) = coord(1,lboun(ibouno)) + iter  
			coordraig(2) = coord(2,lboun(ibouno)) + slope * iter 
	
		else if(iquad.eq.2) then

			coordraig(1) = coord(1,lboun(ibouno)) - iter  
			coordraig(2) = coord(2,lboun(ibouno)) - slope * iter 

		else if(iquad.eq.3) then

			coordraig(1) = coord(1,lboun(ibouno)) - iter  
			coordraig(2) = coord(2,lboun(ibouno)) - slope * iter 

		else if(iquad.eq.4) then

			coordraig(1) = coord(1,lboun(ibouno)) + iter  
			coordraig(2) = coord(2,lboun(ibouno)) + slope * iter 

		end if			
	
	else

		if(iquad.eq.1) then

			coordraig(1) = coord(1,lboun(ibouno)) + iter / slope
			coordraig(2) = coord(2,lboun(ibouno)) + iter 

		else if(iquad.eq.2) then

			coordraig(1) = coord(1,lboun(ibouno)) + iter / slope
			coordraig(2) = coord(2,lboun(ibouno)) + iter 

		else if(iquad.eq.3) then

			coordraig(1) = coord(1,lboun(ibouno)) - iter / slope
			coordraig(2) = coord(2,lboun(ibouno)) - iter 

		else if(iquad.eq.4) then

			coordraig(1) = coord(1,lboun(ibouno)) - iter / slope
			coordraig(2) = coord(2,lboun(ibouno)) - iter 

		end if

	end if

else if(islop.eq.1) then

	coordraig(1) = coord(1,lboun(ibouno))
	coordraig(2) = coord(2,lboun(ibouno)) + iter

else if(islop.eq.-1) then

	coordraig(1) = coord(1,lboun(ibouno))
	coordraig(2) = coord(2,lboun(ibouno)) - iter

end if


end subroutine lpn_marpat


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lpn_marcon(ielemx,ipoinmin1,ipoinmin2)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-----------------------------------------------------------------------
!****f* lpn_marcon
! NAME 
!    lpn_marcon
! DESCRIPTION
!
! Find those points closest to rays
!

! USES
!

! USED BY
!	 lpn_marcom
!***
!-----------------------------------------------------------------------

use def_sim_lpn

implicit none 

integer(4) :: ipoinmin1,ipoinmin2,ielemx

! npoin = number of points arount origin

! two neares points

! xcoord = coord(1,lnele(1,lboun(iboun))
! ycoord = coord(2,lnele(1,lboun(iboun))

xcoord = coord(1,lneig(1,ielemx))
ycoord = coord(2,lneig(1,ielemx))

ipoinmin1 = 1
distmin1 = (xcoord-coordraig(1))*(xcoord-coordraig(1)) + (ycoord-coordraig(2))*(ycoord-coordraig(2))


xcoord = coord(1,lneig(2,ielemx))
ycoord = coord(2,lneig(2,ielemx))

ipoinmin2 = 2
distmin2 = (xcoord-coordraig(1))*(xcoord-coordraig(1)) + (ycoord-coordraig(2))*(ycoord-coordraig(2))


if(distmin2.lt.distmin1) then

	ipoinmind = ipoinmin1
	ipoinmin1 = ipoinmin2
	ipoinmin2 = ipoinmind

	distmind = distmin1
	distmin1 = distmin2
	distmin2 = distmind

end if


do ineig = 3, nneig

	xcoord = coord(1,lneig(ineig,ielemx))
	ycoord = coord(2,lneig(ineig,ielemx))

	distcoord = (xcoord-coordraig(1))*(xcoord-coordraig(1)) + (ycoord-coordraig(2))*(ycoord-coordraig(2))
	

	if(distcoord.lt.distmin1 .AND. distcoord.lt.distmin2) then

		ipoinmin1 = ineig
		distmin1 = distcoord

	else if(distcoord.ge.distmin1 .AND. distcoord.le.distmin2) then

		ipoinmin2 = ineig
		distmin2 = distcoord
		
	end if

end do

!! Nearest points

!! ipoinmin2

!! ipoinmin1

!! For equal distances:
!! 
!!	 6  2  5
!!	 3	   1    => ray through the middle of two neighbors
!!	 7	4  8
!!

if(distmin1.eq.distmin2) then

	if(ipoinmin1.lt.ipoinmin2) then

		ipoinmind = ipoinmin2
		ipoinmin2 = ipoinmin1
		ipoinmin1 = ipoinmind

	end if

end if

!! If distance equal 0

if(distmin1.eq.0) then

	if(ipoinmin1.le.4) then

		if(ltype(lneig(ipoinmin1 + 4,ielemx)).eq.0) then
		
			distmin2 = 1.
			ipoinmin2 = ipoinmin1 + 4
			
		else
		
			distmin2 = 1.

			if(ipoinmin1.eq.1) then

				ipoinmin2 = 8

			else
	
				ipoinmin2 = ipoinmin1 + 3
	
			end if
				
		end if
	
	else 

		if(ltype(lneig(ipoinmin1 - 4,ielemx)).eq.0) then
		
			distmin2 = 1.
			ipoinmin2 = ipoinmin1 - 4
			
		else
		
			distmin2 = 1.

			if(ipoinmin1.eq.8) then

				ipoinmin2 = 1

			else
	
				ipoinmin2 = ipoinmin1 - 3
	
			end if
				
		end if

	end if

end if
		
! end if distance = 0

	
end subroutine lpn_marcon

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lpn_marwat(ielem2,ielem3,ielem4,ielem5,iwatpath)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!-----------------------------------------------------------------------
!****f* codi_pal_neo_ind_micro/lpn_marwat
! NAME 
!    lpn_marwat
! DESCRIPTION
!
!	 Determines whether or not the straight line travels between sea points
!
!
! USES
!    

! USED BY
!	lpn_marcom
!***
!-----------------------------------------------------------------------

use def_sim_lpn

implicit none 

integer(4) :: ielem2,ielem3,ielem4,ielem5,iwatpath

xelem2 = coord(1,ielem2)
yelem2 = coord(2,ielem2)

xelem3 = coord(1,ielem3)
yelem3 = coord(2,ielem3)

xelem4 = coord(1,ielem4)
yelem4 = coord(2,ielem4)

xelem5 = coord(1,ielem5)
yelem5 = coord(2,ielem5)

dist24 = (xelem2 - xelem4)*(xelem2 - xelem4) + (yelem2 - yelem4)*(yelem2 - yelem4)
dist25 = (xelem2 - xelem5)*(xelem2 - xelem5) + (yelem2 - yelem5)*(yelem2 - yelem5)

dist34 = (xelem3 - xelem4)*(xelem3 - xelem4) + (yelem3 - yelem4)*(yelem3 - yelem4)
dist35 = (xelem3 - xelem5)*(xelem3 - xelem5) + (yelem3 - yelem5)*(yelem3 - yelem5)


if(ltype(ielem2).eq.0) then

	if(dist24.le.1.01) then

		if(ltype(ielem4).eq.0) then
		
			iwatpath = 1

			return
			
		end if 
		
	 else if(dist25.le.1.01) then
	 
		if(ltype(ielem5).eq.0) then
		
			iwatpath = 1
			
			return
			
		end if
	
	 end if

end if

if(ltype(ielem3).eq.0) then

	if(dist34.le.1.01) then

		if(ltype(ielem4).eq.0) then
		
			iwatpath = 1

			return
			
		end if 
		
	 else if(dist35.le.1.01) then
	 
		if(ltype(ielem5).eq.0) then
		
			iwatpath = 1
			
			return
			
		end if
	
	 end if
	
end if

iwatpath = 0
	

end subroutine lpn_marwat

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lpn_eledis(ielem1,ielem2,distelem)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!-----------------------------------------------------------------------
!****f* lpn_eledis
! NAME 
!    lpn_eledis
! DESCRIPTION
!
!	  Calculate distance between ielem1 and ielem2
!
! USES
!   
!
! USED BY
!	lpn_marcos
!   lpn_marcom
!***
!-----------------------------------------------------------------------

use def_sim_lpn

implicit none 

integer(4) :: ielem1,ielem2

real*8 :: distelem

xelem1 = coord(1,ielem1)
yelem1 = coord(2,ielem1)

xelem2 = coord(1,ielem2)
yelem2 = coord(2,ielem2)

distelem = sqrt((xelem1 - xelem2)*(xelem1 - xelem2) + (yelem1 - yelem2)*(yelem1 - yelem2))

end subroutine lpn_eledis


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lpn_solite
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!-----------------------------------------------------------------------
!****f* codi_pal_neo_ind_micro/lpn_solite
! NAME 
!    lpn_solite
!
! DESCRIPTION
!	 starts the solution routine
!
! USES
!	lpn_popite
!   lpn_output
!    
!
! USED BY
!	  sim_in_micro
!
!***
!--------------------------------------------------------------------------
!

use      def_sim_lpn

!

implicit none

integer*4 :: jelem

! weigth calculus 

do ielem = 1, nelem

	do ineig = 1, 4

		if(lcount(ineig,ielem).gt.0.) then

			weight(ielem) = weight(ielem) + 1.

		end if

	end do

end do


distranshor0 = 2.01		! horizontal transmission distance (x 50 km = 100 km)
distranshor1 = 4.01		! horizontal transmission distance (x 50 km = 200 km)
distranshor2 = 8.01		! horizontal transmission distance (x 50 km = 400 km)
distranshor3 = 12.01	! horizontal transmission distance (x 50 km = 600 km)
distranshor4 = 16.01	! horizontal transmission distance (x 50 km = 800 km)

do ielem = 1, nelem

	xohor = coord(1,ielem)
	yohor = coord(2,ielem)

	ihor0 = 0
	ihor1 = 0
	ihor2 = 0
	ihor3 = 0
	ihor4 = 0

	do jelem = 1, nelem

		xdhor = coord(1,jelem)
		ydhor = coord(2,jelem)

		disthor = sqrt((xdhor-xohor)*(xdhor-xohor) + (ydhor-yohor)*(ydhor-yohor))

		if(imetodehoritzontalweight.eq.1) then

			if(disthor .le. distranshor0 .AND. jelem.ne.ielem) then		  ! < 100 km

				ihor0 = ihor0 + 1
				iodhor0(ihor0,ielem) = jelem

			end if

			if(disthor .le. distranshor1 .AND. disthor .gt. distranshor0) then		! 200 > and > 100 km

				ihor1 = ihor1 + 1
				iodhor1(ihor1,ielem) = jelem

			end if

			if(disthor .le. distranshor2 .AND. disthor .gt. distranshor1) then		! 400 > and > 200 km

				ihor2 = ihor2 + 1
				iodhor2(ihor2,ielem) = jelem

			end if

			if(disthor .le. distranshor3 .AND. disthor .gt. distranshor2) then		 ! 600 > and > 400 km

				ihor3 = ihor3 + 1
				iodhor3(ihor3,ielem) = jelem

			end if

			if(disthor .le. distranshor4 .AND. disthor .gt. distranshor3) then		 ! 800 > and > 600 km

				ihor4 = ihor4 + 1
				iodhor4(ihor4,ielem) = jelem

			end if
	
		else

			if(disthor .le. distranshor .AND. jelem.ne.ielem) then

				ihor4 = ihor4 + 1
				iodhor4(ihor4,ielem) = jelem
			end if

		end if

	end do

	nhor0(ielem) = ihor0
	nhor1(ielem) = ihor1
	nhor2(ielem) = ihor2
	nhor3(ielem) = ihor3
	nhor4(ielem) = ihor4

	if(ihor0.gt.25 .OR. ihor1.gt.90 .OR. ihor2.gt.210 .OR. ihor3.gt.350 .OR. ihor4.gt.800) then
	 write(*,*) '25, 90, 210, 350, 800',ihor0,ihor1,ihor2,ihor3,ihor4
	 stop
	end if

end do


!
! Calculating probabilities
!


if(imetodehoritzontalweight .eq. 1) then

	if(distranshormax.lt.distranshor2 .AND. distranshormax.gt.distranshor1) then    ! 400 > and > 200 km (case 4.1 => 200 km)

		do ielem = 1, nelem

			nhor2(ielem) = 0
			nhor3(ielem) = 0
			nhor4(ielem) = 0

		end do

		coefvh100 = 0.408/(0.408 + 0.157)
		coefvh200 = 0.157/(0.408 + 0.157)
		coefvh400 = 0.
		coefvh600 = 0.
		coefvh800 = 0.

	else if(distranshormax.lt.distranshor3) then									 ! < 600 km (cas 8.1 => 400 km)

		do ielem = 1, nelem

			nhor3(ielem) = 0
			nhor4(ielem) = 0

		end do

		coefvh200 = 0.565/(0.565 + 0.217)
		coefvh400 = 0.217/(0.565 + 0.217)
		coefvh600 = 0.
		coefvh800 = 0.
																					! < 800 km (case 12.1 => 600 km)
	else if(distranshormax.lt.distranshor4) then

		do ielem = 1, nelem

			nhor4(ielem) = 0

		end do

		coefvh200 = 0.565/(0.565 + 0.217 + 0.174)
		coefvh400 = 0.217/(0.565 + 0.217 + 0.174)
		coefvh600 = 0.174/(0.565 + 0.217 + 0.174)
		coefvh800 = 0.

	else if(distranshormax.gt.distranshor4) then									 ! > 800 km	(case 16.1 => 800 km)

		coefvh200 = 0.565
		coefvh400 = 0.217
		coefvh600 = 0.174
		coefvh800 = 0.044

	else

		write(*,*) "coefvh = cap"
		stop

	end if

end if

!!!	open output files

open(11,file=outfilen)
open(12,file=outfilep)


! Initialize iterations (1 iteration = 1 generation)


do igen = 1,igenmax 

  write(*,*) 'igen = ',igen

  call lpn_popite

  call lpn_output

end do

!! To totally fill the output matrix

do ielem = 1, nelem

	if(indp(ielem).ne.1) then

		write(12,101) coord(1,ielem),coord(2,ielem),igenmax,varpp1(ielem)

		indp(ielem) = 1

	end if

	if(indn(ielem).ne.1) then

		write(11,101) coord(1,ielem),coord(2,ielem),igenmax,varpn1(ielem)

		indn(ielem) = 1

	end if

end do



101 format(2f16.5,i5,1f16.5) 



!! clouse output files

close(11)
close(12)  

end subroutine lpn_solite


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lpn_popite
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-----------------------------------------------------------------------
!**** lpn_popite
! NAME 
!    lpn_popite
!
! DESCRIPTION
!    Compute lattice iterations
!
!
! USES
!	
!   
!    
!
! USED BY
!	lpn_solite
!
!***
!-----------------------------------------------------------------------

use	def_sim_lpn

implicit none

integer*4 :: ielemorigin,icountorigin,ineigorigin,ineigorigin2destiny

integer*4 :: imetodehoritzontal,jelem
integer*4 :: jelemhgv(8),ielemhgv

real*8 :: contpopn,contpopp,weightorigin
real*8 :: transver,transhor
real*8 :: termfreq


!
! Initialization - constants
!

if (icoeffp.eq.1) then

	do ielem = 1, nelem

		coefhgvisit(ielem) = 0.

		if(varpn1(ielem).gt.0) then
		
			coefhgvisit(ielem) = 1.

		else

			xohor = coord(1,ielem)
			yohor = coord(2,ielem)
		
			jelemhgv(1) = ielem-1
			jelemhgv(2) = ielem+1
			jelemhgv(3) = ielem-180
			jelemhgv(4) = ielem+180
			jelemhgv(5) = ielem-2
			jelemhgv(6) = ielem+2
			jelemhgv(7) = ielem-180-180
			jelemhgv(8) = ielem+180+180
			jelemhgv(9) = ielem-179
			jelemhgv(10)= ielem-181
			jelemhgv(11)= ielem+181
			jelemhgv(12)= ielem+179

			do ielemhgv = 1, 8

				jelem = jelemhgv(ielemhgv)

				xdhor = coord(1,jelem)
				ydhor = coord(2,jelem)
				disthor = sqrt((xdhor-xohor)*(xdhor-xohor) + (ydhor-yohor)*(ydhor-yohor))
				
				if(varpn1(jelem).gt.0. .AND. disthor.lt.distranshor) coefhgvisit(ielem) = 1.

			end do

		end if

	end do

end if


imetodehoritzontal = 2 ! horitzontal method 

!
! 1st step (vertical transmission)
!

do ielem = 1, nelem

	if(varpn1(ielem).gt.0. .AND. gamma.gt.0.) then

		transver = gamma * varpn1(ielem)/(varpp1(ielem)+varpn1(ielem))

		if(transver .lt. 1.) then

			varpn2(ielem) = rtermn * (varpn1(ielem) + transver * varpp1(ielem))

			varpp2(ielem) = rtermp * (varpp1(ielem) - transver * varpp1(ielem))

		else

			varpn2(ielem) = rtermn * (varpn1(ielem) + varpp1(ielem))

			varpp2(ielem) = 0.

		end if		

	else

		varpn2(ielem) = rtermn * varpn1(ielem) 

		varpp2(ielem) = rtermp * varpp1(ielem) 

	end if

	if(varpp2(ielem).lt.1.) varpp2(ielem)=0.
	if(varpn2(ielem).lt.1.) varpn2(ielem)=0.
	if(varpn2(ielem).gt.densnmx) varpn2(ielem) = densnmx
	if(varpp2(ielem).gt.denspmx) varpp2(ielem) = denspmx

end do

!
! 2n step (horizontal / oblique transmission)
!

if(coeff.gt.0.) then

  if(imetodehoritzontal.eq.1) then

	do ielem = 1, nelem

		if(varpn2(ielem).gt.0.) then

			transhor = coeff * varpn2(ielem)/(varpp2(ielem)+varpn2(ielem))

			if(transhor .lt. 1.) then

				varpn1(ielem) = varpn2(ielem) + transhor * varpp2(ielem)

				varpp1(ielem) = varpp2(ielem) - transhor * varpp2(ielem)

			else

				varpn1(ielem) = varpn2(ielem) + varpp2(ielem)

				varpp1(ielem) = 0.

			end if		

		else

			varpn1(ielem) = varpn2(ielem) 

			varpp1(ielem) = varpp2(ielem) 

		end if

		if(varpp1(ielem).lt.1.) varpp1(ielem)=0.
		if(varpn1(ielem).lt.1.) varpn1(ielem)=0.
		if(varpn1(ielem).gt.densnmx) varpn1(ielem) = densnmx
		if(varpp1(ielem).gt.denspmx) varpp1(ielem) = denspmx

	end do

 else if(imetodehoritzontal.eq.2) then

	do ielem = 1, nelem

		varpnhor(ielem) = 0.

	end do

	if(imetodehoritzontalweight.eq.1 .AND. distranshormax.le.4.2) then

		do ielem = 1, nelem

			rnhor0 = real(nhor0(ielem))
		
			do ihor0 = 1, nhor0(ielem)

				jelem = iodhor0(ihor0,ielem)

				varpnhor(jelem) = varpn2(ielem) * coefvh100 / rnhor0 + varpnhor(jelem)

			end do

			rnhor1 = real(nhor1(ielem))

			do ihor1 = 1, nhor1(ielem)

				jelem = iodhor1(ihor1,ielem)

				varpnhor(jelem) = varpn2(ielem) * coefvh200 / rnhor1 + varpnhor(jelem)

			end do

		end do

	else if (imetodehoritzontalweight.eq.1 .AND. distranshormax.gt.4.2) then

		do ielem = 1, nelem

			rnhor1 = real(nhor1(ielem))
			rnhor0 = real(nhor0(ielem))

			do ihor0 = 1, nhor0(ielem)

				jelem = iodhor0(ihor0,ielem)

				varpnhor(jelem) = varpn2(ielem) * coefvh200 / (rnhor1 + rnhor0) + varpnhor(jelem)

			end do

			do ihor1 = 1, nhor1(ielem)

				jelem = iodhor1(ihor1,ielem)

				varpnhor(jelem) = varpn2(ielem) * coefvh200 / (rnhor1 + rnhor0) + varpnhor(jelem)

			end do

  			rnhor2 = real(nhor2(ielem))
			
			do ihor2 = 1, nhor2(ielem)

				jelem = iodhor2(ihor2,ielem)

				varpnhor(jelem) = varpn2(ielem) * coefvh400 / rnhor2 + varpnhor(jelem)

			end do

			rnhor3 = real(nhor3(ielem))
			
			do ihor3 = 1, nhor3(ielem)

				jelem = iodhor3(ihor3,ielem)

				varpnhor(jelem) = varpn2(ielem) * coefvh600 / rnhor3 + varpnhor(jelem)

			end do

			rnhor4 = real(nhor4(ielem))
			
			do ihor4 = 1, nhor4(ielem)

				jelem = iodhor4(ihor4,ielem)

				varpnhor(jelem) = varpn2(ielem) * coefvh800 / rnhor4 + varpnhor(jelem)

			end do

		end do

	else 

		do ielem = 1, nelem

			rnhor4 = real(nhor4(ielem))
		
			do ihor4 = 1, nhor4(ielem)

				jelem = iodhor4(ihor4,ielem)

				varpnhor(jelem) = varpn2(ielem) / rnhor4 + varpnhor(jelem)

			end do

		end do

	end if

	do ielem = 1, nelem

		varpnhor(ielem) = coefpnvisit * varpnhor(ielem) + varpn2(ielem) 

		if(varpnhor(ielem).gt.0.) then


			termfreq = coefalfa*(2.*varpnhor(ielem)/(varpp2(ielem)+varpnhor(ielem)) -1.)
			transhor = (coeff + termfreq) * (varpnhor(ielem) + coefpp * coefhgvisit(ielem) * varpp2(ielem))/(varpp2(ielem)+varpnhor(ielem))

			if(transhor .lt. 1.) then

				varpn1(ielem) = varpn2(ielem) + transhor * varpp2(ielem)

				varpp1(ielem) = varpp2(ielem) - transhor * varpp2(ielem)

			else

				varpn1(ielem) = varpn2(ielem) + varpp2(ielem)

				varpp1(ielem) = 0.

			end if		

		else

			varpn1(ielem) = varpn2(ielem) 

			varpp1(ielem) = varpp2(ielem) 

		end if

		if(varpp1(ielem).lt.1.) varpp1(ielem)=0.
		if(varpn1(ielem).lt.1.) varpn1(ielem)=0.
		if(varpn1(ielem).gt.densnmx) varpn1(ielem) = densnmx
		if(varpp1(ielem).gt.denspmx) varpp1(ielem) = denspmx

	end do

 end if

else 

	do ielem = 1, nelem

		varpn1(ielem) = varpn2(ielem) 
		varpp1(ielem) = varpp2(ielem) 

		if(varpp1(ielem).lt.1.) varpp1(ielem)=0.
		if(varpn1(ielem).lt.1.) varpn1(ielem)=0.
		if(varpn1(ielem).gt.densnmx) varpn1(ielem) = densnmx
		if(varpp1(ielem).gt.denspmx) varpp1(ielem) = denspmx

	end do

end if


!
! 3rd step (Advection step - Streaming)
!

do ielem = 1, nelem

	contpopp = 0.
	contpopn = 0.

	do ineig = 1, 4

		if(lcount(ineig,ielem).gt.0) then

			if(lcount(ineig,ielem) .ge. 500) then

				write(*,*) 'trepitja memoria lcount = ', lcount(ineig,ielem)

				stop

			end if
				
			do icount = 1, lcount(ineig,ielem)

				ielemorigin = lelecont(ineig,ielem,icount)

				ineigorigin2destiny = 0

				do ineigorigin = 1, 4

					do icountorigin = 1, lcount(ineigorigin,ielemorigin)

						if(ielem.eq.lelecont(ineigorigin,ielemorigin,icountorigin)) ineigorigin2destiny = ineigorigin

					end do

				end do

				weightorigin = lcount(ineigorigin2destiny,ielemorigin)

				
				if(weightorigin.ne.0 .AND. weight(ielemorigin).ne.0) then

					contpopn = contpopn + (1.-ctnpenv(ielemorigin))*popdens(ielemorigin)*varpn1(ielemorigin)/(weightorigin*weight(ielemorigin))
	
				end if
			
			end do
		
		end if
	
	end do

	varpn2(ielem) = ctnpenv(ielem)*varpn1(ielem) + contpopn

	do ineig = 1, 4

		jelem = lneig(ineig,ielem)

		contpopp = (1. - ctnpepv(jelem))*varpp1(jelem)*popdens(jelem) + contpopp

	end do

	varpp2(ielem) = ctnpepv(ielem)*varpp1(ielem) + contpopp

	if(varpp2(ielem).lt.1.) varpp2(ielem)=0.
	if(varpn2(ielem).lt.1.) varpn2(ielem)=0.


end do					


!
! Final control
!

do ielem = 1, nelem

	varpn1(ielem) = varpn2(ielem)

	varpp1(ielem) = varpp2(ielem)

	if(varpn1(ielem).gt.densnmx) varpn1(ielem) = densnmx

	if(varpn1(ielem).lt.1.) varpn1(ielem) = 0.

	if(varpp1(ielem).gt.denspmx) varpp1(ielem) = denspmx

	if(varpp1(ielem).lt.1.) varpp1(ielem) = 0.

end do

!
 
end subroutine lpn_popite

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lpn_output
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-----------------------------------------------------------------------
!****f* lpn_output
! NAME 
!    lpn_output
!
! DESCRIPTION
!	 output results
!
! USES
!    
!
! USED BY
!	  lpn_solite
!
!***
!--------------------------------------------------------------------------
!

use	def_sim_lpn

implicit none

do ielem = 1, nelem

	if(varpn1(ielem).gt.0.9*densnmx) then
			
		if(indn(ielem).ne.1) then

			write(11,101) coord(1,ielem),coord(2,ielem),igen

			indn(ielem) = 1

		end if

	end if

end do

do ielem = 1, nelem

	if(varpp1(ielem).lt.0.9*denspmx) then
				
		if(indp(ielem).ne.1) then

			write(12,101) coord(1,ielem),coord(2,ielem),igen

			indp(ielem) = 1

		end if
 
	end if

end do


101 format(2f16.5,i5) 


end subroutine lpn_output
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

