Line data    Source code

       1             : !> \file mo_common_read_data.f90
       2             : !> \brief   \copybrief mo_common_read_data
       3             : !> \details \copydetails mo_common_read_data
       4             : 
       5             : !> \brief Common reading routines
       6             : !> \details Routines to read the DEM and landcover files.
       7             : !> \authors Robert Schweppe
       8             : !> \date Jun 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_common
      12             : module mo_common_read_data
      13             :   USE mo_kind, ONLY : i4, dp
      14             :   use mo_message, only: message, error_message
      15             : 
      16             :   IMPLICIT NONE
      17             : 
      18             :   PRIVATE
      19             : 
      20             :   PUBLIC :: read_dem, read_lcover
      21             : 
      22             :   ! ------------------------------------------------------------------
      23             : 
      24             : CONTAINS
      25             : 
      26             :   !    NAME
      27             :   !        read_dem
      28             : 
      29             :   !    PURPOSE
      30             :   !>       \brief TODO: add description
      31             : 
      32             :   !>       \details TODO: add description
      33             : 
      34             :   !    HISTORY
      35             :   !>       \authors Robert Schweppe
      36             : 
      37             :   !>       \date Jun 2018
      38             : 
      39             :   ! Modifications:
      40             : 
      41          13 :   subroutine read_dem
      42             : 
      43             :     use mo_append, only : append
      44             :     use mo_common_constants, only : nodata_dp
      45             :     use mo_common_file, only : file_dem, udem
      46             :     use mo_common_types, only: Grid
      47             :     use mo_common_variables, only : L0_elev, dirMorpho, level0, domainMeta, &
      48             :                                     resolutionHydrology
      49             :     use mo_grid, only : set_domain_indices
      50             :     use mo_read_spatial_data, only : read_header_ascii, read_spatial_data_ascii
      51             :     use mo_string_utils, only : num2str
      52             : 
      53             :     implicit none
      54             : 
      55             :     ! loop variables
      56             :     integer(i4) :: domainID, iDomain
      57             : 
      58             :     ! file name of file to read
      59             :     character(256) :: fName
      60             : 
      61          13 :     real(dp), dimension(:, :), allocatable :: data_dp_2d
      62             : 
      63             :     type(Grid), pointer :: level0_iDomain
      64             : 
      65             : 
      66             :     ! ************************************************
      67             :     ! READ SPATIAL DATA FOR EACH DOMAIN
      68             :     ! ************************************************
      69             :     ! allocate necessary variables at Level0
      70          64 :     allocate(level0(domainMeta%nDomains))
      71             : 
      72          38 :     do iDomain = 1, domainMeta%nDomains
      73          25 :       domainID = domainMeta%indices(iDomain)
      74             : 
      75          25 :       level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
      76             : 
      77             :       ! check whether L0 data is shared
      78          25 :       if (iDomain .gt. 1) then
      79          12 :         if (domainMeta%L0DataFrom(iDomain) < iDomain) then
      80             :           !
      81             :           call message('      Using dem of domain ', &
      82           3 :                   trim(adjustl(num2str(domainMeta%indices(domainMeta%L0DataFrom(iDomain))))), ' for domain: ',&
      83           3 :                   trim(adjustl(num2str(iDomain))), '...')
      84             : 
      85             :           ! DO NOT read L0 data
      86           3 :           cycle
      87             : 
      88             :         end if
      89             :       end if
      90             : 
      91          22 :       call message('      Reading dem for domain: ', trim(adjustl(num2str(domainID))), ' ...')
      92             : 
      93             :       ! Header (to check consistency)
      94          22 :       fName = trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(file_dem))
      95             :       call read_header_ascii(trim(fName), udem, &
      96             :               level0_iDomain%nrows, level0_iDomain%ncols, level0_iDomain%xllcorner, &
      97          22 :               level0_iDomain%yllcorner, level0_iDomain%cellsize, level0_iDomain%nodata_value)
      98             : 
      99             :       ! check for L0 and L1 scale consistency
     100          22 :       if(resolutionHydrology(iDomain) .LT. level0_iDomain%cellsize) then
     101           0 :         call error_message('***ERROR: resolutionHydrology (L1) should be smaller than the input data resolution (L0)', raise=.false.)
     102           0 :         call error_message('          check set-up (in mhm.nml) for domain: ', trim(adjustl(num2str(domainID))), ' ...')
     103             :       end if
     104             : 
     105             :       ! DEM + overall mask creation
     106          22 :       fName = trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(file_dem))
     107             :       call read_spatial_data_ascii(trim(fName), udem, &
     108             :               level0_iDomain%nrows, level0_iDomain%ncols, level0_iDomain%xllcorner, &
     109          22 :               level0_iDomain%yllcorner, level0_iDomain%cellsize, data_dp_2d, level0_iDomain%mask)
     110             : 
     111             :       ! put global nodata value into array (probably not all grid cells have values)
     112     2612684 :       data_dp_2d = merge(data_dp_2d, nodata_dp, level0_iDomain%mask)
     113             :       ! put data in variable
     114          22 :       call append(L0_elev, pack(data_dp_2d, level0_iDomain%mask))
     115             :       ! deallocate arrays
     116          22 :       deallocate(data_dp_2d)
     117             : 
     118     2612675 :       level0_iDomain%nCells = count(level0_iDomain%mask)
     119             : 
     120             :     end do
     121             : 
     122             :     call set_domain_indices(level0, indices=domainMeta%L0DataFrom)
     123             : 
     124          13 :   end subroutine read_dem
     125             : 
     126             :   !    NAME
     127             :   !        read_lcover
     128             : 
     129             :   !    PURPOSE
     130             :   !>       \brief TODO: add description
     131             : 
     132             :   !>       \details TODO: add description
     133             : 
     134             :   !    HISTORY
     135             :   !>       \authors Robert Schweppe
     136             : 
     137             :   !>       \date Jun 2018
     138             : 
     139             :   ! Modifications:
     140             : 
     141          22 :   subroutine read_lcover
     142             : 
     143          13 :     use mo_append, only : append, paste
     144             :     use mo_common_constants, only : nodata_i4
     145             :     use mo_common_file, only : ulcoverclass
     146             :     use mo_common_types, only: Grid
     147             :     use mo_common_variables, only : L0_LCover, LCfilename, dirLCover, level0, domainMeta, nLCoverScene
     148             :     use mo_read_spatial_data, only : read_spatial_data_ascii
     149             :     use mo_string_utils, only : num2str
     150             : 
     151             :     implicit none
     152             : 
     153             :     ! loop variables
     154             :     integer(i4) :: domainID, iDomain, iVar
     155             : 
     156             :     ! file name of file to read
     157             :     character(256) :: fName
     158             : 
     159          22 :     integer(i4), dimension(:, :), allocatable :: data_i4_2d
     160             : 
     161          22 :     integer(i4), dimension(:, :), allocatable :: dataMatrix_i4
     162             : 
     163          22 :     logical, dimension(:, :), allocatable :: mask_2d
     164             : 
     165             :     type(Grid), pointer :: level0_iDomain
     166             : 
     167             : 
     168          63 :     do iDomain = 1, domainMeta%nDomains
     169          41 :       domainID = domainMeta%indices(iDomain)
     170             : 
     171          41 :       level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
     172             : 
     173             :       ! check whether L0 data is shared
     174             :       ! ToDo: check change
     175          41 :       if (domainMeta%L0DataFrom(iDomain) < iDomain) then
     176             :         call message('      Using lcover of domain ', &
     177             :                 trim(adjustl(num2str(domainMeta%L0DataFrom(iDomain)))), ' for domain: ',&
     178           6 :                 trim(adjustl(num2str(domainID))), '...')
     179             :         ! DO NOT read L0 data
     180           6 :         cycle
     181             : 
     182             :       end if
     183             : 
     184          35 :       call message('      Reading lcover for domain: ', trim(adjustl(num2str(domainID))), ' ...')
     185             : 
     186             :       ! LCover read in is realized seperated because of unknown number of scenes
     187         105 :       do iVar = 1, nLCoverScene
     188          70 :         fName = trim(adjustl(dirLCover(iDomain))) // trim(adjustl(LCfilename(iVar)))
     189             :         call read_spatial_data_ascii(trim(fName), ulcoverclass, &
     190             :                 level0_iDomain%nrows, level0_iDomain%ncols, level0_iDomain%xllcorner, &
     191          70 :                 level0_iDomain%yllcorner, level0_iDomain%cellsize, data_i4_2d, mask_2d)
     192             :         ! put global nodata value into array (probably not all grid cells have values)
     193     8203436 :         data_i4_2d = merge(data_i4_2d, nodata_i4, mask_2d)
     194          70 :         call paste(dataMatrix_i4, pack(data_i4_2d, level0_iDomain%mask), nodata_i4)
     195         105 :         deallocate(data_i4_2d)
     196             :       end do
     197          35 :       call append(L0_LCover, dataMatrix_i4)
     198          57 :       deallocate(dataMatrix_i4)
     199             : 
     200             :     end do
     201             : 
     202          22 :   end subroutine read_lcover
     203             : 
     204             : end module mo_common_read_data