mo_common_read_data.f90

Go to the documentation of this file.

!> \file mo_common_read_data.f90
!> \brief   \copybrief mo_common_read_data
!> \details \copydetails mo_common_read_data
 
!> \brief Common reading routines
!> \details Routines to read the DEM and landcover files.
!> \authors Robert Schweppe
!> \date Jun 2018
!> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
!! mHM is released under the LGPLv3+ license \license_note
!> \ingroup f_common
module mo_common_read_data
  USE mo_kind, ONLY : i4, dp
  use mo_message, only: message, error_message
  use mo_read_spatial_data, only : read_spatial_data_nc_or_ascii
 
  IMPLICIT NONE
 
  PRIVATE
 
  PUBLIC :: read_dem, read_lcover
 
  ! ------------------------------------------------------------------
 
CONTAINS
 
  !    NAME
  !        read_dem
 
  !    PURPOSE
  !>       \brief TODO: add description
 
  !>       \details TODO: add description
 
  !    HISTORY
  !>       \authors Robert Schweppe
 
  !>       \date Jun 2018
 
  ! Modifications:
 
  subroutine read_dem
 
    use mo_append, only : append
    use mo_common_constants, only : nodata_dp
    use mo_common_file, only : file_dem, udem
    use mo_common_types, only: grid
    use mo_common_variables, only : l0_elev, dirmorpho, level0, domainmeta, &
                                    resolutionhydrology
    use mo_common_grid, only : set_domain_indices
    use mo_read_spatial_data, only : read_header_ascii, read_spatial_data_ascii
    use mo_string_utils, only : num2str
 
    implicit none
 
    ! loop variables
    integer(i4) :: domainid, idomain
 
    ! file name of file to read
    character(256) :: fname
 
    real(dp), dimension(:, :), allocatable :: data_dp_2d
 
    type(grid), pointer :: level0_idomain
 
 
    ! ************************************************
    ! READ SPATIAL DATA FOR EACH DOMAIN
    ! ************************************************
    ! allocate necessary variables at Level0
    allocate(level0(domainmeta%nDomains))
 
    do idomain = 1, domainmeta%nDomains
      domainid = domainmeta%indices(idomain)
 
      level0_idomain => level0(domainmeta%L0DataFrom(idomain))
 
      ! check whether L0 data is shared
      if (idomain .gt. 1) then
        if (domainmeta%L0DataFrom(idomain) < idomain) then
          !
          call message('      Using dem of domain ', &
                  trim(adjustl(num2str(domainmeta%indices(domainmeta%L0DataFrom(idomain))))), ' for domain: ',&
                  trim(adjustl(num2str(idomain))), '...')
 
          ! DO NOT read L0 data
          cycle
 
        end if
      end if
 
      call message('      Reading dem for domain: ', trim(adjustl(num2str(domainid))), ' ...')
 
      ! Header (to check consistency)
      fname = trim(adjustl(dirmorpho(idomain))) // trim(adjustl(file_dem))
      call read_header_ascii(trim(fname), udem, &
              level0_idomain%nrows, level0_idomain%ncols, level0_idomain%xllcorner, &
              level0_idomain%yllcorner, level0_idomain%cellsize, level0_idomain%nodata_value)
 
      ! check for L0 and L1 scale consistency
      if(resolutionhydrology(idomain) .LT. level0_idomain%cellsize) then
        call error_message('***ERROR: resolutionHydrology (L1) should be smaller than the input data resolution (L0)', raise=.false.)
        call error_message('          check set-up (in mhm.nml) for domain: ', trim(adjustl(num2str(domainid))), ' ...')
      end if
 
      ! DEM + overall mask creation
      fname = trim(adjustl(dirmorpho(idomain))) // trim(adjustl(file_dem))
      call read_spatial_data_nc_or_ascii(trim(fname), udem, &
              level0_idomain%nrows, level0_idomain%ncols, level0_idomain%xllcorner, &
              level0_idomain%yllcorner, level0_idomain%cellsize, data_dp_2d, level0_idomain%mask)
 
      ! put global nodata value into array (probably not all grid cells have values)
      data_dp_2d = merge(data_dp_2d, nodata_dp, level0_idomain%mask)
      ! put data in variable
      call append(l0_elev, pack(data_dp_2d, level0_idomain%mask))
      ! deallocate arrays
      deallocate(data_dp_2d)
 
      level0_idomain%nCells = count(level0_idomain%mask)
 
    end do
 
    call set_domain_indices(level0, indices=domainmeta%L0DataFrom)
 
  end subroutine read_dem
 
  !    NAME
  !        read_lcover
 
  !    PURPOSE
  !>       \brief TODO: add description
 
  !>       \details TODO: add description
 
  !    HISTORY
  !>       \authors Robert Schweppe
 
  !>       \date Jun 2018
 
  ! Modifications:
 
  subroutine read_lcover
 
    use mo_append, only : append, paste
    use mo_common_constants, only : nodata_i4
    use mo_common_file, only : ulcoverclass
    use mo_common_types, only: grid
    use mo_common_variables, only : l0_lcover, lcfilename, dirlcover, level0, domainmeta, nlcoverscene
    use mo_read_spatial_data, only : read_spatial_data_ascii
    use mo_string_utils, only : num2str
 
    implicit none
 
    ! loop variables
    integer(i4) :: domainid, idomain, ivar
 
    ! file name of file to read
    character(256) :: fname
 
    integer(i4), dimension(:, :), allocatable :: data_i4_2d
 
    integer(i4), dimension(:, :), allocatable :: datamatrix_i4
 
    logical, dimension(:, :), allocatable :: mask_2d
 
    type(grid), pointer :: level0_idomain
 
 
    do idomain = 1, domainmeta%nDomains
      domainid = domainmeta%indices(idomain)
 
      level0_idomain => level0(domainmeta%L0DataFrom(idomain))
 
      ! check whether L0 data is shared
      ! ToDo: check change
      if (domainmeta%L0DataFrom(idomain) < idomain) then
        call message('      Using lcover of domain ', &
                trim(adjustl(num2str(domainmeta%L0DataFrom(idomain)))), ' for domain: ',&
                trim(adjustl(num2str(domainid))), '...')
        ! DO NOT read L0 data
        cycle
 
      end if
 
      call message('      Reading lcover for domain: ', trim(adjustl(num2str(domainid))), ' ...')
 
      ! LCover read in is realized seperated because of unknown number of scenes
      do ivar = 1, nlcoverscene
        fname = trim(adjustl(dirlcover(idomain))) // trim(adjustl(lcfilename(ivar)))
        call read_spatial_data_nc_or_ascii(trim(fname), ulcoverclass, &
                level0_idomain%nrows, level0_idomain%ncols, level0_idomain%xllcorner, &
                level0_idomain%yllcorner, level0_idomain%cellsize, data_i4_2d, mask_2d)
        ! put global nodata value into array (probably not all grid cells have values)
        data_i4_2d = merge(data_i4_2d, nodata_i4, mask_2d)
        call paste(datamatrix_i4, pack(data_i4_2d, level0_idomain%mask), nodata_i4)
        deallocate(data_i4_2d)
      end do
      call append(l0_lcover, datamatrix_i4)
      deallocate(datamatrix_i4)
 
    end do
 
  end subroutine read_lcover
 
end module mo_common_read_data