!  Program Name:
!  Author(s)/Contact(s):
!  Abstract:
!  History Log:
! 
!  Usage:
!  Parameters: <Specify typical arguments passed>
!  Input Files:
!        <list file names and briefly describe the data they include>
!  Output Files:
!        <list file names and briefly describe the information they include>
! 
!  Condition codes:
!        <list exit condition or error codes returned >
!        If appropriate, descriptive troubleshooting instructions or
!        likely causes for failures could be mentioned here with the
!        appropriate error code
! 
!  User controllable options: <if applicable>

! This subrouting includs the data structure and tools used for NHDPlus network mapping.
module module_UDMAP

  use config_base, only: nlst
#ifdef MPP_LAND
  use module_mpp_land, only: my_id, local_startx_rt, local_starty_rt,  &
       local_endx_rt,local_endy_rt, left_id, right_id, down_id, up_id, mpp_collect_1d_int_mem, &
       IO_id , numprocs
  use module_mpp_land, only: mpp_land_bcast_int, mpp_land_bcast_real8_1d, mpp_land_bcast_int1, mpp_land_bcast_int8
  
  use module_mpp_land, only: sum_int1d, global_rt_nx, global_rt_ny, write_IO_rt_int, MPP_LAND_COM_INTEGER
  
  use MODULE_mpp_ReachLS, only : updatelinkv, ReachLS_write_io, com_write1dInt, &
       com_decomp1dInt, pack_decomp_int, pack_decomp_real8, &
       com_decomp1dint8, pack_decomp_int8, com_write1dInt8
  use module_hydro_stop, only:HYDRO_stop
  use hashtable

  use iso_fortran_env, only: int64
#endif

  implicit none

#ifndef MPP_LAND
    integer, parameter :: numprocs=1
#endif

#include <netcdf.inc>

type userDefineMapping
    integer, allocatable, dimension(:) :: grid_i, grid_j
    real, allocatable, dimension(:) :: weight, nodeArea, cellArea
    integer :: ngrids
    integer(kind=int64) :: myid
!   for bucket model definition
    real, allocatable, dimension(:) :: cellWeight
    integer, allocatable, dimension(:) :: cell_i, cell_j
    integer :: ncell
end type userDefineMapping

TYPE ( userDefineMapping ), allocatable, DIMENSION (:) :: LUDRSL

integer(kind=int64), allocatable, dimension(:) :: bufid
real*8 , allocatable, dimension(:) :: bufw
integer :: LNUMRSL  ! number of local links
integer :: ter_rt_flag
real*8, allocatable, dimension(:) :: basns_area
integer :: gnpid, lnsize
integer, allocatable, dimension(:) :: bufi,bufj

contains
    subroutine UDMP_ini(nlinksl,ixrt,jxrt,rtmask, OVRTSWCRT, SUBRTSWCRT,cell_area)
!This is the driver for user defined mapping file funciton application.
        integer :: ixrt, jxrt, OVRTSWCRT, SUBRTSWCRT, nlinksl
        integer, intent(in), dimension(ixrt,jxrt):: rtmask
        integer :: npid    !local variable.
        real,dimension(:,:) :: cell_area
        ter_rt_flag = 0
        if(OVRTSWCRT .eq. 1 .or. SUBRTSWCRT .eq. 1) then
            ter_rt_flag = 1
        endif
        call readUDMP(ixrt,jxrt,npid,nlinksl)
        call UDMP2LOCAL(npid,ixrt,jxrt,rtmask, ter_rt_flag)
        call getUDMP_area(cell_area)
    end subroutine UDMP_ini

    subroutine readUDMP(ixrt,jxrt,npid, nlinksl)
        implicit none
        integer :: i,j,Ndata, did, Npid, nlinksl, k, m, kk
        integer,allocatable,dimension(:) ::  nprocs_map, lnsizes, istart
        integer(kind=int64), allocatable, dimension(:) :: g1bufid, gbufid, linkid, bufid_tmp, bufidflag
        integer :: ix_bufid, ii, ixrt,jxrt
        integer, allocatable, dimension(:) :: gbufi,gbufj,bufsize
        real*8 , allocatable, dimension(:) :: gbufw
        
        did = 1
        call get_dimension(trim(nlst(did)%UDMAP_FILE), ndata, npid) 

