!  subroutine to determine the four flow directions
! for four flow direction system
! fp(i,1) is the northern nodal point of cell i
! fp(i,2) is the eastern nodal point of cell i
! fp(i,3) is the southern nodal point of cell i
! fp(i,4) is the western nodal point of cell i

! for eight flow direction system
! fp(i,1) is the northern nodal point of cell i
! fp(i,2) is the northeastern nodal point of cell i
! fp(i,3) is the eastern nodal point of cell i
! fp(i,4) is the southeastern nodal point of cell i
! fp(i,5) is the southern nodal point of cell i
! fp(i,6) is the southwestern nodal point of cell i
! fp(i,7) is the western nodal point of cell i
! fp(i,8) is the northwestern nodal point of cell i
!  by Chen H.X., HKUST Feb. 2014 latest revison 25 Feb. 2014

	subroutine flodir(grd,pth,ncol,nrow,celsiz,nodat,pf,pf1,ctr,imax,temp,u,infil,param,header,u1,fp,outflow)
	implicit none
	integer grd,pth,i,m,ncol,nrow,ctr,imax,u,u1
	double precision param(6),nodat,celsiz,cns,cew,temp(pth)
	double precision east,west,north,south
	double precision pf(imax),pf1(grd),nodats
	integer fp(imax,8),cellrecord(grd),cellnow
	character*14 header(6)
	character*255 infil
	logical outflow(imax)
	
	cellrecord=0
!  	
	open(u,file=infil,status='old',err=23)
	do 200, m=1,6
	read(u,*) header(m),param(m)
  200	continue
! set default value of nodat & celsiz for use with GIS ascii files   
  	nodat=-9999.d0
  	celsiz=-10.d0
	do 210, m=1,6
	if (trim(header(m)).eq.'ncols') ncol=int(param(m))
	if (trim(header(m)).eq.'nrows') nrow=int(param(m))
	if (trim(header(m)).eq.'cellsize') celsiz=param(m)
	if (trim(header(m)).eq.'NODATA_value') nodat=param(m)
	if (trim(header(m)).eq.'nodata_value') nodat=param(m)
	if (trim(header(m)).eq.'cols:') ncol=int(param(m))
	if (trim(header(m)).eq.'rows:') nrow=int(param(m))
	if (trim(header(m)).eq.'east:') east=param(m)
	if (trim(header(m)).eq.'west:') west=param(m)
	if (trim(header(m)).eq.'north:') north=param(m)
	if (trim(header(m)).eq.'south:') south=param(m)
  210	continue
  	if (celsiz.le.0) then
  	  cew=abs(east-west)/ncol
  	  cns=abs(north-south)/nrow
  	  if (cew.eq.cns) then
  	    celsiz=cew
  	  else
  	    celsiz=sqrt(cew*cns)
  	    write(*,*) 'Rectangular cells ',cew, ' X ', cns
  	    write(u1,*) 'Rectangular cells ',cew, ' X ', cns
  	  end if
  	end if
  	if (ncol*nrow .gt. grd) then
  	  write(*,*) 'Grid file exceeds array size'
	  write (*,*) '--> ',trim(infil)
  	  write(*,*) 'Check intialization file row and column values.'
  	  write(u1,*) 'Grid file exceeds array size'
	  write (u1,*) '--> ',trim(infil)
  	  write(u1,*) 'Check intialization file row and column values.'
	  pause 'Press RETURN to exit'
	  close(u)
	  close(u1)
	  stop
  	end if
  	nodats=nodat
	ctr=0
	do 120, m=1,nrow
!  next sequence of lines read data in but skips no_data values
!  count maintained by ctr should coincide with node numbers from GIS
!  pf1() keeps track of positions of nodata values so that results
!  can be written out in grid format.
  	read(u,*,end=125) (temp(i), i=1,ncol)
	do 250, i=1,ncol
	  pf1(i+(m-1)*ncol)=temp(i)
	  if(temp(i).ne.nodats) then
	    ctr=ctr+1
	    cellrecord(i+(m-1)*ncol)=ctr
	   
	    if (ctr>imax) then
  	      write(*,*) 'Number of data cells exceeds array size'
	      write (*,*) '--> ',trim(infil)
  	      write(*,*) 'Check imax value in intialization file.'
  	      write(u1,*) 'Number of data cells exceeds array size'
	      write (u1,*) '--> ',trim(infil)
  	      write(u1,*) 'Check imax value in intialization file.'
	      pause 'Press RETURN to exit'
  	      close(u)
	      close(u1)
	    end if
	    pf(ctr)=temp(i)
	  end if
  250	continue
  120	continue
  
! determine the four flow directions of each cell
do m=1,nrow
do i=1,ncol
cellnow=cellrecord(i+(m-1)*ncol)
! cellnow should not be zero
if (cellnow/=0) then
    if (m==1) then
        fp(cellnow,1)=0
        fp(cellnow,2)=0
        fp(cellnow,8)=0
    else
        fp(cellnow,1)=cellrecord(i+(m-2)*ncol)
        if (i==1) then
        fp(cellnow,8)=0
        else
        fp(cellnow,8)=cellrecord(i-1+(m-2)*ncol)
        end if
        if (i==ncol) then
        fp(cellnow,2)=0
        else
        fp(cellnow,2)=cellrecord(i+1+(m-2)*ncol)
        end if
    end if

    if (m==nrow) then
        fp(cellnow,4)=0
        fp(cellnow,5)=0
        fp(cellnow,6)=0
    else
        fp(cellnow,5)=cellrecord(i+m*ncol)
        if (i==1) then
        fp(cellnow,6)=0
        else
        fp(cellnow,6)=cellrecord(i-1+m*ncol)
        end if
        if (i==ncol) then
        fp(cellnow,4)=0
        else
        fp(cellnow,4)=cellrecord(i+1+m*ncol)
        end if
    end if

    if (i==1) then
        fp(cellnow,6)=0
        fp(cellnow,7)=0
        fp(cellnow,8)=0
    else
        fp(cellnow,7)=cellrecord(i-1+(m-1)*ncol)
        if (m==1) then
        fp(cellnow,8)=0
        else
        fp(cellnow,8)=cellrecord(i-1+(m-2)*ncol)
        end if
        if (m==nrow) then
        fp(cellnow,6)=0
        else
        fp(cellnow,6)=cellrecord(i-1+m*ncol)
        end if
    end if

    if (i==ncol) then
        fp(cellnow,2)=0
        fp(cellnow,3)=0
        fp(cellnow,4)=0
    else
        fp(cellnow,3)=cellrecord(i+1+(m-1)*ncol)
        if (m==1) then 
        fp(cellnow,2)=0
        else
        fp(cellnow,2)=cellrecord(i+1+(m-2)*ncol)
        end if
        if (m==nrow) then
        fp(cellnow,4)=0
        else
        fp(cellnow,4)=cellrecord(i+1+m*ncol)
        end if
    end if
! if (cellnow/=0)
end if
end do
end do
  
!  	imax=ctr
  125	close(u)
	return
   23	continue
   	write (*,*) '*** Error opening input file ***'
	write (*,*) '--> ',trim(infil)
	write (*,*) 'Check file name and location'
   	write (u1,*) '*** Error opening input file ***'
	write (u1,*) '--> ',trim(infil)
	write (u1,*) 'Check file name and location'
	pause 'Press RETURN to exit'
	close(u)
	close(u1)
	stop '-10'
	end
