LCOV - code coverage report
Current view: top level - MPR - mo_read_wrapper.f90 (source / functions) Hit Total Coverage
Test: mHM coverage Lines: 113 126 89.7 %
Date: 2023-08-24 16:49:59 Functions: 2 2 100.0 %

          Line data    Source code
       1             : !> \file mo_read_wrapper.f90
       2             : !> \brief \copybrief mo_read_wrapper
       3             : !> \details \copydetails mo_read_wrapper
       4             : 
       5             : !> \brief Wrapper for all reading routines.
       6             : !> \details This module is to wrap up all reading routines.
       7             : !! The general written reading routines are used to store now the read data into global variables.
       8             : !> \authors Juliane Mai, Matthias Zink
       9             : !> \date Jan 2013
      10             : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
      11             : !! mHM is released under the LGPLv3+ license \license_note
      12             : !> \ingroup f_mpr
      13             : MODULE mo_read_wrapper
      14             : 
      15             :   USE mo_kind, ONLY : i4, dp
      16             :   use mo_common_constants, only : nodata_dp, nodata_i4
      17             :   use mo_message, only: message, error_message
      18             : 
      19             :   IMPLICIT NONE
      20             : 
      21             :   PUBLIC :: read_data            ! reads all available data
      22             : 
      23             : CONTAINS
      24             : 
      25             :   ! ------------------------------------------------------------------
      26             : 
      27             :   !    NAME
      28             :   !        read_data
      29             : 
      30             :   !    PURPOSE
      31             :   !>       \brief Reads data.
      32             : 
      33             :   !>       \details The namelists are already read by read_config call.
      34             :   !>       All LUTs are read from their respective directory and information within those
      35             :   !>       files are shared across all domains to be modeled.
      36             : 
      37             :   !    INTENT(IN), OPTIONAL
      38             :   !>       \param[in] "type(period), dimension(:), optional :: LAIPer"
      39             : 
      40             :   !    HISTORY
      41             :   !>       \authors Juliane Mai & Matthias Zink
      42             : 
      43             :   !>       \date Feb 2013
      44             : 
      45             :   ! Modifications:
      46             :   ! Luis Samaniego               Feb 2013 - rotate fdir variable to the new coordinate system
      47             :   ! Rohini Kumar                 Aug 2013 - name changed from "L0_LAI" to "L0_LCover_LAI"
      48             :   ! Rohini Kumar                 Aug 2013 - added dirSoil_LUT and dirGeology_LUT, and changed to
      49             :   !                                         read datapaths and variables made accordingly
      50             :   ! Rohini Kumar                 Aug 2013 - added iFlag_LAI_data_format to handle LAI options, and changed
      51             :   !                                         within the code made accordingly
      52             :   ! Rohini  Kumar                Sep 2013 - read input data for routing processes according to process_matrix flag
      53             :   ! Matthias Zink                Mar 2014 - added inflow gauge
      54             :   ! Kumar & Schroen              Apr 2014 - added check for consistency of L0 and L1 spatial resolution
      55             :   ! Stephan Thober               Jun 2014 - added perform_mpr for omitting L0 read
      56             :   ! Matthias Cuntz & Juliane Mai Nov 2014 - LAI input from daily, monthly or yearly files
      57             :   ! Stephan Thober               Aug 2015 - moved routing related variables and routines to mRM
      58             :   ! Rohini Kumar                 Mar 2016 - options to handle different soil databases
      59             :   ! Matthias Zink                Mar 2014 - added subroutine for consistency check
      60             :   ! Stephan Thober               Nov 2016 - moved processMatrix to common variables
      61             :   ! Rohini Kumar                 Dec 2016 - option to handle monthly mean gridded fields of LAI
      62             :   ! Robert Schweppe              Jun 2018 - refactoring and reformatting
      63             : 
      64          26 :   subroutine read_data(LAIPer)
      65             : 
      66             :     use mo_append, only : append, paste
      67             :     use mo_constants, only : YearMonths
      68             :     use mo_common_read_data, only : read_dem, read_lcover
      69             :     use mo_common_types, only: period, Grid
      70             :     use mo_common_variables, only : dirCommonFiles, dirMorpho, &
      71             :                                     global_parameters, level0, domainMeta, processMatrix
      72             :     use mo_mpr_file, only : file_aspect, file_geolut, file_hydrogeoclass, &
      73             :                             file_laiclass, file_lailut, file_slope, file_soil_database, file_soil_database_1, &
      74             :                             file_soilclass, uaspect, ugeolut, uhydrogeoclass, ulaiclass, ulailut, uslope, usoilclass
      75             :     use mo_mpr_global_variables, only : GeoUnitKar, &
      76             :                                         GeoUnitList, L0_asp, L0_geoUnit, L0_gridded_LAI, L0_slope, L0_soilId, LAILUT, &
      77             :                                         LAIUnitList, iFlag_soilDB, nGeoUnits, nLAI, nLAIclass, nSoilHorizons_mHM, soilDB, &
      78             :                                         timeStep_LAI_input, LAIBoundaries
      79             :     use mo_prepare_gridded_lai, only : prepare_gridded_daily_LAI_data, prepare_gridded_mean_monthly_LAI_data
      80             :     use mo_read_latlon, only : read_latlon
      81             :     use mo_read_lut, only : read_geoformation_lut, read_lai_lut
      82             :     use mo_read_spatial_data, only : read_spatial_data_ascii
      83             :     use mo_soil_database, only : read_soil_LUT
      84             :     use mo_string_utils, only : num2str
      85             :     use mo_timer, only : timer_get, timer_start, &
      86             :                          timer_stop
      87             : 
      88             :     implicit none
      89             : 
      90             :     type(period), dimension(:), intent(in), optional :: LAIPer
      91             : 
      92             :     ! loop variables
      93             :     integer(i4) :: domainID, iDomain, iVar, iHorizon, iMon, itimer, ll
      94             : 
      95             :     ! dummy variable
      96             :     integer(i4) :: nH
      97             : 
      98             :     ! file unit of file to read
      99             :     integer(i4) :: nunit
     100             : 
     101             :     ! file name of file to read
     102             :     character(256) :: fName
     103             : 
     104          13 :     real(dp), dimension(:, :), allocatable :: data_dp_2d
     105             : 
     106          13 :     integer(i4), dimension(:, :), allocatable :: data_i4_2d
     107             : 
     108          13 :     integer(i4), dimension(:, :), allocatable :: dataMatrix_i4
     109             : 
     110          13 :     logical, dimension(:, :), allocatable :: mask_2d
     111             : 
     112          13 :     integer(i4), dimension(:), allocatable :: dummy_i4
     113             : 
     114             :     type(Grid), pointer :: level0_iDomain
     115             : 
     116             :     real(dp), parameter :: slope_minVal = 0.01_dp
     117             : 
     118             :     real(dp), parameter :: aspect_minVal = 1.00_dp
     119             : 
     120             : 
     121          13 :     call message('  Reading data ...')
     122          13 :     itimer = 1
     123          13 :     call timer_start(itimer)
     124             : 
     125          13 :     call message('    Reading dem and lcover ...')
     126          13 :     call read_dem()
     127          13 :     call read_lcover()
     128             : 
     129             :     ! ************************************************
     130             :     ! READ LOOKUP TABLES
     131             :     ! ************************************************
     132             :     !
     133             :     ! Soil LUT
     134          13 :     if(iFlag_soilDB .eq. 0) then
     135          13 :       fName = trim(adjustl(dirCommonFiles)) // trim(adjustl(file_soil_database))
     136           0 :     else if(iFlag_soilDB .eq. 1) then
     137           0 :       fName = trim(adjustl(dirCommonFiles)) // trim(adjustl(file_soil_database_1))
     138             :     end if
     139          13 :     call read_soil_LUT(trim(fName))
     140             : 
     141             :     ! Geological formation LUT
     142          13 :     fName = trim(adjustl(dirCommonFiles)) // trim(adjustl(file_geolut))
     143          13 :     call read_geoformation_lut(trim(fName), ugeolut, nGeoUnits, GeoUnitList, GeoUnitKar)
     144             : 
     145             :     ! LAI LUT
     146          13 :     if(timeStep_LAI_input .EQ. 0) then
     147          12 :       fName = trim(adjustl(dirCommonFiles)) // trim(adjustl(file_lailut))
     148          12 :       call read_lai_lut(trim(fName), ulailut, nLAIclass, LAIUnitList, LAILUT)
     149             :     end if
     150             : 
     151          38 :     domains: do iDomain = 1, domainMeta%nDomains
     152          25 :       domainID = domainMeta%indices(iDomain)
     153             : 
     154          25 :       level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
     155             : 
     156          25 :       call message('    Reading data for domain: ', trim(adjustl(num2str(domainID))), ' ...')
     157             :       ! check whether L0 data is shared
     158             :       ! ToDo: check change
     159          25 :       if (domainMeta%L0DataFrom(iDomain) < iDomain) then
     160             :         !
     161             :         call message('      Using data of domain ', &
     162           0 :                 trim(adjustl(num2str(domainMeta%indices(domainMeta%L0DataFrom(iDomain))))), ' for domain: ',&
     163           3 :                 trim(adjustl(num2str(domainID))), '...')
     164             :         ! DO NOT read L0 data
     165           3 :         cycle
     166             : 
     167             :       end if
     168             : 
     169          22 :       itimer = 2
     170          22 :       call timer_start(itimer)
     171             : 
     172             :       ! read slope and aspect - datatype real
     173          66 :       nVars_real : do iVar = 1, 2
     174          22 :         select case (iVar)
     175             :         case(1) ! slope
     176          22 :           call message('      Reading slope ...')
     177          22 :           fName = trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(file_slope))
     178          22 :           nunit = uslope
     179             :         case(2) ! aspect
     180          22 :           call message('      Reading aspect ...')
     181          22 :           fName = trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(file_aspect))
     182          66 :           nunit = uaspect
     183             :         end select
     184             : 
     185             :         ! reading
     186             :         call read_spatial_data_ascii(trim(fName), nunit, &
     187             :                 level0_iDomain%nrows, level0_iDomain%ncols, level0_iDomain%xllcorner, &
     188          44 :                 level0_iDomain%yllcorner, level0_iDomain%cellsize, data_dp_2d, mask_2d)
     189             :         ! put global nodata value into array (probably not all grid cells have values)
     190     5225368 :         data_dp_2d = merge(data_dp_2d, nodata_dp, mask_2d)
     191             :         ! put data in variable
     192             :         select case (iVar)
     193             :         case(1) ! slope
     194          22 :           call append(L0_slope, pack(data_dp_2d, level0_iDomain%mask))
     195             : 
     196             :         case(2) ! aspect
     197          44 :           call append(L0_asp, pack(data_dp_2d, level0_iDomain%mask))
     198             :         end select
     199             :         ! deallocate arrays
     200          66 :         deallocate(data_dp_2d, mask_2d)
     201             : 
     202             :       end do nVars_real
     203             : 
     204             :       ! read datatype integer
     205             :       ! ***** CHANGE IS MADE HERE TO ACCOMODATE READ SEVERAL TYPES OF SOIL DATABASES
     206             :       ! change from the earlier code where everything was done in the DO loop (**see below **)
     207             :       ! here everything is more explicit and no do loop appears as was the case in previous version
     208             : 
     209             :       ! read soilID in both options
     210          22 :       call message('      Reading soil ids ...')
     211             : 
     212          22 :       nH = 1 !> by default; when iFlag_soilDB = 0
     213          22 :       if(iFlag_soilDB .eq. 1) nH = nSoilHorizons_mHM
     214             :       ! modified way to read multiple horizons specific soil class
     215          44 :       do iHorizon = 1, nH
     216          22 :         if(iFlag_soilDB .eq. 0) then
     217          22 :           fName = trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(file_soilclass))
     218           0 :         else if(iFlag_soilDB .eq. 1) then
     219           0 :           write(fName, 172) iHorizon
     220             :           172             format('soil_class_horizon_', i2.2, '.asc')
     221           0 :           fName = trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(fName))
     222             :         end if
     223             :         call read_spatial_data_ascii(trim(fName), usoilclass, &
     224             :                 level0_iDomain%nrows, level0_iDomain%ncols, level0_iDomain%xllcorner, &
     225          22 :                 level0_iDomain%yllcorner, level0_iDomain%cellsize, data_i4_2d, mask_2d)
     226             :         ! put global nodata value into array (probably not all grid cells have values)
     227     2612684 :         data_i4_2d = merge(data_i4_2d, nodata_i4, mask_2d)
     228          22 :         call paste(dataMatrix_i4, pack(data_i4_2d, level0_iDomain%mask), nodata_i4)
     229          44 :         deallocate(data_i4_2d)
     230             :       end do
     231          22 :       call append(L0_soilId, dataMatrix_i4)
     232          22 :       deallocate(dataMatrix_i4)
     233             : 
     234             :       ! read geoUnit
     235          22 :       fName = trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(file_hydrogeoclass))
     236             :       ! reading and transposing
     237             :       call read_spatial_data_ascii(trim(fName), uhydrogeoclass, &
     238             :               level0_iDomain%nrows, level0_iDomain%ncols, level0_iDomain%xllcorner, &
     239          22 :               level0_iDomain%yllcorner, level0_iDomain%cellsize, data_i4_2d, mask_2d)
     240             :       ! put global nodata value into array (probably not all grid cells have values)
     241     2612684 :       data_i4_2d = merge(data_i4_2d, nodata_i4, mask_2d)
     242          22 :       call append(L0_geoUnit, pack(data_i4_2d, level0_iDomain%mask))
     243          22 :       deallocate(data_i4_2d, mask_2d)
     244             : 
     245             :       ! LAI values
     246          22 :       call message('      Reading LAI ...')
     247             :       select case (timeStep_LAI_input)
     248             :       case(1) ! long term mean monthly gridded fields
     249           0 :         call prepare_gridded_mean_monthly_LAI_data(iDomain, level0_iDomain%nrows, level0_iDomain%ncols, level0_iDomain%mask)
     250             : 
     251             :       case(0) ! long term mean monthly values per class with LUT
     252             :         ! only set if not yet allocated (e.g. domain 1)
     253          21 :         if (.not. allocated(LAIBoundaries)) then
     254          12 :           nLAI = int(YearMonths, i4)
     255          12 :           allocate(LAIBoundaries(nLAI+1))
     256         492 :           LAIBoundaries = [(iMon, iMon=1, nLAI+1)]
     257             :         end if
     258             : 
     259          21 :         fName = trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(file_laiclass))
     260             :         ! reading and transposing
     261             :         call read_spatial_data_ascii(trim(fName), ulaiclass, &
     262             :                 level0_iDomain%nrows, level0_iDomain%ncols, level0_iDomain%xllcorner, &
     263          21 :                 level0_iDomain%yllcorner, level0_iDomain%cellsize, data_i4_2d, mask_2d)
     264             :         ! put global nodata value into array (probably not all grid cells have values)
     265     2487834 :         data_i4_2d = merge(data_i4_2d, nodata_i4, mask_2d)
     266     2487855 :         allocate(dummy_i4(count(level0_iDomain%mask)))
     267          21 :         dummy_i4 = pack(data_i4_2d, level0_iDomain%mask)
     268          21 :         deallocate(data_i4_2d, mask_2d)
     269             : 
     270             :         call check_consistency_lut_map(dummy_i4, LAIUnitList, file_laiclass)
     271             : 
     272     2487876 :         allocate(data_dp_2d(count(level0_iDomain%mask), nLAI))
     273         273 :         do iMon = 1, nLAI
     274             :           ! determine LAIs per month
     275        2793 :           do ll = 1, size(LAILUT, dim = 1)
     276   111171372 :             data_dp_2d(:, iMon) = merge(LAILUT(ll, iMon), data_dp_2d(:, iMon), dummy_i4(:) .EQ. LAIUnitList(ll))
     277             :           end do
     278             :         end do
     279          21 :         call append(L0_gridded_LAI, data_dp_2d(:, :))
     280          21 :         deallocate(dummy_i4, data_dp_2d)
     281             :         ! correction for 0 LAI values to avoid numerical instabilities
     282    24161913 :         L0_gridded_LAI(:, :) = merge(1.00E-10_dp, L0_gridded_LAI(:, :), L0_gridded_LAI(:, :) .LT. 1.00E-10_dp)
     283    24161934 :         L0_gridded_LAI(:, :) = merge(30.0_dp, L0_gridded_LAI(:, :), L0_gridded_LAI(:, :) .GT. 30.0_dp)
     284             :       case(-3 : -1) ! daily, monthly or yearly gridded fields (time-series)
     285             :         call prepare_gridded_daily_LAI_data(iDomain, level0_iDomain%nrows, level0_iDomain%ncols, level0_iDomain%mask, &
     286          22 :         LAIPer(iDomain))
     287             : 
     288             :       end select
     289             : 
     290             :       ! read lat lon coordinates of each domain
     291          22 :       call message('      Reading latitude/logitude ...')
     292          22 :       call read_latlon(iDomain, "lon_l0", "lat_l0", "level0", level0_iDomain)
     293             : 
     294          22 :       call timer_stop(itimer)
     295          57 :       call message('    in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
     296             : 
     297             :     end do domains
     298          13 :     itimer = 1
     299          13 :     call timer_stop(itimer)
     300          13 :     call message('  in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
     301             : 
     302             :     !Soil
     303             :     ! determine name od soil class definition file based on input option iFlag_soilDB
     304          13 :     if(iFlag_soilDB .eq. 0) then
     305          13 :       fName = file_soil_database
     306           0 :     else if(iFlag_soilDB .eq. 1) then
     307           0 :       fName = file_soil_database_1
     308             :     end if
     309             : 
     310             :     ! If you are getting segmentation fault with intel on large file (e.g. here on soils), then
     311             :     ! make sure to include following into your Marefile: INTEL_EXCLUDE  := mo_read_wrapper.f90
     312           0 :     call check_consistency_lut_map(reshape(L0_soilId, (/ size(L0_soilId, 1) * size(L0_soilId, 2) /)), &
     313          26 :             soilDB%id(:), fName)
     314             :     ! Geology
     315          13 :     call check_consistency_lut_map(L0_geoUnit, GeoUnitList, file_hydrogeoclass, dummy_i4)
     316             : 
     317             :     ! deactivate parameters of non existing geological classes in study domain for optimization
     318             :     ! loop over geological units in look up list
     319         143 :     do iVar = 1, size(GeoUnitList, 1)
     320             :       ! check if unit appears in geological map (dummy_i4 is unique number in L0_geoUnit)
     321         689 :       if (.not. ANY(dummy_i4 .EQ. GeoUnitList(iVar))) then
     322             :         ! deactivate optimization flag (dim=4 from global_parameters)
     323          39 :         global_parameters(processMatrix(9, 3) - processMatrix(9, 2) + iVar, 4) = 0
     324          39 :         call message('***WARNING: Geological unit ', trim(adjustl(num2str(GeoUnitList(iVar)))))
     325          39 :         call message('            is not appearing in study domain.')
     326             :       end if
     327             :     end do
     328             : 
     329          13 :     deallocate(dummy_i4) ! is allocated in subroutine check_consistency_lut_map - geology
     330             : 
     331             :     !----------------------------------------------------------------
     332             :     ! Correction for slope and aspect -- min value set above
     333             :     !----------------------------------------------------------------
     334             :     ! keep the colons (:) in the statements because of Intel's reallocation lhs problem
     335      972963 :     L0_slope(:) = merge(slope_minVal, L0_slope(:), (L0_slope(:)  .lt.  slope_minVal))
     336      972963 :     L0_asp(:) = merge(aspect_minVal, L0_asp(:), (L0_asp(:)    .lt. aspect_minVal))
     337             : 
     338          13 :   end subroutine read_data
     339             : 
     340             :   ! ------------------------------------------------------------------
     341             : 
     342             :   !    NAME
     343             :   !        check_consistency_lut_map
     344             : 
     345             :   !    PURPOSE
     346             :   !>       \brief Checks if classes in input maps appear in look up tables.
     347             : 
     348             :   !>       \details Determines wether a class appearing in the morphological input
     349             :   !>       is occuring in the respective look up table. mHM breaks if inconsistencies
     350             :   !>       are discovered.
     351             : 
     352             :   !    INTENT(IN)
     353             :   !>       \param[in] "integer(i4), dimension(:) :: data"        map of study domain
     354             :   !>       \param[in] "integer(i4), dimension(:) :: lookuptable" look up table corresponding to map
     355             :   !>       \param[in] "character(*) :: filename"                 name of the lut file - ERORR warn
     356             : 
     357             :   !    INTENT(OUT), OPTIONAL
     358             :   !>       \param[out] "integer(i4), dimension(:), optional :: unique_values" array of unique values in dataone
     359             : 
     360             :   !    HISTORY
     361             :   !>       \authors Matthias Zink
     362             : 
     363             :   !>       \date Nov 2016
     364             : 
     365             :   ! Modifications:
     366             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     367             : 
     368          94 :   subroutine check_consistency_lut_map(data, lookuptable, filename, unique_values)
     369             : 
     370          13 :     use mo_orderpack, only : unista
     371             :     use mo_string_utils, only : num2str
     372             : 
     373             :     implicit none
     374             : 
     375             :     ! map of study domain
     376             :     integer(i4), dimension(:), intent(in) :: data
     377             : 
     378             :     ! look up table corresponding to map
     379             :     integer(i4), dimension(:), intent(in) :: lookuptable
     380             : 
     381             :     ! name of the lut file - ERORR warn
     382             :     character(*), intent(in) :: filename
     383             : 
     384             :     ! array of unique values in dataone
     385             :     integer(i4), dimension(:), allocatable, intent(out), optional :: unique_values
     386             : 
     387             :     integer(i4) :: n_unique_elements
     388             : 
     389             :     integer(i4) :: ielement
     390             : 
     391          47 :     integer(i4), dimension(:), allocatable :: temp
     392             : 
     393             : 
     394         141 :     allocate(temp(size(data, 1)))
     395             :     ! copy L0_geoUnit because subroutine unista overwrites the nunit entries of the
     396             :     ! input array with the unique array values
     397     2872399 :     temp = data
     398             :     ! retrieve unique values of data
     399          47 :     call unista(temp, n_unique_elements)
     400             :     ! check if unit exists in look up table
     401         824 :     do ielement = 1, n_unique_elements
     402         777 :       if ( temp(ielement) == nodata_i4 ) call error_message( &
     403             :         '***ERROR: Class ', trim(adjustl(num2str(temp(ielement)))), &
     404             :         ' was searched in ', trim(adjustl(filename)), &
     405             :         ' which indicates a masking problem!' &
     406           0 :       )
     407      310520 :       if (.not. ANY(lookuptable .EQ. temp(ielement))) then
     408           0 :         call error_message('***ERROR: Class ', trim(adjustl(num2str(temp(ielement)))), ' is missing', raise=.false.)
     409           0 :         call error_message('          in input file ', trim(adjustl(filename)), ' ...')
     410             :       end if
     411             :     end do
     412             : 
     413             :     ! pass unique values if optional argument unique_values is given
     414          47 :     if (present(unique_values)) then
     415          39 :       allocate(unique_values(n_unique_elements))
     416         117 :       unique_values(:) = temp(1 : n_unique_elements)
     417             :     end if
     418          47 :     deallocate(temp)
     419             : 
     420          94 :   end subroutine check_consistency_lut_map
     421             : 
     422             : 
     423             : END MODULE mo_read_wrapper

Generated by: LCOV version 1.16