#ifdef MPP_LAND
        gnpid = npid
        allocate (lnsizes(numprocs))
        if(my_id .eq. io_id) then
           allocate (istart(numprocs))
           allocate (nprocs_map(ndata))
           allocate(gbufi(ndata))
           allocate(gbufj(ndata))
           call get1d_int(trim(nlst(did)%UDMAP_FILE),"i_index",gbufi)
           call get1d_int(trim(nlst(did)%UDMAP_FILE),"j_index",gbufj)
        endif
           call get_nprocs_map(ixrt,jxrt,gbufi,gbufj,nprocs_map,ndata)

        if(my_id .eq. io_id) then
           lnsizes = 0 
           do i =1 , ndata
               if(nprocs_map(i) .gt. 0) then
                  lnsizes(nprocs_map(i)) = lnsizes(nprocs_map(i)) + 1 
               endif
           enddo
        endif
        call mpp_land_bcast_int(numprocs,lnsizes)

     if(my_id .eq. io_id ) then
        kk = 0
        do i = 1, numprocs 
           kk = kk + lnsizes(i) 
        end do
     end if

      if(my_id .eq. IO_id) then
          ii = 1
          do i = 1, numprocs
             istart(i) = ii
             if(lnsizes(i) .gt. 0) then
                ii = lnsizes(i) + ii
             else
                istart(i) = -999
             endif
          end do
      endif

      if(lnsizes(my_id+1) .gt. 0)  allocate(bufi(lnsizes(my_id+1) ))
      call pack_decomp_int(gbufi, ndata, nprocs_map, lnsizes, istart,bufi)
      if(my_id .eq. io_id) then 
           if(allocated(gbufi))  deallocate(gbufi)
      endif

      
      if(lnsizes(my_id+1) .gt. 0) allocate(bufj(lnsizes(my_id+1) ))
      call pack_decomp_int(gbufj, ndata, nprocs_map, lnsizes, istart,bufj)
      if(my_id .eq. io_id)  then 
         if(allocated(gbufj)) deallocate(gbufj)
      endif


! check bufid
!      check  polyid and linkid
        allocate(linkid(nlinksl))
        if(my_id .eq. io_id) then
            call get1d_int64(trim(nlst(did)%route_link_f),"link",linkid)
            allocate(gbufid(npid))
            call get1d_int64(trim(nlst(did)%UDMAP_FILE),"polyid",gbufid)
        endif
#ifdef MPP_LAND
       if(nlinksl .gt. 0) then
          call mpp_land_bcast_int8(nlinksl,linkid)
       endif
       call com_decomp1dInt8(gbufid,npid,bufid_tmp,ix_bufid)
#endif
       if(ix_bufid .gt. 0) then
          allocate(bufidflag(ix_bufid))
          bufidflag = -999
       endif

       ! The following loops are replaced by a hashtable-based algorithm
       !  do i = 1, ix_bufid
       !           do j = 1, nlinksl
       !                if(bufid_tmp(i) .eq. linkid(j)) then
       !                   bufidflag(i) = bufid_tmp(i)
       !                   goto 102
       !                endif
       !           end do
       ! 102       continue
       !        end do

       block
         type(hash_t) :: hash_table
         integer(kind=int64) :: val,it
         logical :: found
         
         call hash_table%set_all_idx(bufid_tmp, ix_bufid)
         do it = 1, nlinksl
            call hash_table%get(linkid(it), val, found)
            if(found .eqv. .true.) then
               bufidflag(val) = bufid_tmp(val)
            end if
         end do
         call hash_table%clear()
       end block

#ifdef MPP_LAND
      call com_write1dInt8(bufidflag,ix_bufid,gbufid,npid)
