program landuse
implicit none

integer, parameter :: jptype = 5, jplb1 = 30, jplb2 = 400, jpsup = 40, jplf = 1600000

character :: filename*256
real      :: zbs2(jplb2), zsup(jpsup), fractions(jplf,jptype)
integer   :: ibs1(jplb1), ibs2(jplb2)
integer   :: ibpw = 0, ibpwio = 0, imdi
integer   :: landtype(jplf)
integer   :: i, j, jtype, ilf, ierr, nlon, nlat

call get_command_argument(1, filename)

call ASIMHM(ibpw, ibpwio, imdi)

call ASIMHO(1, filename, 0, ierr)
if (ierr /= 0) then
   write(0,*) 'LANDUSE: CANNOT OPEN ',trim(filename),' ASIMHO RETURNS:', ierr
   call abort
endif

do jtype = 1, jptype        ! 1: Water; 2: Ice; 3: Bare soil; 4: Low vegetation; 5: Forest.

   ibs1 = imdi
   ibs1(9:11) = [194, 105, 900+jtype]

   ilf = jplf
   call GETFD(1, ibs1, jplb1, ibs2, jplb2, zbs2, jplb2, zsup, jpsup, fractions(:,jtype), ilf, 0, ierr)
   if (ierr < 0) then
      write(0,*) 'LANDUSE: ERROR IN GETFD; IERR=',ierr
      write(0,'(30(1x,i3))') ibs1
      call abort
   endif

enddo

call ASIMHC(1, ierr)

! The largest fraction in the array of 5 (along the second dimension)
! determines the (overriding) type of land use.

landtype(1:ilf) = maxloc(fractions(1:ilf,:), dim=2)

nlon = ibs2(7)
nlat = ibs2(9)

do i = 1, nlat
   do j = 1, nlon
      print*, j, i, landtype((i-1)*nlon+j)
   enddo
   print*
enddo

end program landuse