c	The Program can code the drainage basin by Pfafstetter coding system (3 levels).
c	Input: acctxt, dirtxt,(dem,dis)
c	Output: basin0, basin1, basin2, basin3
c	output 2: mainstem grid, distance-altitue graph
c	Contact with: (TANG, Qiuhong) tangqh@iis.u-tokyo.ac.jp


      character*5 ncols,nrows
      character*14 xllcorner,yllcorner
      character*14 cellsize
      character*14 nodata
      integer x0,y0
      integer s, znodata,nc,nr
      integer acc(5000,5000), dir(5000,5000)
      integer maxtmp, outi, outj, max1i,max1j,surround(8,3),max1
      integer suboutlet(4,5),mintmp
      integer suboutlet2(4,5),suboutlet3(4,5)
      integer outlet(5000,2),id,basin(5000,5000)
      integer basin2(5000,5000),basin3(5000,5000)
      integer i,j,k,m,n,tmpi,tmpj
      integer outi2, outj2, stopi2, stopj2, avoidi3, avoidj3
      integer stopi3, stopj3,outi3, outj3
      common /acc/dir
	integer dem(500,500),dis(500,500),ooi,ooj
	real	rivlen
	common dem,dis,ooi,ooj,rivlen

c	river length (0.0: not use, 5464. km for yellow river)      
	rivlen = 5464.
c      read input file (dem,dis) acc and dir file from ArcGIS
      open(1, file= 'wsdem.asc' ,status='old')
      read(1,'(A5, i18)') ncols,nc
      read(1,'(A5, i18)') nrows,nr
      read(1,'(A14, i)') xllcorner,x0
      read(1,'(A14, i)') yllcorner,y0
      read(1,'(A14, i)') cellsize,s
      read(1,'(A14, i)') nodata,znodata
      do i=1,nr        
      read(1,*) (dem(i,j), j=1,nc)
      end do           
      close (1)

      open(1, file= 'dis.asc' ,status='old')
      read(1,'(A5, i18)') ncols,nc
      read(1,'(A5, i18)') nrows,nr
      read(1,'(A14, i)') xllcorner,x0
      read(1,'(A14, i)') yllcorner,y0
      read(1,'(A14, i)') cellsize,s
      read(1,'(A14, i)') nodata,znodata
      do i=1,nr        
      read(1,*) (dis(i,j), j=1,nc)
      end do           
      close (1)
		      	      
      open(1, file= 'acc.asc' ,status='old')
      read(1,'(A5, i18)') ncols,nc
      read(1,'(A5, i18)') nrows,nr
      read(1,'(A14, i)') xllcorner,x0
      read(1,'(A14, i)') yllcorner,y0
      read(1,'(A14, i)') cellsize,s
      read(1,'(A14, i)') nodata,znodata
      do i=1,nr        
      read(1,*) (acc(i,j), j=1,nc)
      end do           
      close (1)
      
      open(1, file= 'dir.asc' ,status='old')
      read(1,'(A5, i18)') ncols,nc
      read(1,'(A5, i18)') nrows,nr
      read(1,'(A14, i)') xllcorner,x0
      read(1,'(A14, i)') yllcorner,y0
      read(1,'(A14, i)') cellsize,s
      read(1,'(A14, i)') nodata,znodata
      do i=1,nr        
      read(1,*) (dir(i,j), j=1,nc)
      end do           
      close (1)     
      
c       search the outlet (Max of accumulation)
      maxtmp=0
      do i=1,nr
      	do j=1,nc
      		if (acc(i,j).gt.maxtmp) then
      		maxtmp= acc (i,j)
      		outi=i
      		outj=j
      		end if
      	end do
      end do
	ooi=outi
	ooj=outj
      print *, 'nc,nr,x0,y0,size',nc,nr,x0,y0,s
      print *, 'outlet acc =', maxtmp, ' outi=',outi, ' outj=', outj
      
      id=1
      outlet(1,1)=outi
      outlet(1,2)=outj
      do while (upgrids(acc,dir,basin,outlet,id).ne.1)
      end do
      
      open(1,file='basin0',status='unknown')
      write(1,'(A, i)') ncols,nc
      write(1,'(A, i)') nrows,nr
      write(1,'(A, i)') xllcorner,x0
      write(1,'(A, i)') yllcorner,y0
      write(1,'(A, i)') cellsize,s
      write(1,'(A, i)') nodata,0      
      do i=1,nr
      	do j=1,nc
      	write (1,'(i5,$)') basin(i,j)
      	basin(i,j) = 0
      	end do
      write (1,*)
      end do
      close(1)
      
      do i=1,5000
      	outlet(i,1)=0
      	outlet(i,2)=0
      end do      
