routing__D8.F90

Go to the documentation of this file.

module routing__d8
 
  use iso_c_binding, only : c_short, c_int, c_float, c_double, c_bool
  use constants_and_conversions
  use data_catalog
  use data_catalog_entry
  use exceptions
  use grid
  use fstring
  implicit none
 
  private
 
  public :: routing_d8_initialize, d8_undetermined,                                 &
            sum_of_upslope_cells, number_of_upslope_connections,                    &
            get_cell_index, get_target_index, get_sort_order
 
  type (data_catalog_entry_t), pointer :: pd8_flowdir    ! data catalog object => D8 flow direction grid
 
  integer (c_int), allocatable     :: target_row(:,:)
  integer (c_int), allocatable     :: target_col(:,:)
  logical (c_bool), allocatable    :: is_downslope_target_marked(:,:)
  integer (c_int), allocatable     :: sum_of_upslope_cells(:,:)
  integer (c_int), allocatable     :: number_of_upslope_connections(:,:)
 
  !> @TODO remove redundant data elements; row1d, col1d, etc. are now also stored in the model data structure.
 
  integer (c_int), allocatable     :: row2d(:,:)
  integer (c_int), allocatable     :: col2d(:,:)
  integer (c_int), allocatable     :: row1d(:)
  integer (c_int), allocatable     :: col1d(:)
 
  integer (c_int), allocatable     :: row_index(:)
  integer (c_int), allocatable     :: column_index(:)
  integer (c_int), allocatable     :: sort_order_l(:)
  integer (c_int), allocatable     :: target_index_l(:)
 
  integer (c_int), parameter       :: d8_east         = 1
  integer (c_int), parameter       :: d8_southeast    = 2
  integer (c_int), parameter       :: d8_south        = 4
  integer (c_int), parameter       :: d8_southwest    = 8
  integer (c_int), parameter       :: d8_west         = 16
  integer (c_int), parameter       :: d8_northwest    = 32
  integer (c_int), parameter       :: d8_north        = 64
  integer (c_int), parameter       :: d8_northeast    = 128
 
  integer (c_int), parameter       :: d8_undetermined = -999
 
  ! idea here is to create an array of row and column numbers, then use the
  ! PACK statement to create a vector of row-column numbers in the same order that
  ! the rest of the vectors are packed.
 
  ! the vector of row-column numbers may then be used to make the required D8 routing connections.
  ! there should not be anumber_of_rows difference between ROW_INDEX, COL_INDEX, and SORT_ORDER_l 
 
contains
 
!
! concept:
!          ^
! +-----+--|--+-----+
! |  -> |     |  <- |    with cell indices of  1   2   3      d8 flow dirs of    1  64  16
! +-----+--^--+-----+                              4                                64
!       |  |  |
!       +-----+
!
! solution must be ordered upslope to downslope: cell 1, 3, 4, then 2 as one possible sort order
!
! iteration_index            cell_index       target_index
! -----------------         ------------     --------------
!         1                       1                2
!         2                       3                2
!         3                       4                2
!         4                       2              -9999
!
! function 'get_cell_index' simply looks up the appropriate cell index given a predefined
! cell sort order.
!
  elemental function get_cell_index( iteration_index )  result( cell_index )
 
    integer (c_int), intent(in)    :: iteration_index
    integer (c_int)                :: cell_index
 
    ! if D8 flow routing has not been initialized, report back the iteration index
    cell_index = iteration_index
 
    if ( allocated( sort_order_l ) )   cell_index = sort_order_l( iteration_index )
 
 
  end function get_cell_index
 
!---------------------------------------------------------------------------------------------------
 
elemental function get_target_index( iteration_index )  result( target_index )
 
  integer (c_int), intent(in)    :: iteration_index
  integer (c_int)                :: target_index
 
   if ( allocated( target_index_l ) ) then
 
     target_index = target_index_l( iteration_index )
 
  else
 
    target_index = -9999
 
  endif
 
end function get_target_index
 