#endif
      if(ix_bufid .gt. 0) then
          if(allocated(bufidflag)) deallocate(bufidflag)
          if(allocated(bufid_tmp)) deallocate(bufid_tmp)
      endif
      if(allocated(linkid)) deallocate(linkid)
      if(my_id .eq. io_id) then
          allocate(bufsize(npid))
          allocate(g1bufid(ndata))
          call get1d_int(trim(nlst(did)%UDMAP_FILE),"overlaps",bufsize)
          g1bufid = -999
          i = 1
          do k = 1, npid
               do j = 1, bufsize(k)
                 g1bufid(i) = gbufid(k)
                 i = i + 1 
               end do
          enddo
          if(allocated(bufsize))  deallocate(bufsize)
      endif


      if(my_id .eq. io_id) then 
           if(allocated(gbufid)) deallocate(gbufid)
      endif


      if(lnsizes(my_id+1) .gt. 0) allocate(bufid(lnsizes(my_id+1) ))
      call pack_decomp_int8(g1bufid, ndata, nprocs_map, lnsizes, istart,bufid)
      if(my_id .eq. io_id) then 
            if(allocated(g1bufid)) deallocate(g1bufid)
      endif


      if(my_id .eq. io_id) then
          allocate(gbufw(ndata))
          call get1d_real8(trim(nlst(did)%UDMAP_FILE),"regridweight",gbufw)
      endif
      if(lnsizes(my_id+1) .gt. 0) allocate(bufw(lnsizes(my_id+1) ))
      call pack_decomp_real8(gbufw, ndata, nprocs_map, lnsizes, istart,bufw)
      if(my_id .eq. io_id) then 
          if(allocated(gbufw))     deallocate(gbufw)
      endif


        if(my_id .eq. io_id) then
           if(allocated(nprocs_map)) deallocate (nprocs_map)
           if(allocated(istart)) deallocate (istart)
        endif
        lnsize = lnsizes(my_id + 1)
        if(allocated(lnsizes)) deallocate(lnsizes)
#else
       call hydro_stop("FATAL ERROR in UDMP : sequential not defined.")
#endif

    end subroutine readUDMP

    subroutine UDMP2LOCAL(npid,ix,jx,rtmask, ter_rt_flag)
        implicit none
        integer :: i,j,k, ngrids, ix,jx, starti,startj, endi,endj, ii,jj, npid, kk
        integer, intent(in), dimension(ix,jx) :: rtmask
        integer, dimension(lnsize) :: lndflag,gridflag , tmpgridflag
        integer :: ter_rt_flag, m, c


!   find ngrids is 0 so that we need to mapping from subsurface runoff.
#ifdef MPP_LAND
        if(left_id .ge. 0) then
           starti = local_startx_rt  + 1
        else
           starti = local_startx_rt 
        endif
        if(down_id .ge. 0) then
           startj = local_starty_rt  + 1
        else
           startj = local_starty_rt 
        endif
        if(right_id .ge. 0) then
           endi = local_startx_rt + ix -2
        else
           endi = local_startx_rt + ix -1
        endif
        if(up_id .ge. 0) then
           endj = local_starty_rt + jx -2
        else
           endj = local_starty_rt + jx -1
        endif
#else
        starti = 1
        startj = 1
        endi = ix
        endj = jx
#endif
        gridflag = 0
        lndflag = 0
      
#ifdef MPP_LAND
        k = 0
        do i = 1, lnsize
           if(bufid(i) .gt. 0) then
                if(bufi(i) .ge. starti .and. bufj(i) .ge. startj .and. &
                    bufi(i) .le. endi   .and. bufj(i) .le. endj) then
                    if(k .eq. 0) then
                       k = 1
                    else
                       if(bufid(i) .ne. bufid(i-1)) k = k + 1
                    endif
                    lndflag(k) = lndflag(k) + 1
                    if(ter_rt_flag .eq. 1) then
                        if(rtmask(bufi(i)-local_startx_rt+1,bufj(i)-local_starty_rt+1) .ge. 0) then
                             gridflag(k) = gridflag(k) + 1 
                        endif
                    endif
                 endif
           endif
        end do