c     search for the the four largest tributary (intersections with main stem) 
c     level 1

c     get outlet of sub-basin (in array suboutlet, no order)      
      call getsubout(acc,dir,outi,outj,-1,-1,suboutlet)
c	  according to suboutlet, search sub-basin 1-9    
      call getsub (acc,dir,outi,outj,-1,-1,-1,-1,suboutlet,basin) 
c	  output to ascii file      
      open(1,file='basin1',status='unknown')
      write(1,'(A, i)') ncols,nc
      write(1,'(A, i)') nrows,nr
      write(1,'(A, i)') xllcorner,x0
      write(1,'(A, i)') yllcorner,y0
      write(1,'(A, i)') cellsize,s
      write(1,'(A, i)') nodata,0        
      do i=1,nr
      	do j=1,nc
      	write (1,'(i5,$)') basin(i,j)
c      	basin(i,j) = 0
      	end do
      write (1,*)
      end do
      close(1)
      
      do i=1,5000
      	outlet(i,1)=0
      	outlet(i,2)=0
      end do  
      
      print *, 'level 1 OK!'  
            
c     level 2
	  do i=1,9
	  	call getsub2(acc,dir,outi,outj,-1,-1,-1,-1,
     *	  	suboutlet,suboutlet2,basin2,i)
c		level 3
		if (i.le.4) then
			outi3=suboutlet(i,2)
			outj3=suboutlet(i,3)
			stopi3=-1
			stopj3=-1
			avoidi3=-1
			avoidj3=-1
		else
				k=i-4
			if (k.eq.1) then
				outi3=outi
				outj3=outj
				stopi3=suboutlet(k,4)
				stopj3=suboutlet(k,5)	
				avoidi3=-1
				avoidj3=-1				
			else if (k.eq.5) then
				outi3=suboutlet(k-1,4)
				outj3=suboutlet(k-1,5)
				stopi3=-1
				stopj3=-1
				avoidi3=suboutlet(k-1,2)
				avoidj3=suboutlet(k-1,3)											
			else
		  		outi3=suboutlet(k-1,4)
		  		outj3=suboutlet(k-1,5)
		  		stopi3=suboutlet(k,4)
		  		stopj3=suboutlet(k,5)	
				avoidi3=suboutlet(k-1,2)
				avoidj3=suboutlet(k-1,3)						  				
			end if		
		end if	
		do j=1,9
c		if (i.eq.5.and.j.eq.7) then	
c		print *, 'input: ',outi3,outj3,stopi3,stopj3,avoidi3,avoidj3		
			call getsub2(acc,dir,outi3,outj3,stopi3,stopj3,avoidi3,avoidj3,
     *		suboutlet2,suboutlet3,basin3,j)
c     	end if     
		end do				 	
	  end do
	  
c	  output to ascii file      
      open(1,file='basin2',status='unknown')
      open(2,file='basin3',status='unknown')
      write(1,'(A, i)') ncols,nc
      write(1,'(A, i)') nrows,nr
      write(1,'(A, i)') xllcorner,x0
      write(1,'(A, i)') yllcorner,y0
      write(1,'(A, i)') cellsize,s
      write(1,'(A, i)') nodata,0  
      
      write(2,'(A, i)') ncols,nc
      write(2,'(A, i)') nrows,nr
      write(2,'(A, i)') xllcorner,x0
      write(2,'(A, i)') yllcorner,y0
      write(2,'(A, i)') cellsize,s
      write(2,'(A, i)') nodata,0             
      do i=1,nr
      	do j=1,nc
      	write (1,'(i5,$)') basin(i,j)*10+basin2(i,j)
      	if (basin3(i,j).eq.-1) then
      	write (2,'(i5,$)') basin(i,j)*100+basin2(i,j)*10+0
      	else
      	write (2,'(i5,$)') basin(i,j)*100+basin2(i,j)*10+basin3(i,j)
      	end if
      	basin(i,j) = 0
      	end do
      write (1,*)
      write (2,*) 
      end do
      close(1)
      close(2)
      
      do i=1,5000
      	outlet(i,1)=0
      	outlet(i,2)=0
      end do  
      
      print *, 'level 2,3 OK!'    
             
         
      end
      
