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