! decide how many mapping land grids on current domain
!       tmpgridflag = gridflag
#ifdef MPP_LAND
!       call mpp_collect_1d_int_mem(npid,tmpgridflag)
#endif

! decide how many user defined links on current domain
        kk = k
        LNUMRSL = 0 
        do k = 1, lnsize
           if(lndflag(k) .gt. 0) LNUMRSL = LNUMRSL + 1
        enddo


        if(LNUMRSL .gt. 0) then 

               allocate(LUDRSL(LNUMRSL))
               allocate( basns_area(LNUMRSL) )
        else
! When MPI is performed,for every subdomain in each process, all the links 
! are listed and if there is no link in the subdomain then it is calling 
! cleanBuf (memory cleaning purposes), this used to print a warning 
! that is not necessary for the user to see it, therefore it is been commented out here
!               write(6,*) "Warning: no routing links found."
               call cleanBuf()
               return
        endif

        kk = 0
        do k = 1, lnsize 
           if( bufid(k) .ge. 0 ) then
             if (bufi(k) .ge. starti .and. bufj(k) .ge. startj .and. &
                    bufi(k) .le. endi   .and. bufj(k) .le. endj ) then
                 if(kk .eq. 0) then
                       kk = 1
                 else
                       if(bufid(k) .ne. bufid(k-1)) kk = kk + 1
                 endif
                 LUDRSL(kk)%myid = bufid(k) 
                 LUDRSL(kk)%ngrids = -999
                 if(gridflag(kk) .gt. 0) then
                   LUDRSL(kk)%ngrids = gridflag(kk)
                   if(.not. allocated(LUDRSL(kk)%weight) ) then
                         allocate( LUDRSL(kk)%weight(LUDRSL(kk)%ngrids ))
                         allocate( LUDRSL(kk)%grid_i(LUDRSL(kk)%ngrids ))
                         allocate( LUDRSL(kk)%grid_j(LUDRSL(kk)%ngrids ))
                         allocate( LUDRSL(kk)%nodeArea(LUDRSL(kk)%ngrids ))
                   endif
                 endif
!  define bucket variables
                 LUDRSL(kk)%ncell = lndflag(kk) 
                 if(.not. allocated(LUDRSL(kk)%cellweight) ) then
                     allocate( LUDRSL(kk)%cellweight(LUDRSL(kk)%ncell))
                     allocate( LUDRSL(kk)%cell_i(LUDRSL(kk)%ncell))
                     allocate( LUDRSL(kk)%cell_j(LUDRSL(kk)%ncell))
                     allocate( LUDRSL(kk)%cellArea(LUDRSL(kk)%ncell))
                 endif
             endif
           endif
        enddo


! maping grid_i, grid_j and weight
        kk = 0
        m  = 1
        c  = 1
        do i = 1, lnsize 
               if( (bufid(i) .ge. 0)  ) then 
                   if(bufi(i) .ge. starti .and. bufj(i) .ge. startj .and. &
                      bufi(i) .le. endi   .and. bufj(i) .le. endj) then
                      if(kk .eq. 0) then
                         kk = 1
                      else
                         if(bufid(i) .ne. bufid(i-1)) then 
                             kk = kk + 1
                             m  = 1
                             c  = 1
                         endif
                      endif

                      if(LUDRSL(kk)%ngrids .gt. 0) then 
                          if(rtmask(bufi(i)-local_startx_rt+1,bufj(i)-local_starty_rt+1) .ge. 0) then
                             LUDRSL(kk)%grid_i(m) = bufi(i) - local_startx_rt+1
                             LUDRSL(kk)%grid_j(m) = bufj(i) - local_starty_rt+1
                             LUDRSL(kk)%weight(m) = bufw(i) 
                             m  = m  + 1
                          endif
                      endif
!! begin define bucket variables
                          LUDRSL(kk)%cell_i(c) = bufi(i) - local_startx_rt+1
                          LUDRSL(kk)%cell_j(c) = bufj(i) - local_starty_rt+1
                          LUDRSL(kk)%cellWeight(c) = bufw(i)
                          c  = c  + 1