c	  get level 2 sub- basins
	  subroutine getsub2(acc,dir,outi,outj,stopi,stopj,avoidi,avoidj,
     * suboutlet,suboutlet2,basin2,j)
      integer avoidi,avoidj
	  integer outi,outj,outi2,outj2,stopi,stopj,avoidi3,avoidj3
	  integer suboutlet(4,5),suboutlet2(4,5)
	  integer basin2(5000,5000)
      integer acc(5000,5000), dir(5000,5000)
      integer i,j, stopi2,stopj2
      i=j	  
c		for sub-basin		4 sub-basins	  
	  if (i.le.4) then
		outi2=suboutlet(i,2)
		outj2=suboutlet(i,3)
		call getsubout(acc,dir,outi2,outj2,-1,-1,suboutlet2)		
		call getsub (acc,dir,outi2,outj2,-1,-1,-1,-1,suboutlet2,basin2)	
	  else
	  i=i-4
c		for inter-basin	  5 inter-basins
	  	if (i.eq.1) then
	  		outi2=outi
	  		outj2=outj
	  		stopi2=suboutlet(i,4)
	  		stopj2=suboutlet(i,5)
	  		avoidi3=avoidi
	  		avoidj3=avoidj	  		
	  	else if (i.eq.5) then
	  		outi2=suboutlet(i-1,4)
	  		outj2=suboutlet(i-1,5)
	  		stopi2=stopi
	  		stopj2=stopj 
	  		avoidi3=suboutlet(i-1,2)
	  		avoidj3=suboutlet(i-1,3)	  		 	
	  	else
	  		outi2=suboutlet(i-1,4)
	  		outj2=suboutlet(i-1,5)
	  		stopi2=suboutlet(i,4)
	  		stopj2=suboutlet(i,5)
	  		avoidi3=suboutlet(i-1,2)
	  		avoidj3=suboutlet(i-1,3)		  				  	
	  	end if 		  	
		call getsubout(acc,dir,outi2,outj2,stopi2,stopj2,suboutlet2)
		call getsub (acc,dir,outi2,outj2,avoidi3,avoidj3
     *		,stopi2,stopj2,suboutlet2,basin2)
	  end if
	  
	  end	        
      
c	  get 9 sub-basin =4 sub-basin + 5 inter-basin, saved in array
      subroutine getsub(acc,dir,outi,outj,avoidi,avoidj
     *     ,stopi,stopj,suboutlet,basin)
      integer outi, outj
	  integer suboutlet(4,5),surround(8,3)      
      integer basin(5000,5000)
      integer acc(5000,5000), dir(5000,5000)  
      integer max1,tmpi,tmpj  
      integer i,j,k,m,n,id
      integer outlet(5000,2)
      integer avoidi,avoidj
      integer stopi,stopj,space
      integer id0
      
      tmpi=outi
      tmpj=outj
      max1=1
      id0=1
      do i=1, 4
      	do j=2,5
	      	if (suboutlet(i,j).eq.0.and.id0.eq.1) then
	      	id0=0      	
	      	end if
      	end do
      end do      
      if (id0.eq.0) then
      
      do while (max1.ge.1.and.
     *      (.not.(max1i.eq.stopi.and.max1j.eq.stopj)))      	
      	call record(acc,dir,tmpi,tmpj,max1,max1i,max1j,surround)   
         do i=1,8      	    		
		  if (.not.(surround(i,2).eq.avoidi.and.surround(i,3).eq.avoidj)) then 
		    do m=1,5000
		      	outlet(m,1)=0
		      	outlet(m,2)=0
		    end do 			   		
		    outlet(1,1)=surround(i,2)
		    outlet(1,2)=surround(i,3)			    
		    do while (upgrids(acc,dir,basin,outlet,-1).ne.1)
		    end do  
		  end if  
		    if (.not.(tmpi.eq.stopi.and.tmpj.eq.stopj)) then
		  		basin(tmpi,tmpj)=-1		  		
		  	end if		  	  
			tmpi=max1i
	      	tmpj=max1j  	  
         end do
	  end do 	      		    
      
      else
	  
      id=1     
      do while (max1.ge.1.and.
     *      (.not.(max1i.eq.stopi.and.max1j.eq.stopj)))      	
      	call record(acc,dir,tmpi,tmpj,max1,max1i,max1j,surround)         	   	
      	k=0
      	do i=1, 8
      		do j= 1,4
      		    if (surround(i,2).eq.suboutlet(j,2).and.
     *	surround(i,3).eq.suboutlet(j,3)) then
					k=j
					n=i
        	    end if        	    
      		end do
      	end do       
		   if (k.ne.0) then
