mo_read_wrapper.f90

Go to the documentation of this file.

!> \file mo_read_wrapper.f90
!> \brief \copybrief mo_read_wrapper
!> \details \copydetails mo_read_wrapper
 
!> \brief Wrapper for all reading routines.
!> \details This module is to wrap up all reading routines.
!! The general written reading routines are used to store now the read data into global variables.
!> \authors Juliane Mai, Matthias Zink
!> \date Jan 2013
!> \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_mpr
MODULE mo_read_wrapper
 
  USE mo_kind, ONLY : i4, dp
  use mo_common_constants, only : nodata_dp, nodata_i4
  use mo_message, only: message, error_message
 
  IMPLICIT NONE
 
  PUBLIC :: read_data            ! reads all available data
 
CONTAINS
 
  ! ------------------------------------------------------------------
 
  !    NAME
  !        read_data
 
  !    PURPOSE
  !>       \brief Reads data.
 
  !>       \details The namelists are already read by read_config call.
  !>       All LUTs are read from their respective directory and information within those
  !>       files are shared across all domains to be modeled.
 
  !    INTENT(IN), OPTIONAL
  !>       \param[in] "type(period), dimension(:), optional :: LAIPer"
 
  !    HISTORY
  !>       \authors Juliane Mai & Matthias Zink
 
  !>       \date Feb 2013
 
  ! Modifications:
  ! Luis Samaniego               Feb 2013 - rotate fdir variable to the new coordinate system
  ! Rohini Kumar                 Aug 2013 - name changed from "L0_LAI" to "L0_LCover_LAI"
  ! Rohini Kumar                 Aug 2013 - added dirSoil_LUT and dirGeology_LUT, and changed to
  !                                         read datapaths and variables made accordingly
  ! Rohini Kumar                 Aug 2013 - added iFlag_LAI_data_format to handle LAI options, and changed
  !                                         within the code made accordingly
  ! Rohini  Kumar                Sep 2013 - read input data for routing processes according to process_matrix flag
  ! Matthias Zink                Mar 2014 - added inflow gauge
  ! Kumar & Schroen              Apr 2014 - added check for consistency of L0 and L1 spatial resolution
  ! Stephan Thober               Jun 2014 - added perform_mpr for omitting L0 read
  ! Matthias Cuntz & Juliane Mai Nov 2014 - LAI input from daily, monthly or yearly files
  ! Stephan Thober               Aug 2015 - moved routing related variables and routines to mRM
  ! Rohini Kumar                 Mar 2016 - options to handle different soil databases
  ! Matthias Zink                Mar 2014 - added subroutine for consistency check
  ! Stephan Thober               Nov 2016 - moved processMatrix to common variables
  ! Rohini Kumar                 Dec 2016 - option to handle monthly mean gridded fields of LAI
  ! Robert Schweppe              Jun 2018 - refactoring and reformatting
 
  subroutine read_data(LAIPer)
 
    use mo_append, only : append, paste
    use mo_constants, only : yearmonths
    use mo_common_read_data, only : read_dem, read_lcover
    use mo_common_types, only: period, grid
    use mo_common_variables, only : dircommonfiles, dirmorpho, &
                                    global_parameters, level0, domainmeta, processmatrix
    use mo_mpr_file, only : file_aspect, file_geolut, file_hydrogeoclass, &
                            file_laiclass, file_lailut, file_slope, file_soil_database, file_soil_database_1, &
                            file_soilclass, uaspect, ugeolut, uhydrogeoclass, ulaiclass, ulailut, uslope, usoilclass
    use mo_mpr_global_variables, only : geounitkar, &
                                        geounitlist, l0_asp, l0_geounit, l0_gridded_lai, l0_slope, l0_soilid, lailut, &
                                        laiunitlist, iflag_soildb, ngeounits, nlai, nlaiclass, nsoilhorizons_mhm, soildb, &
                                        timestep_lai_input, laiboundaries
    use mo_prepare_gridded_lai, only : prepare_gridded_daily_lai_data, prepare_gridded_mean_monthly_lai_data
    use mo_read_latlon, only : read_latlon
    use mo_read_lut, only : read_geoformation_lut, read_lai_lut
    use mo_read_spatial_data, only : read_spatial_data_nc_or_ascii
    use mo_soil_database, only : read_soil_lut
    use mo_string_utils, only : num2str
    use mo_timer, only : timer_get, timer_start, &
                         timer_stop
 
    implicit none
 
    type(period), dimension(:), intent(in), optional :: laiper
 
    ! loop variables
    integer(i4) :: domainid, idomain, ivar, ihorizon, imon, itimer, ll
 
    ! dummy variable
    integer(i4) :: nh
 
    ! file unit of file to read
    integer(i4) :: nunit
 
    ! file name of file to read
    character(256) :: fname
 
    real(dp), dimension(:, :), allocatable :: data_dp_2d
 
    integer(i4), dimension(:, :), allocatable :: data_i4_2d
 
    integer(i4), dimension(:, :), allocatable :: datamatrix_i4
 
    logical, dimension(:, :), allocatable :: mask_2d
 
    integer(i4), dimension(:), allocatable :: dummy_i4
 
    type(grid), pointer :: level0_idomain
 
    real(dp), parameter :: slope_minval = 0.01_dp
 
    real(dp), parameter :: aspect_minval = 1.00_dp
 
 
    call message('  Reading data ...')
    itimer = 1
    call timer_start(itimer)
 
    call message('    Reading dem and lcover ...')
    call read_dem()
    call read_lcover()
 
    ! ************************************************
    ! READ LOOKUP TABLES
    ! ************************************************
    !
    ! Soil LUT
    if(iflag_soildb .eq. 0) then
      fname = trim(adjustl(dircommonfiles)) // trim(adjustl(file_soil_database))
    else if(iflag_soildb .eq. 1) then
      fname = trim(adjustl(dircommonfiles)) // trim(adjustl(file_soil_database_1))
    end if
    call read_soil_lut(trim(fname))
 
    ! Geological formation LUT
    fname = trim(adjustl(dircommonfiles)) // trim(adjustl(file_geolut))
    call read_geoformation_lut(trim(fname), ugeolut, ngeounits, geounitlist, geounitkar)
 
    ! LAI LUT
    if(timestep_lai_input .EQ. 0) then
      fname = trim(adjustl(dircommonfiles)) // trim(adjustl(file_lailut))
      call read_lai_lut(trim(fname), ulailut, nlaiclass, laiunitlist, lailut)
    end if
 
    domains: do idomain = 1, domainmeta%nDomains
      domainid = domainmeta%indices(idomain)
 
      level0_idomain => level0(domainmeta%L0DataFrom(idomain))
 
      call message('    Reading data for domain: ', trim(adjustl(num2str(domainid))), ' ...')
      ! check whether L0 data is shared
      ! ToDo: check change
      if (domainmeta%L0DataFrom(idomain) < idomain) then
        !
        call message('      Using data of domain ', &
                trim(adjustl(num2str(domainmeta%indices(domainmeta%L0DataFrom(idomain))))), ' for domain: ',&
                trim(adjustl(num2str(domainid))), '...')
        ! DO NOT read L0 data
        cycle
 
      end if
 
      itimer = 2
      call timer_start(itimer)
 
      ! read slope and aspect - datatype real
      nvars_real : do ivar = 1, 2
        select case (ivar)
        case(1) ! slope
          call message('      Reading slope ...')
          fname = trim(adjustl(dirmorpho(idomain))) // trim(adjustl(file_slope))
          nunit = uslope
        case(2) ! aspect
          call message('      Reading aspect ...')
          fname = trim(adjustl(dirmorpho(idomain))) // trim(adjustl(file_aspect))
          nunit = uaspect
        end select
 
        ! reading
        call read_spatial_data_nc_or_ascii(trim(fname), nunit, &
                level0_idomain%nrows, level0_idomain%ncols, level0_idomain%xllcorner, &
                level0_idomain%yllcorner, level0_idomain%cellsize, data_dp_2d, mask_2d)
        ! put global nodata value into array (probably not all grid cells have values)
        data_dp_2d = merge(data_dp_2d, nodata_dp, mask_2d)
        ! put data in variable
        select case (ivar)
        case(1) ! slope
          call append(l0_slope, pack(data_dp_2d, level0_idomain%mask))
 
        case(2) ! aspect
          call append(l0_asp, pack(data_dp_2d, level0_idomain%mask))
        end select
        ! deallocate arrays
        deallocate(data_dp_2d, mask_2d)
 
      end do nvars_real
 
      ! read datatype integer
      ! ***** CHANGE IS MADE HERE TO ACCOMODATE READ SEVERAL TYPES OF SOIL DATABASES
      ! change from the earlier code where everything was done in the DO loop (**see below **)
      ! here everything is more explicit and no do loop appears as was the case in previous version
 
      ! read soilID in both options
      call message('      Reading soil ids ...')
 
      nh = 1 !> by default; when iFlag_soilDB = 0
      if(iflag_soildb .eq. 1) nh = nsoilhorizons_mhm
      ! modified way to read multiple horizons specific soil class
      do ihorizon = 1, nh
        if(iflag_soildb .eq. 0) then
          fname = trim(adjustl(dirmorpho(idomain))) // trim(adjustl(file_soilclass))
        else if(iflag_soildb .eq. 1) then
          write(fname, 172) ihorizon
          172             format('soil_class_horizon_', i2.2, '.asc')
          fname = trim(adjustl(dirmorpho(idomain))) // trim(adjustl(fname))
        end if
        call read_spatial_data_nc_or_ascii(trim(fname), usoilclass, &
                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_soilid, datamatrix_i4)
      deallocate(datamatrix_i4)
 
      ! read geoUnit
      fname = trim(adjustl(dirmorpho(idomain))) // trim(adjustl(file_hydrogeoclass))
      ! reading and transposing
      call read_spatial_data_nc_or_ascii(trim(fname), uhydrogeoclass, &
              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 append(l0_geounit, pack(data_i4_2d, level0_idomain%mask))
      deallocate(data_i4_2d, mask_2d)
 
      ! LAI values
      call message('      Reading LAI ...')
      select case (timestep_lai_input)
      case(1) ! long term mean monthly gridded fields
        call prepare_gridded_mean_monthly_lai_data(idomain, level0_idomain%nrows, level0_idomain%ncols, level0_idomain%mask)
 
      case(0) ! long term mean monthly values per class with LUT
        ! only set if not yet allocated (e.g. domain 1)
        if (.not. allocated(laiboundaries)) then
          nlai = int(yearmonths, i4)
          allocate(laiboundaries(nlai+1))
          laiboundaries = [(imon, imon=1, nlai+1)]
        end if
 
        fname = trim(adjustl(dirmorpho(idomain))) // trim(adjustl(file_laiclass))
        ! reading and transposing
        call read_spatial_data_nc_or_ascii(trim(fname), ulaiclass, &
                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)
        allocate(dummy_i4(count(level0_idomain%mask)))
        dummy_i4 = pack(data_i4_2d, level0_idomain%mask)
        deallocate(data_i4_2d, mask_2d)
 
        call check_consistency_lut_map(dummy_i4, laiunitlist, file_laiclass)
 
        allocate(data_dp_2d(count(level0_idomain%mask), nlai))
        do imon = 1, nlai
          ! determine LAIs per month
          do ll = 1, size(lailut, dim = 1)
            data_dp_2d(:, imon) = merge(lailut(ll, imon), data_dp_2d(:, imon), dummy_i4(:) .EQ. laiunitlist(ll))
          end do
        end do
        call append(l0_gridded_lai, data_dp_2d(:, :))
        deallocate(dummy_i4, data_dp_2d)
        ! correction for 0 LAI values to avoid numerical instabilities
        l0_gridded_lai(:, :) = merge(1.00e-10_dp, l0_gridded_lai(:, :), l0_gridded_lai(:, :) .LT. 1.00e-10_dp)
        l0_gridded_lai(:, :) = merge(30.0_dp, l0_gridded_lai(:, :), l0_gridded_lai(:, :) .GT. 30.0_dp)
      case(-3 : -1) ! daily, monthly or yearly gridded fields (time-series)
        call prepare_gridded_daily_lai_data(idomain, level0_idomain%nrows, level0_idomain%ncols, level0_idomain%mask, &
        laiper(idomain))
 
      end select
 
      ! read lat lon coordinates of each domain
      call message('      Reading latitude/logitude ...')
      call read_latlon(idomain, "lon_l0", "lat_l0", "level0", level0_idomain)
 
      call timer_stop(itimer)
      call message('    in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
 
    end do domains
    itimer = 1
    call timer_stop(itimer)
    call message('  in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
 
    !Soil
    ! determine name od soil class definition file based on input option iFlag_soilDB
    if(iflag_soildb .eq. 0) then
      fname = file_soil_database
    else if(iflag_soildb .eq. 1) then
      fname = file_soil_database_1
    end if
 
    ! If you are getting segmentation fault with intel on large file (e.g. here on soils), then
    ! make sure to include following into your Marefile: INTEL_EXCLUDE  := mo_read_wrapper.f90
    call check_consistency_lut_map(reshape(l0_soilid, (/ size(l0_soilid, 1) * size(l0_soilid, 2) /)), &
            soildb%id(:), fname)
    ! Geology
    call check_consistency_lut_map(l0_geounit, geounitlist, file_hydrogeoclass, dummy_i4)
 
    ! deactivate parameters of non existing geological classes in study domain for optimization
    ! loop over geological units in look up list
    do ivar = 1, size(geounitlist, 1)
      ! check if unit appears in geological map (dummy_i4 is unique number in L0_geoUnit)
      if (.not. any(dummy_i4 .EQ. geounitlist(ivar))) then
        ! deactivate optimization flag (dim=4 from global_parameters)
        global_parameters(processmatrix(9, 3) - processmatrix(9, 2) + ivar, 4) = 0
        call message('***WARNING: Geological unit ', trim(adjustl(num2str(geounitlist(ivar)))))
        call message('            is not appearing in study domain.')
      end if
    end do
 
    deallocate(dummy_i4) ! is allocated in subroutine check_consistency_lut_map - geology
 
    !----------------------------------------------------------------
    ! Correction for slope and aspect -- min value set above
    !----------------------------------------------------------------
    ! keep the colons (:) in the statements because of Intel's reallocation lhs problem
    l0_slope(:) = merge(slope_minval, l0_slope(:), (l0_slope(:)  .lt.  slope_minval))
    l0_asp(:) = merge(aspect_minval, l0_asp(:), (l0_asp(:)    .lt. aspect_minval))
 
  end subroutine read_data
 
  ! ------------------------------------------------------------------
 
  !    NAME
  !        check_consistency_lut_map
 
  !    PURPOSE
  !>       \brief Checks if classes in input maps appear in look up tables.
 
  !>       \details Determines wether a class appearing in the morphological input
  !>       is occuring in the respective look up table. mHM breaks if inconsistencies
  !>       are discovered.
 
  !    INTENT(IN)
  !>       \param[in] "integer(i4), dimension(:) :: data"        map of study domain
  !>       \param[in] "integer(i4), dimension(:) :: lookuptable" look up table corresponding to map
  !>       \param[in] "character(*) :: filename"                 name of the lut file - ERORR warn
 
  !    INTENT(OUT), OPTIONAL
  !>       \param[out] "integer(i4), dimension(:), optional :: unique_values" array of unique values in dataone
 
  !    HISTORY
  !>       \authors Matthias Zink
 
  !>       \date Nov 2016
 
  ! Modifications:
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  subroutine check_consistency_lut_map(data, lookuptable, filename, unique_values)
 
    use mo_orderpack, only : unista
    use mo_string_utils, only : num2str
 
    implicit none
 
    ! map of study domain
    integer(i4), dimension(:), intent(in) :: data
 
    ! look up table corresponding to map
    integer(i4), dimension(:), intent(in) :: lookuptable
 
    ! name of the lut file - ERORR warn
    character(*), intent(in) :: filename
 
    ! array of unique values in dataone
    integer(i4), dimension(:), allocatable, intent(out), optional :: unique_values
 
    integer(i4) :: n_unique_elements
 
    integer(i4) :: ielement
 
    integer(i4), dimension(:), allocatable :: temp
 
 
    allocate(temp(size(data, 1)))
    ! copy L0_geoUnit because subroutine unista overwrites the nunit entries of the
    ! input array with the unique array values
    temp = data
    ! retrieve unique values of data
    call unista(temp, n_unique_elements)
    ! check if unit exists in look up table
    do ielement = 1, n_unique_elements
      if ( temp(ielement) == nodata_i4 ) call error_message( &
        '***ERROR: Class ', trim(adjustl(num2str(temp(ielement)))), &
        ' was searched in ', trim(adjustl(filename)), &
        ' which indicates a masking problem!' &
      )
      if (.not. any(lookuptable .EQ. temp(ielement))) then
        call error_message('***ERROR: Class ', trim(adjustl(num2str(temp(ielement)))), ' is missing', raise=.false.)
        call error_message('          in input file ', trim(adjustl(filename)), ' ...')
      end if
    end do
 
    ! pass unique values if optional argument unique_values is given
    if (present(unique_values)) then
      allocate(unique_values(n_unique_elements))
      unique_values(:) = temp(1 : n_unique_elements)
    end if
    deallocate(temp)
 
  end subroutine check_consistency_lut_map
 
 
END MODULE mo_read_wrapper