Line data    Source code

       1             : !> \file    mo_mrm_river_head.f90
       2             : !> \brief   \copybrief mo_mrm_river_head
       3             : !> \details \copydetails mo_mrm_river_head
       4             : 
       5             : !> \brief   River head calculation
       6             : !> \details Enables river - groundwater interaction in mRM.
       7             : !> \authors Lennart Schueler
       8             : !> \date    Jul 2018
       9             : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
      10             : !! mHM is released under the LGPLv3+ license \license_note
      11             : !> \ingroup f_mrm
      12             : module mo_mrm_river_head
      13             :   use mo_common_variables,     only : level0, domainMeta
      14             :   use mo_mrm_global_variables, only : L0_L11_remap, L11_bankfull_runoff_in, &
      15             :                                       L0_slope, L0_channel_elevation
      16             :   use mo_kind,                 only : i4, dp
      17             :   use mo_common_constants,     only : nodata_dp
      18             :   use mo_append,               only : append
      19             : 
      20             :   implicit none
      21             : 
      22             :   private
      23             : 
      24             :   public ::  init_masked_zeros_l0
      25             :   public :: calc_channel_elevation
      26             :   public :: calc_river_head
      27             : 
      28             :   contains
      29             : 
      30             :   !> \brief allocates memory for L0 variable
      31           0 :   subroutine init_masked_zeros_l0(iDomain, data)
      32             :     integer(i4), intent(in) :: iDomain
      33             :     real(dp), dimension(:), allocatable, intent(inout) :: data
      34           0 :     real(dp), dimension(:), allocatable :: dummy_1D
      35             : 
      36           0 :     allocate(dummy_1D(level0(iDomain)%nCells))
      37           0 :     dummy_1D = nodata_dp
      38           0 :     call append(data, dummy_1D)
      39           0 :     deallocate(dummy_1D)
      40           0 :   end subroutine init_masked_zeros_l0
      41             : 
      42             :   !> \brief calculates the channel elevation from the bankfull river discharge
      43           0 :   subroutine calc_channel_elevation()
      44           0 :     use mo_common_constants, only : nodata_i4
      45             :     use mo_common_variables, only : domainMeta, L0_elev
      46             :     use mo_mrm_global_variables, only : L0_fDir, L0_fAcc, L0_channel_depth
      47           0 :     real(dp), dimension(:,:), allocatable :: channel_dpth
      48           0 :     real(dp), dimension(:,:), allocatable :: channel_elev
      49           0 :     real(dp), dimension(:,:), allocatable :: slope
      50           0 :     real(dp), dimension(:,:), allocatable :: elev0
      51           0 :     integer(i4), dimension(:,:), allocatable :: fDir0
      52           0 :     integer(i4), dimension(:,:), allocatable :: fAcc0
      53             :     real(dp) n ! Manning's roughness coefficient
      54             :     integer(i4) :: nrows0, ncols0
      55             :     integer(i4) :: s0, e0
      56             :     integer(i4) i, j, k
      57             :     integer(i4) iDomain
      58             : 
      59           0 :     n = .045_dp ! m^-1/3 s from Sutanudjaja et al. 2011
      60             : 
      61           0 :     do iDomain = 1, domainMeta%nDomains
      62           0 :       nrows0 = level0(iDomain)%nrows
      63           0 :       ncols0 = level0(iDomain)%ncols
      64           0 :       s0 = level0(iDomain)%iStart
      65           0 :       e0 = level0(iDomain)%iEnd
      66           0 :       allocate(channel_dpth(nrows0, ncols0))
      67           0 :       allocate(channel_elev(nrows0, ncols0))
      68           0 :       allocate(elev0(nrows0, ncols0))
      69           0 :       allocate(fDir0(nrows0, ncols0))
      70           0 :       allocate(fAcc0(nrows0, ncols0))
      71           0 :       allocate(slope(nrows0, ncols0))
      72           0 :       channel_dpth(:,:) = nodata_dp
      73           0 :       channel_elev(:,:) = nodata_dp
      74           0 :       slope(:,:) = nodata_dp
      75             : 
      76           0 :       elev0(:,:) = unpack(L0_elev(s0:e0), level0(iDomain)%mask, nodata_dp)
      77           0 :       fDir0(:,:) = unpack(L0_fDir(s0:e0), level0(iDomain)%mask, nodata_i4)
      78           0 :       fAcc0(:,:) = unpack(L0_fAcc(s0:e0), level0(iDomain)%mask, nodata_i4)
      79             : 
      80           0 :       do k = 1, level0(iDomain)%nCells
      81           0 :         i = level0(iDomain)%CellCoor(k, 1)
      82           0 :         j = level0(iDomain)%CellCoor(k, 2)
      83           0 :         if (fAcc0(i,j) > 1) then
      84           0 :           slope(i,j) = calc_slope(iDomain, elev0, fDir0, i, j)
      85           0 :           channel_dpth(i,j) = ((n * &
      86           0 :   L11_bankfull_runoff_in(L0_L11_remap(iDomain)%lowres_id_on_highres(i,j))) / &
      87           0 :   (4.8 * slope(i,j)**.5))**.6
      88           0 :           channel_elev(i,j) = elev0(i,j) - channel_dpth(i,j)
      89             :         end if
      90             :       end do
      91             : 
      92             :       call append(L0_channel_depth, pack(channel_dpth(:,:), &
      93           0 :                   level0(iDomain)%mask))
      94             :       call append(L0_channel_elevation, pack(channel_elev(:,:), &
      95           0 :                   level0(iDomain)%mask))
      96             :       call append(L0_slope, pack(slope(:,:), &
      97           0 :                   level0(iDomain)%mask))
      98             : 
      99           0 :       deallocate(channel_dpth)
     100           0 :       deallocate(channel_elev)
     101           0 :       deallocate(elev0)
     102           0 :       deallocate(fDir0)
     103           0 :       deallocate(fAcc0)
     104           0 :       deallocate(slope)
     105             :     end do
     106           0 :   end subroutine calc_channel_elevation
     107             : 
     108             :   !> \brief calculates the river head
     109           0 :   subroutine calc_river_head(iDomain, L11_Qmod, river_head)
     110             :     integer(i4), intent(in) :: iDomain
     111             :     real(dp), dimension(:), intent(in) :: L11_Qmod
     112             :     real(dp), dimension(:), allocatable, intent(inout) :: river_head
     113             :     real(dp) :: n ! Manning's roughness coefficient
     114             :     integer(i4) :: s0, e0
     115             :     integer(i4) i, j, k, L11_ind
     116             : 
     117           0 :     n = .045_dp ! m^-1/3 s from Sutanudjaja et al. 2011
     118             : 
     119           0 :     s0 = level0(iDomain)%iStart
     120           0 :     e0 = level0(iDomain)%iEnd
     121             : 
     122           0 :     do k = s0, e0
     123           0 :       i = level0(iDomain)%CellCoor(k - s0 + 1, 1)
     124           0 :       j = level0(iDomain)%CellCoor(k - s0 + 1, 2)
     125           0 :       if (i >= 0 .and. i < 99999 .and. j  >= 0 .and. j < 99999) then
     126           0 :         L11_ind = L0_L11_remap(iDomain)%lowres_id_on_highres(i,j)
     127             :         ! TODO L11_Qmid(L11_ind) causes IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
     128           0 :         river_head(k) = L0_channel_elevation(k) + &
     129           0 :             (n * L11_Qmod(L11_ind) / L11_bankfull_runoff_in(L11_ind) / &
     130           0 :             L0_slope(k)**.5)**.6
     131             :       end if
     132             :     end do
     133             : 
     134           0 :   end subroutine calc_river_head
     135             : 
     136             :   !> \brief calculates domain slope
     137           0 :   function calc_slope(iDomain, elev0, fDir0, i, j) result(slope)
     138           0 :     use mo_common_variables, only: iFlag_cordinate_sys
     139             :     use mo_mrm_net_startup,      only: cellLength, moveDownOneCell
     140             :     integer(i4), intent(in) :: iDomain
     141             :     integer(i4), intent(in) :: i, j
     142             :     integer(i4), intent(in), dimension(:,:), allocatable :: fDir0
     143             :     real(dp), intent(in), dimension(:,:), allocatable :: elev0
     144           0 :     real(dp) :: slope, length
     145             :     integer(i4) :: i_down, j_down
     146             : 
     147           0 :     call cellLength(iDomain, fDir0(i, j), i, j, &
     148           0 :                     iFlag_cordinate_sys, length)
     149           0 :     i_down = i
     150           0 :     j_down = j
     151           0 :     call moveDownOneCell(fDir0(i, j), i_down, j_down)
     152             : 
     153           0 :     slope = (elev0(i,j) - elev0(i_down, j_down)) / length
     154             :     ! TODO: as soon as current gfortran compiler is available on EVE,
     155             :     ! use ieee_isnan from ieee_arithmetic module, instead of
     156             :     ! slope /= slope
     157           0 :     if(slope < 0.0001_dp .OR. slope > 20._dp .OR. slope /= slope) then
     158           0 :       slope = 0.0001_dp
     159             :     end if
     160           0 :   end function calc_slope
     161             : 
     162             : end module mo_mrm_river_head