c	it is a section of main stem and tributary		   
			id=id+1
			    do m=1,5000
			      	outlet(m,1)=0
			      	outlet(m,2)=0
			    end do 		
		      	outlet(1,1)=suboutlet(k,2)
		      	outlet(1,2)=suboutlet(k,3)
		      	do while (upgrids(acc,dir,basin,outlet,id).ne.1)  
		      	end do   
		    id=id+1
		      	do i=1, 8
		if (i.ne.n.and.surround(i,2).ne.0.and.surround(i,3).ne.0) then
			outlet(1,1)=surround(i,2)
			outlet(1,2)=surround(i,3)
			do while (upgrids(acc,dir,basin,outlet,id).ne.1)  
			end do 		
		end if
		    	end do			      	   
		   else	
c	it is main stem		   
		   	do i=1, 8
				if (surround(i,2).ne.0.and.surround(i,3).ne.0) then   				
	  if (.not.(surround(i,2).eq.avoidi.and.surround(i,3).eq.avoidj)) then 
	    do m=1,5000
	      	outlet(m,1)=0
	      	outlet(m,2)=0
	    end do 			   		
	    outlet(1,1)=surround(i,2)
	    outlet(1,2)=surround(i,3)				    
	    do while (upgrids(acc,dir,basin,outlet,id).ne.1)
	    end do  
	  end if      		
				end if
			end do
		   end if  
		if (.not.(tmpi.eq.stopi.and.tmpj.eq.stopj)) basin(tmpi,tmpj)=id
		tmpi=max1i
      	tmpj=max1j    	
      end do
c	deal with last grid of the basin  
	   if (.not.(tmpi.eq.stopi.and.tmpj.eq.stopj)) basin(tmpi,tmpj)=id             
       
       end if
      end
      
c     get outlet of 4 sub - basin x ,y
      subroutine getsubout(acc,dir,outi,outj,stopi,stopj,suboutlet)
      integer outi, outj, max2, mintmp
      integer suboutlet(4,5),surround(8,3)
      integer acc(5000,5000), dir(5000,5000),mainstem(126,221)  
      integer tmpi,tmpj  
	integer dem(500,500),dis(500,500),ooi,ooj
	real	rivlen
	common dem,dis,ooi,ooj,rivlen
c	  stopi, stopj is for inter-basin search, 
c     when meeting the end of interbasin, stop search     
      integer stopi,stopj
      integer i,j,k,m,n,l
      real	klen
	
	maxlen= -1.0
	do i=1,500
		do j=1,500
			if (dis(i,j).gt.maxlen) maxlen=dis(i,j)
		enddo
	enddo

	if (rivlen.ne.0.) klen=rivlen/maxlen
	if (rivlen.eq.0.) klen=1.0
      print *,'klen=',klen*1000.
	mainstem=0

	open (100,file='dis_alt.txt',status='unknown')

      do i=1,4
      	do j=1,5
      		suboutlet(i,j)=0
      	end do
      end do 
      tmpi=outi
      tmpj=outj      
      space=0      
      call record(acc,dir,tmpi,tmpj,max1,max1i,max1j,surround)
	if(outi.eq.ooi.and.outj.eq.ooj.and.stopi.eq.-1.and.stopj.eq.-1)then
		mainstem(tmpi,tmpj)=1
		mainstem(max1i,max1j)=1
		write(100,*) dis(tmpi,tmpj)*klen,dem(tmpi,tmpj)
		write(100,*) dis(max1i,max1j)*klen,dem(max1i,max1j)
	endif
      do while (max1.ge.1.and.
     * (.not.(max1i.eq.stopi.and.max1j.eq.stopj)))      
      	tmpi=max1i
      	tmpj=max1j    
      	space=space+1	
      call record(acc,dir,tmpi,tmpj,max1,max1i,max1j,surround) 
	if(outi.eq.ooi.and.outj.eq.ooj.and.stopi.eq.-1.and.stopj.eq.-1)then
		mainstem(max1i,max1j)=1
		write(100,*) dis(max1i,max1j)*klen,dem(max1i,max1j)
	endif
      mintmp=min(suboutlet(1,1),suboutlet(2,1),suboutlet(3,1),
     * suboutlet(4,1))
		max2 = -1
		do i=1,8
			if (surround(i,1).gt.max2) then
			max2= surround(i,1)
			k=i
			end if
		end do 		