!! end define bucket variables 
                   endif
                endif
        end do

        call cleanBuf()

#else
        call hydro_stop("FATAL ERROR in UDMP: Sequential not work.")
#endif
   
    end subroutine UDMP2LOCAL

    subroutine cleanBuf()
        if(allocated(bufi))  deallocate(bufi)
        if(allocated(bufj))  deallocate(bufj)
        if(allocated(bufw))  deallocate(bufw)
        if(allocated(bufid))  deallocate(bufid)
    end subroutine cleanBuf

     subroutine get_dimension(fileName, ndata,npid)
            implicit none
            character(len=*) fileName
            integer ncid , iret, ndata,npid, dimid
#ifdef MPP_LAND
            if(my_id .eq. IO_id) then
#endif
               iret = nf_open(fileName, NF_NOWRITE, ncid)
               if (iret /= 0) then
                  write(*,'("FATAL ERROR: Problem opening mapping file: ''", A, "''")') &
                       trim(fileName)
                  call hydro_stop("In get_dimension() - Problem opening mapping file.")
               endif
        
               iret = nf_inq_dimid(ncid, "polyid", dimid)
        
               if (iret /= 0) then
                  print*, "nf_inq_dimid:  polyid"
                  call hydro_stop("In get_dimension() - nf_inq_dimid:  polyid")
               endif
           
               iret = nf_inq_dimlen(ncid, dimid, npid)
           
               iret = nf_inq_dimid(ncid, "data", dimid)
               if (iret /= 0) then
                          print*, "nf_inq_dimid:  data"
                          call hydro_stop("In get_file_dimension() - nf_inq_dimid:  data")
               endif
        
               iret = nf_inq_dimlen(ncid, dimid, ndata)
               iret = nf_close(ncid)
#ifdef MPP_LAND
            endif
            call mpp_land_bcast_int1(ndata)
            call mpp_land_bcast_int1(npid)