!---------------------------------------------------------------------------------------------------
 
  elemental function get_sort_order( cell_index )  result( iteration_index )
 
    integer (c_int), intent(in)    :: cell_index
    integer (c_int)                :: iteration_index
 
    ! [ LOCALS ]
    integer (c_int)  :: current_cell_index
    logical (c_bool) :: found_match
 
    found_match = false
 
    if ( allocated( sort_order_l ) ) then
 
      do iteration_index=1, ubound( sort_order_l, 1 )
 
         current_cell_index = sort_order_l( iteration_index )
         if ( current_cell_index == cell_index ) then
           found_match = true
          exit
        endif
 
      enddo
 
    endif
 
    ! if this function is called and routing is turned off, the 'sort_order'
    ! will be the same as the cell_index
    if ( .not. found_match )  iteration_index = cell_index
 
  end function get_sort_order
 
!---------------------------------------------------------------------------------------------------
 
  subroutine routing_d8_initialize( lActive, sort_order )
 
    logical (c_bool), intent(in)    :: lactive(:,:)
    integer (c_int), intent(inout)  :: sort_order(:)
 
    ! [ LOCALS ]
    integer (c_int)                 :: number_of_cols
    integer (c_int)                 :: number_of_rows
    integer (c_int)                 :: istat
    integer (c_int)                 :: column_num
    integer (c_int)                 :: row_num
    integer (c_int)                 :: iteration_index
    integer (c_int)                 :: icount
    character (len=256)                  :: sbuf
 
    integer (c_int) :: col_lbound, col_ubound
    integer (c_int) :: row_lbound, row_ubound
    integer (c_int) :: cell_index, target_index
    integer (c_int) :: iunitnum
    integer (c_int) :: target_row_num, target_col_num
 
    type (general_grid_t), pointer  :: ptempgrid
 
    col_lbound = lbound(lactive, 1)
    col_ubound = ubound(lactive, 1)
 
    row_lbound = lbound(lactive, 2)
    row_ubound = ubound(lactive, 2)
 
    icount = 0
 
    number_of_cols = ubound(lactive, 1)
    number_of_rows = ubound(lactive, 2)
 
    allocate( target_row( number_of_cols, number_of_rows ), stat=istat )
    call assert( istat==0, "Problem allocating memory", __file__, __line__ )
 
    allocate( target_col( number_of_cols, number_of_rows ), stat=istat )
    call assert( istat==0, "Problem allocating memory", __file__, __line__ )
 
 
    allocate( is_downslope_target_marked( number_of_cols, number_of_rows ), stat=istat )
    call assert( istat==0, "Problem allocating memory", __file__, __line__ )
 
    allocate( col2d( number_of_cols, number_of_rows ), stat=istat )
    call assert( istat==0, "Problem allocating memory", __file__, __line__ )
 
    allocate( row2d( number_of_cols, number_of_rows ), stat=istat )
    call assert( istat==0, "Problem allocating memory", __file__, __line__ )
 
    allocate( col1d( count( lactive) ), stat=istat )
    call assert( istat==0, "Problem allocating memory", __file__, __line__ )
 
    allocate( row1d( count( lactive) ), stat=istat )
    call assert( istat==0, "Problem allocating memory", __file__, __line__ )
 
    allocate( sum_of_upslope_cells( number_of_cols, number_of_rows ), stat=istat )
    call assert( istat==0, "Problem allocating memory", __file__, __line__ )
 
    allocate( number_of_upslope_connections( number_of_cols, number_of_rows ), stat=istat )
    call assert( istat==0, "Problem allocating memory", __file__, __line__ )
 
    allocate( row_index( count( lactive) ), stat=istat )
    call assert( istat==0, "Problem allocating memory", __file__, __line__ )
 
    allocate( column_index( count( lactive) ), stat=istat )
    call assert( istat==0, "Problem allocating memory", __file__, __line__ )
 
    allocate( target_index_l( count( lactive) ), stat=istat )
    call assert( istat==0, "Problem allocating memory", __file__, __line__ )
 
    allocate( sort_order_l( count( lactive) ), stat=istat )
    call assert( istat==0, "Problem allocating memory", __file__, __line__ )
 
    ! locate the data structure associated with the gridded flow direction entries
    pd8_flowdir => dat%find("FLOW_DIRECTION")
    if ( .not. associated(pd8_flowdir) ) &
        call die("A FLOW_DIRECTION grid must be supplied in order to make use of this option.", __file__, __line__)
 
    call pd8_flowdir%getvalues()
 
    ptempgrid=>grid_create( bnds%iNumCols, bnds%iNumRows, bnds%fX_ll, bnds%fY_ll, &
      bnds%fX_ur, bnds%fY_ur, datatype_int )
 
    ptempgrid%iData = pd8_flowdir%pGrdBase%iData
 
    call grid_writearcgrid("D8_Flow_Direction_Grid.asc", ptempgrid)
 
    call routing_d8_assign_downstream_row_col( lactive )
 
    do row_num=row_lbound, row_ubound
      do column_num=col_lbound, col_ubound
 
        col2d( column_num, row_num ) = column_num   ! right-most index should vary slowest
        row2d( column_num, row_num ) = row_num
 
      enddo
    enddo
 
    ! create 1D vectors of column and row numbers for full grid
    col1d = pack( col2d, lactive )
    row1d = pack( row2d, lactive )
 
    call routing_d8_determine_solution_order( lactive )
 
    ptempgrid%iData = target_row
    call grid_writearcgrid("D8_TARGET_ROW_Grid.asc", ptempgrid)
 
    ptempgrid%iData = target_col
    call grid_writearcgrid("D8_TARGET_COL_Grid.asc", ptempgrid)
 
    open( newunit=iunitnum, file=trim(output_directory_name)//trim("D8_routing_table.txt"),   &
      iostat=istat, status="REPLACE")
 
    write(iunitnum,*) "Sort_Order"//tab//"CELL_INDEX"//tab//"TARGET_INDEX"//tab//"From_COL"    &
      // tab//"From_ROW"//tab//"To_COL"//tab//"To_ROW"//tab//"D8_flowdir"//tab                     &
      //"Num_Adjacent_Upslope_Connections"//tab//"Sum_of_Upslope_Contributing_Cells"
 
     ! solution order has been determined; remaining code simply writes a summary to a
     ! file for further analysis
    do iteration_index = 1, ubound(col1d,1)
 
      sbuf = ""
 
      cell_index = get_cell_index( iteration_index )
      target_index = get_target_index( iteration_index )
 
      sort_order(iteration_index) = cell_index
 
      associate(                                                               &
        rownum => row1d( cell_index ),                                         &
        colnum => col1d( cell_index ),                                         &
!        Target_row => ROW1D( target_index ),                                   &
!        Target_col => COL1D( target_index ),                                   &
        d8_flowdir => pd8_flowdir%pGrdBase%iData( col1d( cell_index ),         &
                        row1d( cell_index ) ),                                 &
        num_adjacent => number_of_upslope_connections( col1d( cell_index ),    &
                        row1d( cell_index ) ),                                 &
        sum_upslope => sum_of_upslope_cells( col1d( cell_index ),              &
                        row1d( cell_index ) )  )
 
        write(sbuf,*)  iteration_index, tab, cell_index, tab, target_index
        write(sbuf,*) trim(sbuf)//tab//ascharacter( colnum )//tab              &
          //ascharacter( rownum )
 
        if ( target_index > 0 .and. target_index <= ubound(col1d,1) ) then
          target_row_num = row1d( target_index )
          target_col_num = col1d( target_index )
        else
          target_row_num = -9999
          target_col_num = -9999
        endif
 
        write(sbuf,*) trim(sbuf)//tab//ascharacter( target_col_num )//tab    &
          //ascharacter( target_row_num )//tab//ascharacter( d8_flowdir )//tab &
          //ascharacter( num_adjacent )//tab//ascharacter( sum_upslope )
 
        write(iunitnum,*)  trim(sbuf)
 
      end associate
 
    enddo
 
    close ( iunitnum )
 
  end subroutine routing_d8_initialize
 