c give the x,y of suboutlet, 1-4 order: from down to upstream		
		if (max2.ge.mintmp.and.surround(k,2).ne.0) then
			l=0
			do i=1,4 
				if (suboutlet(i,1).eq.mintmp.and.l.eq.0) then
					j=i
					l=l+1
				end if
			end do
			if (j.ne.4) then
				do m=j,3
					suboutlet(m,1)=suboutlet(m+1,1)
					suboutlet(m,2)=suboutlet(m+1,2)
					suboutlet(m,3)=suboutlet(m+1,3)
					suboutlet(m,4)=suboutlet(m+1,4)
					suboutlet(m,5)=suboutlet(m+1,5)				
				end do			
			end if
			suboutlet(4,1)=max2
			suboutlet(4,2)=surround(k,2)
			suboutlet(4,3)=surround(k,3)
			suboutlet(4,4)=tmpi
			suboutlet(4,5)=tmpj	
		end if
      end do      

	if(outi.eq.ooi.and.outj.eq.ooj.and.stopi.eq.-1.and.stopj.eq.-1)then
c	  output to ascii file      
      open(1,file='mainstem',status='unknown')
      write(1,'(A, i)') 'ncols',221
      write(1,'(A, i)') 'nrows',126
      write(1,'(A, i)') 'xllcorner',-457500
      write(1,'(A, i)') 'yllcorner',-1500500
      write(1,'(A, i)') 'cellsize',10000
      write(1,'(A, i)') 'NODATA_value',0        
      do i=1,126
      	do j=1,221
      	write (1,'(i5,$)') mainstem(i,j)
c      	basin(i,j) = 0
      	end do
      write (1,*)
      end do
      close(1)
	endif
	close(100) 
	          
      end
            
c     record the grid in main stem and largest teibutary 
      subroutine record(acc,dir,maini,mainj,
     * maxtmp1,max1i,max1j,tmp)
      	integer maini,mainj,maxtmp1,max1i,max1j
      	integer acctmp,gridvalue
      	integer acc(5000,5000), dir(5000,5000)
      	integer tmp(8,3),i,j,k
            	