#endif
            return
     end subroutine get_dimension

       subroutine get1d_real8(fileName,var_name,out_buff)
          implicit none
          integer :: ivar, iret,varid,ncid
          real*8 out_buff(:)
          character(len=*), intent(in) :: var_name
          character(len=*), intent(in) :: fileName

          iret = nf_open(trim(fileName), NF_NOWRITE, ncid)
          if (iret .ne. 0) then
            print*,"failed to open the netcdf file: ",trim(fileName)
            call hydro_stop("In get1d_real8() - failed to open the netcdf file.")
            return
          endif
          ivar = nf_inq_varid(ncid,trim(var_name),  varid)
          if(ivar .ne. 0) then
               write(6,*) "Read Variable Error file: ",trim(fileName)
               write(6,*) "Read Error: could not find ",trim(var_name)
               call hydro_stop("In get1d_real8() - failed to read netcdf varialbe name. ")
          end if
          iret = nf_get_var_double(ncid, varid, out_buff)
          iret = nf_close(ncid)
      end subroutine get1d_real8

       subroutine get1d_int(fileName,var_name,out_buff)
          implicit none
          integer :: ivar, iret,varid,ncid
          integer out_buff(:)
          character(len=*), intent(in) :: var_name
          character(len=*), intent(in) :: fileName

          iret = nf_open(trim(fileName), NF_NOWRITE, ncid)
          if (iret .ne. 0) then
            print*,"FATAL ERROR: Failed to open the netcdf file: ",trim(fileName)
            call hydro_stop("In get1d_int() -  Failed to open the netcdf file")
            return
          endif
          ivar = nf_inq_varid(ncid,trim(var_name),  varid)
          if(ivar .ne. 0) then
               write(6,*) "Read Variable Error file: ",trim(fileName)
               write(6,*) "Read Error: could not find ",trim(var_name)
               call hydro_stop("In get1d_int() - failed to read netcdf variable name.")
          end if
          iret = nf_get_var_int(ncid, varid, out_buff)
          iret = nf_close(ncid)
      end subroutine get1d_int

    subroutine get1d_int64(fileName,var_name,out_buff)
        implicit none
        integer :: ivar, iret,varid,ncid
        integer(kind=int64) out_buff(:)
        character(len=*), intent(in) :: var_name
        character(len=*), intent(in) :: fileName

        iret = nf_open(trim(fileName), NF_NOWRITE, ncid)
        if (iret .ne. 0) then
            print*,"FATAL ERROR: Failed to open the netcdf file: ",trim(fileName)
            call hydro_stop("In get1d_int() -  Failed to open the netcdf file")
            return
        endif
        ivar = nf_inq_varid(ncid,trim(var_name),  varid)
        if(ivar .ne. 0) then
            write(6,*) "Read Variable Error file: ",trim(fileName)
            write(6,*) "Read Error: could not find ",trim(var_name)
            call hydro_stop("In get1d_int() - failed to read netcdf variable name.")
        end if
        iret = nf_get_var_int64(ncid, varid, out_buff)
        iret = nf_close(ncid)
    end subroutine get1d_int64

      subroutine getUDMP_area(cell_area)
         implicit none
         integer i,j,k, m
         real, dimension(:,:) :: cell_area
         do k  = 1, LNUMRSL
            if(LUDRSL(k)%ngrids .gt. 0) then
                do m = 1, LUDRSL(k)%ngrids
                    LUDRSL(k)%nodeArea(m) = cell_area(LUDRSL(k)%grid_i(m),LUDRSL(k)%grid_j(m)) 
                enddo
            endif
                do m = 1, LUDRSL(k)%ncell
                    LUDRSL(k)%cellArea(m) = cell_area(LUDRSL(k)%cell_i(m),LUDRSL(k)%cell_j(m)) 
                enddo
           
            basns_area(k) = 0 
            do m = 1, LUDRSL(k)%ncell
                    basns_area(k) = basns_area(k) + &
                          cell_area(LUDRSL(k)%cell_i(m),LUDRSL(k)%cell_j(m)) * LUDRSL(k)%cellWeight(m) 
            enddo
            
         end do
      end subroutine getUDMP_area
    
      subroutine get_basn_area_nhd(inOut)
         implicit none
         real, dimension(:) :: inOut
         real, dimension(gnpid) :: buf
#ifdef MPP_LAND
         call updateLinkV(basns_area, inOut)
#else
         inOut = basns_area
#endif

      
      end subroutine get_basn_area_nhd

      subroutine get_nprocs_map(ix,jx,bufi,bufj,nprocs_map,ndata)
          implicit none
          integer,dimension(:)  :: bufi, bufj,nprocs_map 
!          integer, allocatable, dimension(:) ::  lbufi,lbufj, lmap
          integer  :: ndata, lsize, ix,jx
          integer, dimension(ix,jx) :: mask
          integer, allocatable,dimension(:,:) :: gmask

        integer :: i,j,k, starti,startj, endi,endj, ii,jj, npid, kk
#ifdef MPP_LAND
           
          mask = my_id + 1
          if(my_id .eq. IO_id) allocate(gmask(global_rt_nx, global_rt_ny))

          call MPP_LAND_COM_INTEGER(mask,IX,JX,99)
          call write_IO_rt_int(mask, gmask) 

          if(my_id .eq. io_id ) then
             nprocs_map = -999
             do i = 1, ndata
                  if( (bufi(i) .gt. 0 .and. bufi(i) .le. global_rt_nx) .and.  &
                     (bufj(i) .gt. 0 .and. bufj(i) .le. global_rt_ny) ) then
                     nprocs_map(i) = gmask(bufi(i), bufj(i))
                     if( gmask(bufi(i), bufj(i)) .lt. 0) then
                         write(6,*) "mapping error in gmask : ", bufi(i) ,bufj(i)
                     endif
                  else
                      write(6,*) "no mapping for i,j : ", bufi(i) ,bufj(i)
                  endif
             end do

             if(allocated(gmask)) deallocate(gmask)
          endif
#else
        call hydro_stop("FATAL ERROR in UDMP: Sequential not work.")
#endif


      end subroutine get_nprocs_map


end module module_UDMAP