module params
implicit none
integer, parameter :: ncols=21601, nrows=10801
real*8, parameter   :: xllcorner=-180.00833333334, yllcorner=-90.008333333335
real*8, parameter   :: cellsize= 0.01666666667
real*4, parameter   :: nodata = -9999
end module

program readetopo
use params
implicit none
!
real*4, allocatable :: height(:,:)
real*8              :: posx, posy
integer             :: i, j
character(len=80)   :: filename

! for histgram
real, parameter     :: hmin=-10000, hmax=10000
integer, parameter  :: nclass=400
integer :: counter(nclass), countout(2)
real    :: class(nclass, 2)
real :: max, min
!
filename='../Binary/etopo1_ice_g.flt'
allocate( height(1:ncols, 1:nrows) )

do i=1, nclass
 class(i,1) = hmin + (i-1)*(hmax-hmin)/nclass
 class(i,2) = hmin + (i  )*(hmax-hmin)/nclass
end do

open(file=filename, unit=10, status='old', &
	access='direct', recl=4*ncols, form='unformatted')
do j=1, nrows
  read(10, rec=j) (height(i, j), i=1, ncols)
end do

close(10)
call find_maxmin( height, ncols, nrows, max, min)
write(0,*) 'height max =', max, 'height min=', min

call mkhist(height, ncols, nrows, counter, nclass, countout, hmin, hmax)

do i=1, nclass
 write(*,*) (class(i,1)+ class(i,2))/2, counter(i)
end do
end

subroutine  find_maxmin( height, ncols, nrows, max, min )
implicit none
integer :: ncols, nrows
real*4  :: height(ncols, nrows), max, min

integer:: i, j

max = height(1,1)
min = height(1,1)

do i=1, ncols
 do j=1, nrows
  if( height(i,j) < min ) then
   min=height(i,j)
  end if
  if(height(i,j) > max ) then
   max = height(i,j)
  end if
 end do
end do

end


subroutine mkhist(height, ncols, nrows, counter, nclass, &
					countout, hmin, hmax)
implicit none
integer :: ncols, nrows, nclass, counter(nclass), countout(2)
real :: height(ncols, nrows)
real :: hmax, hmin
integer :: i,j,k
real :: dh

counter = 0
countout = 0

dh = (hmax-hmin)/nclass

do i=1, ncols
  do j=1, nrows
   k = (height(i,j)- hmin)/dh+1
   if(k>= 1 .and. k<= nclass ) then
    counter(k) = counter(k) + 1
   else if(k<0) then
    countout(1) = countout(1)+1
   else
    countout(2) = countout(2)+1
   end if
  end do
end do

end

subroutine position( i, j, posx, posy )
use params
implicit none
integer :: i,j
real*8  :: posx, posy

posx = xllcorner+(i-1)*cellsize+cellsize/2
posy = yllcorner+(j-1)*cellsize+cellsize/2

end