c	check the acc file
	do i=1,8
	do j=1,3
		tmp(i,j)=0
	end do	
	end do
    
      i=1      
      if (dir(maini,mainj-1).eq.1) then
      	tmp(i,1)=acc(maini,mainj-1)
      	tmp(i,2)=maini
      	tmp(i,3)=mainj-1
      	i=i+1
      end if
      if (dir(maini-1,mainj-1).eq.2) then
      	tmp(i,1)=acc(maini-1,mainj-1)
      	tmp(i,2)=maini-1
      	tmp(i,3)=mainj-1 
      	i=i+1     
      end if      
      if (dir(maini-1,mainj).eq.4) then
      	tmp(i,1)=acc(maini-1,mainj)
      	tmp(i,2)=maini-1
      	tmp(i,3)=mainj 
      	i=i+1 
      end if 
      if (dir(maini-1,mainj+1).eq.8) then
      	tmp(i,1)=acc(maini-1,mainj+1)
      	tmp(i,2)=maini-1
      	tmp(i,3)=mainj+1
      	i=i+1 
      end if       
      if (dir(maini,mainj+1).eq.16) then
      	tmp(i,1)=acc(maini,mainj+1)
      	tmp(i,2)=maini
      	tmp(i,3)=mainj+1
      	i=i+1 
      end if       
      if (dir(maini+1,mainj+1).eq.32) then
      	tmp(i,1)=acc(maini+1,mainj+1)
      	tmp(i,2)=maini+1
      	tmp(i,3)=mainj+1
      	i=i+1 
      end if       
      if (dir(maini+1,mainj).eq.64) then
      	tmp(i,1)=acc(maini+1,mainj)
      	tmp(i,2)=maini+1
      	tmp(i,3)=mainj
      	i=i+1 
      end if        
      if (dir(maini+1,mainj-1).eq.128) then
      	tmp(i,1)=acc(maini+1,mainj-1)
      	tmp(i,2)=maini+1
      	tmp(i,3)=mainj-1
      	i=i+1 
      end if 
      
      maxtmp1 =-1
      acctmp =0   
      do j=1,i-1
      	acctmp=acctmp+tmp(j,1)+1
      	if (tmp(j,1).gt.maxtmp1) then
      		maxtmp1=tmp(j,1)
      		max1i=tmp(j,2)
      		max1j=tmp(j,3)
      	end if      	
      end do    
      k=0
      do j=1,i-1
      	if (tmp(j,1).eq.maxtmp1.and.k.eq.0) then      	
      		tmp(j,1)=0
      		tmp(j,2)=0
      		tmp(j,3)=0
      		k=1
      	end if      	
      end do 
             
      if (acctmp.ne.acc(maini,mainj)) then
      	print *, maini, mainj, acctmp, ' error!'
      	stop
      end if
      
      return
      end      
      
c     delineate the river basin that water flows to the given point
      FUNCTION upgrids(acc,dir,basin,outlet,id)
      	integer outlet(5000,2),id
      	integer acc(5000,5000), dir(5000,5000), basin(5000,5000)
      	integer tmpij(5000,2)
      	integer maini,mainj
      	
      do i=1, 5000
      	do j=1,2
      	tmpij(i,j)=0
      	end do
      end do    
        	
      j=1
      do i=1, 5000
      	if  (outlet(i,1).gt.0.and.outlet(i,1).gt.0) then    
	    basin(outlet(i,1),outlet(i,2))=id            	
c	check the dir file
		maini=outlet(i,1)
		mainj=outlet(i,2)
	      if (dir(maini,mainj-1).eq.1) then
		basin(maini,mainj-1)=id
		tmpij(j,1)=maini
		tmpij(j,2)=mainj-1
		j=j+1
	      end if
	      if (dir(maini-1,mainj-1).eq.2) then
		basin(maini-1,mainj-1)=id
		tmpij(j,1)=maini-1
		tmpij(j,2)=mainj-1
		j=j+1
	      end if      
	      if (dir(maini-1,mainj).eq.4) then
		basin(maini-1,mainj)=id
		tmpij(j,1)=maini-1
		tmpij(j,2)=mainj
		j=j+1
	      end if 
	      if (dir(maini-1,mainj+1).eq.8) then
		basin(maini-1,mainj+1)=id
		tmpij(j,1)=maini-1
		tmpij(j,2)=mainj+1
		j=j+1
	      end if       
	      if (dir(maini,mainj+1).eq.16) then
		basin(maini,mainj+1)=id
		tmpij(j,1)=maini
		tmpij(j,2)=mainj+1
		j=j+1
	      end if       
	      if (dir(maini+1,mainj+1).eq.32) then
		basin(maini+1,mainj+1)=id
		tmpij(j,1)=maini+1
		tmpij(j,2)=mainj+1
		j=j+1
	      end if       
	      if (dir(maini+1,mainj).eq.64) then
		basin(maini+1,mainj)=id
		tmpij(j,1)=maini+1
		tmpij(j,2)=mainj
		j=j+1
	      end if        
	      if (dir(maini+1,mainj-1).eq.128) then
		basin(maini+1,mainj-1)=id
		tmpij(j,1)=maini+1
		tmpij(j,2)=mainj-1
		j=j+1
	      end if 
	end if	      
      end do
      
      do i=1, 5000
      	do k=1,2
      	outlet(i,k)=tmpij(i,k)
      	end do
      end do  
c      print *,k
      upgrids=j
      
      return
      end