!--------------------------------------------------------------------------------------------------
 
  function routing_d8_get_index( iCol, iRow )   result( cell_index )
 
    integer (c_int), intent(in)   :: icol
    integer (c_int), intent(in)   :: irow
    integer (c_int)               :: cell_index
 
    ! [ LOCALS ]
    integer (c_int)   :: iindex
    logical (c_bool)  :: lfound
 
    iindex        = -9999
    cell_index    = -9999
    lfound        = false
 
    ! iterate over 1D vector of column numbers
    do iindex = lbound(col1d,1), ubound(col1d,1)
 
      if( col1d( iindex ) == icol  .and.  row1d( iindex ) == irow ) then
        cell_index = iindex
        lfound = true
        exit
      endif
 
    enddo
 
!    if ( .not. lFound ) call warn("Did not find matching column and row number; col=" &
!      //asCharacter(iCol)//"; row="//asCharacter(iRow), __FILE__, __LINE__ )
 
  end function routing_d8_get_index
 
!--------------------------------------------------------------------------------------------------
 
  subroutine routing_d8_assign_downstream_row_col( lActive )
 
    logical (c_bool), intent(in)    :: lActive(:,:)
 
    ! [ LOCALS ]
    integer (c_int) :: row_num
    integer (c_int) :: column_num
    integer (c_int) :: col_lbound, col_ubound
    integer (c_int) :: row_lbound, row_ubound
 
    col_lbound = lbound(lactive, 1)
    col_ubound = ubound(lactive, 1)
 
    row_lbound = lbound(lactive, 2)
    row_ubound = ubound(lactive, 2)
 
    associate( dir => pd8_flowdir%pGrdBase%iData )
 
        do row_num=row_lbound, row_ubound
        do column_num=col_lbound, col_ubound
 
          select case ( dir( column_num, row_num ) )
 
            case ( d8_east )
 
              target_col(column_num, row_num) = column_num + 1
              target_row(column_num, row_num) = row_num
 
            case ( d8_southeast )
 
              target_col(column_num, row_num) = column_num + 1
              target_row(column_num, row_num) = row_num + 1
 
            case ( d8_south )
 
              target_col(column_num, row_num) = column_num
              target_row(column_num, row_num) = row_num + 1
 
            case ( d8_southwest )
 
              target_col(column_num, row_num) = column_num - 1
              target_row(column_num, row_num) = row_num + 1
 
            case ( d8_west )
 
              target_col(column_num, row_num) = column_num - 1
              target_row(column_num, row_num) = row_num
 
            case ( d8_northwest )
 
              target_col(column_num, row_num) = column_num - 1
              target_row(column_num, row_num) = row_num - 1
 
            case ( d8_north )
 
              target_col(column_num, row_num) = column_num
              target_row(column_num, row_num) = row_num - 1
 
            case ( d8_northeast )
 
              target_col(column_num, row_num) = column_num + 1
              target_row(column_num, row_num) = row_num - 1
 
            case default
 
              target_col(column_num, row_num) = d8_undetermined
              target_row(column_num, row_num) = d8_undetermined
 
          end select
 
        enddo
 
      enddo
 
    end associate
 
 
  end subroutine routing_d8_assign_downstream_row_col
 
!--------------------------------------------------------------------------------------------------
 
  subroutine routing_d8_determine_solution_order( lActive )
 
    use grid
 
    logical (c_bool), intent(in)    :: lActive(:,:)
 
    ! [ LOCALS ]
    integer (c_int)  :: row_num
    integer (c_int)  :: column_num
    integer (c_int)  :: iColsrch
    integer (c_int)  :: iRowsrch
    integer (c_int)  :: iNumberOfChangedCells
    integer (c_int)  :: col_lbound, col_ubound
    integer (c_int)  :: row_lbound, row_ubound
    integer (c_int)  :: iUpslopeSum, iUpslopeConnections
    logical (c_bool) :: are_there_unmarked_upslope_cells
    logical (c_bool) :: lCircular
    integer (c_int)  :: iNumberRemaining
    integer (c_int)  :: indx, k, iCount
    integer (c_int)  :: num_cells_marked_this_iteration
    integer (c_int)  :: iPasses
    integer (c_int)  :: iPassesWithoutChange
    type (GENERAL_GRID_T), pointer  :: pTempGrid
 
    col_lbound = lbound(lactive, 1)
    col_ubound = ubound(lactive, 1)
 
    row_lbound = lbound(lactive, 2)
    row_ubound = ubound(lactive, 2)
 
    is_downslope_target_marked = false
    sum_of_upslope_cells = 0_c_int
 
    indx = 0
    ipasses = 0
    ipasseswithoutchange = 0
 
    ptempgrid=>grid_create( bnds%iNumCols, bnds%iNumRows, bnds%fX_ll, bnds%fY_ll, &
      bnds%fX_ur, bnds%fY_ur, datatype_int )
 
 
main_loop: do
 
      inumberofchangedcells = 0_c_int
      num_cells_marked_this_iteration = 0_c_int
      ipasses = ipasses + 1
 
      ! iterate over entire model domain
        do row_num=row_lbound, row_ubound
        do column_num=col_lbound, col_ubound
 
            if ( .not. lactive(column_num, row_num) ) cycle
 
          ! we have already determined the identiy of the downslope cell; skip this one
            if ( is_downslope_target_marked(column_num, row_num) ) cycle
 
          inumberofchangedcells = inumberofchangedcells + 1
            iupslopesum = 0_c_int
            iupslopeconnections = 0_c_int
            are_there_unmarked_upslope_cells = false
          lcircular = false
 
          ! search the 8 cells immediately adjacent to the current cell
    local_search: do irowsrch=max( row_num-1, row_lbound),min( row_num+1, row_ubound)
                  do icolsrch=max( column_num-1, col_lbound),min( column_num+1, col_ubound)
 
                            ! if adjacent cell points to current cell, determine what to do with runoff
                            if ( ( target_col(icolsrch, irowsrch) == column_num ) &
                                .and. ( target_row(icolsrch, irowsrch) == row_num ) ) then
 
                      ! if adjacent cell falls outside the area of active cells,
                      ! ignore it and move on
                      if ( .not. lactive( icolsrch, irowsrch ) )  cycle
 
                      ! if the target of the current cell points back at the adjacent
                      ! cell, mark current cell as having a "circular" connection
                      if ( ( target_col( column_num, row_num ) == icolsrch )                &
                        .and. ( target_row( column_num, row_num ) == irowsrch ) )  lcircular = true
 
                      ! if the adjacent cell is marked (that is, has no unresolved
                      ! upslope contributions), add the 1 to the number of connections
                      ! and add the adjacent cells' sum of upslope cells to the
                      ! current cells' running sum of upslope cells
                              if(   is_downslope_target_marked( icolsrch, irowsrch ) ) then
 
                        iupslopesum = iupslopesum + sum_of_upslope_cells(icolsrch, irowsrch)
                        iupslopeconnections = iupslopeconnections + 1
 
                      ! add number of upslope cells and connections to temporary variables
                      ! if we get to the end of the search and find no unmarked adjacent cells,
                      ! we have determined the number of upslope contributing cells and will
                      ! be able to mark the current cell
                              else
 
                        are_there_unmarked_upslope_cells = true
 
                              endif
 
                            endif
 
                    enddo
                  enddo local_search
 
                  ! OK this is the end of the local search area; did we uncover anumber_of_rows unmarked cells
                  ! contributing flow to the current cell? if so, ignore and move on. otherwise,
                  ! mark current cell as marked, update stats, and continue with next cell
 
                  if ( .not. are_there_unmarked_upslope_cells ) then
!                  .or. (iUpslopeConnections == 1 .and. lCircular ) ) then
 
                  num_cells_marked_this_iteration = num_cells_marked_this_iteration + 1
                  indx = indx + 1
 
                  ! need to add one to the sum already tabulated to account for the
                  ! current cell as well
                    sum_of_upslope_cells( column_num, row_num ) = iupslopesum + 1
                    number_of_upslope_connections( column_num, row_num ) = iupslopeconnections
                    is_downslope_target_marked( column_num, row_num ) = true
                  column_index( indx )  = column_num
                  row_index( indx )     = row_num
                  sort_order_l( indx )    = routing_d8_get_index( column_num, row_num )
                  target_index_l( indx )  = routing_d8_get_index( target_col( column_num, row_num ), &
                                                                                   target_row( column_num, row_num ) )
 
                  ! print *, "   ===> assigning order number: "
                  ! print *, "        indx                =",indx
                  ! print *, "        sort_order                =", SORT_ORDER_l( indx )
                  ! print *, "        target_index              =", TARGET_INDEX_l( indx )
                  ! print *, "        target_index(SORT_ORDER_l =", TARGET_INDEX_l( SORT_ORDER_l( indx ) )
 
                  if ( lcircular )  target_index_l( sort_order_l( indx ) ) = d8_undetermined
 
                ! if cells are not being marked, record cell as D8_UNDETERMINED
                ! and move on
                elseif ( ipasseswithoutchange > 10 ) then
 
                  num_cells_marked_this_iteration = num_cells_marked_this_iteration + 1
                  indx = indx + 1
 
                  sum_of_upslope_cells( column_num, row_num ) = iupslopesum
                  number_of_upslope_connections( column_num, row_num ) = iupslopeconnections
                  is_downslope_target_marked( column_num, row_num ) = true
                  column_index( indx ) = column_num
                  row_index( indx ) = row_num
                  sort_order_l( indx ) = routing_d8_get_index( column_num, row_num )
 
                  target_index_l( sort_order_l( indx ) ) = d8_undetermined
 
                endif
 
          enddo
        enddo
 
!      if (iNumberOfChangedCells == 0) exit main_loop
 
      inumberremaining = count( lactive ) - count( is_downslope_target_marked )
 
      if ( inumberremaining==0 )   exit main_loop
 
      print *, "### determining solution order... ", count( is_downslope_target_marked ), " cells marked so far " &
        //"out of ", count( lactive ), " active cells."
 
 
      if ( num_cells_marked_this_iteration==0 ) then
 
        ipasseswithoutchange = ipasseswithoutchange + 1
 
        where( is_downslope_target_marked .or.  ( .not. lactive ) )
          ptempgrid%iData = -1
        elsewhere
          ptempgrid%iData = pd8_flowdir%pGrdBase%iData
        end where
 
        write(*,"(/,1x,'Summary of remaining unmarked cells')")
 
        ! loop over possible (legal) values of the flow direction grid
        do k=0,128
          icount=count( .not. is_downslope_target_marked                         &
            .and. pd8_flowdir%pGrdBase%iData==k .and. lactive )
          if( icount > 0 ) then
            write(*,fmt="(3x,i8,' unmarked grid cells have flowdir value: ',i8)") icount, k
          end if
        end do
 
        write(*,fmt="(3x,a)") repeat("-",60)
        write(*,fmt="(3x,i8,' Total cells with nonzero flow " &
          //"direction values')") count( pd8_flowdir%pGrdBase%iData > 0 )
 
 
        call grid_writearcgrid("iteration"//ascharacter(ipasses)// &
          "problem_gridcells.asc", ptempgrid)
 
      end if
 
 
    enddo main_loop
 
    ptempgrid%iData = sum_of_upslope_cells
    call grid_writearcgrid("D8_Flow_Routing__Upslope_Contributing_Area.asc", ptempgrid)
 
    ptempgrid%iData = number_of_upslope_connections
    call grid_writearcgrid("D8_Flow_Routing__Number_of_Adjacent_Upslope_Connections.asc", ptempgrid)
 
    call grid_destroy( ptempgrid )
 
  end subroutine routing_d8_determine_solution_order
 
  !------------------------------------------------------------------------------------------------
 
  ! subroutine calculate_routing_D8( indx )
  !
  !     if ( (    TARGET_INDEX_l( indx ) >= lbound( SORT_ORDER_l, 1) ) &
  !       .and. ( TARGET_INDEX_l( indx ) <= ubound( SORT_ORDER_l, 1) ) ) then
  !
  !       this%runon( TARGET_INDEX_l( index ) ) = this%runoff( SORT_ORDER_l( index ) )
  !
  !     endif
  !
  !   enddo
  !
  ! end subroutine calculate_routing_D8
 
end module routing__d8