LCOV - code coverage report
Current view: top level - mRM - mo_mrm_read_data.f90 (source / functions) Hit Total Coverage
Test: mHM coverage Lines: 96 143 67.1 %
Date: 2024-04-15 17:48:09 Functions: 3 5 60.0 %

          Line data    Source code
       1             : !> \file mo_mrm_read_data.f90
       2             : !> \brief \copybrief mo_mrm_read_data
       3             : !> \details \copydetails mo_mrm_read_data
       4             : 
       5             : !> \brief mRM reading routines
       6             : !> \details This module contains all routines to read mRM data from file.
       7             : !> \authors Stephan Thober
       8             : !> \date Aug 2015
       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_mrm
      12             : module mo_mrm_read_data
      13             :   use mo_kind, only : i4, dp
      14             :   use mo_message, only: message, error_message
      15             : 
      16             :   implicit none
      17             : 
      18             :   public :: mrm_read_L0_data
      19             :   public :: mrm_read_discharge
      20             :   public :: mrm_read_total_runoff
      21             :   public :: mrm_read_bankfull_runoff
      22             :   private
      23             : contains
      24             :   ! ------------------------------------------------------------------
      25             : 
      26             :   !    NAME
      27             :   !        mrm_read_L0_data
      28             : 
      29             :   !    PURPOSE
      30             :   !>       \brief read L0 data from file
      31             : 
      32             :   !>       \details With the exception of L0_mask, L0_elev, and L0_LCover, all
      33             :   !>       L0 variables are read from file. The former three are only read if they
      34             :   !>       are not provided as variables.
      35             : 
      36             :   !    INTENT(IN)
      37             :   !>       \param[in] "logical :: do_reinit"
      38             :   !>       \param[in] "logical :: do_readlatlon"
      39             :   !>       \param[in] "logical :: do_readlcover"
      40             : 
      41             :   !    HISTORY
      42             :   !>       \authors Juliane Mai, Matthias Zink, and Stephan Thober
      43             : 
      44             :   !>       \date Aug 2015
      45             : 
      46             :   ! Modifications:
      47             :   ! Stephan Thober Sep 2015 - added L0_mask, L0_elev, and L0_LCover
      48             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
      49             :   ! Stephan Thober Jun 2018 - including varying celerity functionality
      50             : 
      51          12 :   subroutine mrm_read_L0_data(do_reinit, do_readlatlon, do_readlcover)
      52             : 
      53             :     use mo_append, only : append
      54             :     use mo_common_constants, only : nodata_i4
      55             :     use mo_common_read_data, only : read_dem, read_lcover
      56             :     use mo_common_types, only: Grid
      57             :     use mo_common_variables, only : L0_LCover, dirMorpho, level0, domainMeta, processMatrix
      58             :     use mo_mpr_file, only: file_slope, uslope
      59             :     use mo_mrm_file, only : file_facc, file_fdir, &
      60             :                             file_gaugeloc, ufacc, ufdir, ugaugeloc
      61             :     use mo_mrm_global_variables, only : L0_InflowGaugeLoc, L0_fAcc, L0_fDir, L0_gaugeLoc, domain_mrm
      62             :     use mo_read_latlon, only : read_latlon
      63             :     use mo_read_spatial_data, only : read_spatial_data_ascii
      64             :     use mo_string_utils, only : num2str
      65             : 
      66             :     implicit none
      67             : 
      68             :     logical, intent(in) :: do_reinit
      69             : 
      70             :     logical, intent(in) :: do_readlatlon
      71             : 
      72             :     logical, intent(in) :: do_readlcover
      73             : 
      74             :     integer(i4) ::domainID, iDomain
      75             : 
      76             :     integer(i4) :: iVar
      77             : 
      78             :     integer(i4) :: iGauge
      79             : 
      80             :     character(256) :: fname
      81             : 
      82             :     integer(i4) :: nunit
      83             : 
      84          12 :     integer(i4), dimension(:, :), allocatable :: data_i4_2d
      85             : 
      86          12 :     integer(i4), dimension(:, :), allocatable :: dataMatrix_i4
      87             : 
      88          12 :     logical, dimension(:, :), allocatable :: mask_2d
      89             : 
      90          12 :     logical, dimension(:, :), allocatable :: mask_global
      91             : 
      92             :     type(Grid), pointer :: level0_iDomain => null()
      93             : 
      94             : 
      95             :     ! ************************************************
      96             :     ! READ SPATIAL DATA FOR EACH DOMAIN
      97             :     ! ************************************************
      98          12 :     if (do_reinit) then
      99           0 :       call read_dem()
     100             :     end if
     101             : 
     102          12 :     if (do_readlcover .and. processMatrix(8, 1) .eq. 1) then
     103           9 :       call read_lcover()
     104           3 :     else if (do_readlcover .and. ((processMatrix(8, 1) .eq. 2) .or. (processMatrix(8, 1) .eq. 3))) then
     105           0 :       allocate(dataMatrix_i4(count(mask_global), 1))
     106           0 :       dataMatrix_i4 = nodata_i4
     107           0 :       call append(L0_LCover, dataMatrix_i4)
     108             :       ! free memory
     109           0 :       deallocate(dataMatrix_i4)
     110             :     end if
     111             : 
     112          36 :     do iDomain = 1, domainMeta%nDomains
     113          24 :       domainID = domainMeta%indices(iDomain)
     114             : 
     115          24 :       level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
     116             : 
     117             :       ! ToDo: check if change is correct
     118             :       ! check whether L0 data is shared
     119          24 :       if (domainMeta%L0DataFrom(iDomain) < iDomain) then
     120             :         !
     121             :         call message('      Using data of domain ', &
     122           3 :                 trim(adjustl(num2str(domainMeta%indices(domainMeta%L0DataFrom(iDomain))))), ' for domain: ',&
     123           3 :                 trim(adjustl(num2str(domainID))), '...')
     124           3 :         cycle
     125             :         !
     126             :       end if
     127             :       !
     128          21 :       call message('      Reading data for domain: ', trim(adjustl(num2str(domainID))), ' ...')
     129             : 
     130          25 :       if (do_readlatlon) then
     131             :         ! read lat lon coordinates of each domain
     132           4 :         call read_latlon(iDomain, "lon_l0", "lat_l0", "level0", level0_iDomain)
     133             :       end if
     134             : 
     135             :       ! read fAcc, fDir, gaugeLoc
     136         117 :       do iVar = 1, 4
     137          21 :         select case (iVar)
     138             :         case(1) ! flow accumulation
     139          21 :           fName = trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(file_facc))
     140          21 :           nunit = ufacc
     141             :         case(2) ! flow direction
     142          21 :           fName = trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(file_fdir))
     143          21 :           nunit = ufdir
     144             :         case(3) ! location of gauging stations
     145          21 :           fName = trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(file_gaugeloc))
     146          21 :           nunit = ugaugeloc
     147             :        case(4)
     148          21 :           fName = trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(file_slope))
     149         105 :           nunit = uslope
     150             :        end select
     151             : 
     152          84 :        if (iVar .le. 2) then
     153             :           !
     154             :           ! reading and transposing
     155             :           call read_spatial_data_ascii(trim(fName), nunit, &
     156             :                level0_iDomain%nrows, level0_iDomain%ncols, level0_iDomain%xllcorner, &
     157          42 :                level0_iDomain%yllcorner, level0_iDomain%cellsize, data_i4_2d, mask_2d)
     158             : 
     159             :           ! put global nodata value into array (probably not all grid cells have values)
     160     4975668 :           data_i4_2d = merge(data_i4_2d, nodata_i4, mask_2d)
     161             :        end if
     162             : 
     163          84 :        if (iVar .eq. 3) then
     164          21 :           if (domain_mrm(iDomain)%nGauges .ge. 1_i4) then
     165             :              !
     166             :              ! reading and transposing
     167             :              call read_spatial_data_ascii(trim(fName), nunit, &
     168             :                   level0_iDomain%nrows, level0_iDomain%ncols, level0_iDomain%xllcorner, &
     169          15 :                   level0_iDomain%yllcorner, level0_iDomain%cellsize, data_i4_2d, mask_2d)
     170             : 
     171             :              ! put global nodata value into array (probably not all grid cells have values)
     172     1805742 :              data_i4_2d = merge(data_i4_2d, nodata_i4, mask_2d)
     173             : 
     174             :           else
     175             :              ! set to nodata, but not ommit because data association between arrays and domains might break
     176      682092 :              data_i4_2d = merge(nodata_i4, nodata_i4, level0_iDomain%mask)
     177             :           end if
     178             :        end if
     179             : 
     180             :         ! put data into global L0 variable
     181             :         select case (iVar)
     182             :         case(1) ! flow accumulation
     183          21 :           call append(L0_fAcc, pack(data_i4_2d, level0_iDomain%mask))
     184             :         case(2) ! flow direction
     185             :           ! rotate flow direction and any other variable with directions
     186             :           ! NOTE: ONLY when ASCII files are read
     187          21 :           call rotate_fdir_variable(data_i4_2d)
     188             :           ! append
     189          21 :           call append(L0_fDir, pack(data_i4_2d, level0_iDomain%mask))
     190             :         case(3) ! location of evaluation and inflow gauging stations
     191             :           ! evaluation gauges
     192             :           ! Input data check
     193          37 :           do iGauge = 1, domain_mrm(iDomain)%nGauges
     194             :             ! If gaugeId is found in gauging location file?
     195      209111 :             if (.not. any(data_i4_2d .EQ. domain_mrm(iDomain)%gaugeIdList(iGauge))) then
     196             :               call error_message('***ERROR: Gauge ID "', trim(adjustl(num2str(domain_mrm(iDomain)%gaugeIdList(iGauge)))), &
     197           0 :                       '" not found in ', raise=.false.)
     198             :               call error_message('          Gauge location input file: ', &
     199           0 :                       trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(file_gaugeloc)))
     200             :             end if
     201             :           end do
     202             : 
     203          21 :           call append(L0_gaugeLoc, pack(data_i4_2d, level0_iDomain%mask))
     204             : 
     205             :           ! inflow gauges
     206             :           ! if no inflow gauge for this subdomain exists still matirx with dim of subdomain has to be paded
     207          21 :           if (domain_mrm(iDomain)%nInflowGauges .GT. 0_i4) then
     208             :             ! Input data check
     209           2 :             do iGauge = 1, domain_mrm(iDomain)%nInflowGauges
     210             :               ! If InflowGaugeId is found in gauging location file?
     211       55318 :               if (.not. any(data_i4_2d .EQ. domain_mrm(iDomain)%InflowGaugeIdList(iGauge))) then
     212             :                 call error_message('***ERROR: Inflow Gauge ID "', &
     213             :                         trim(adjustl(num2str(domain_mrm(iDomain)%InflowGaugeIdList(iGauge)))), &
     214           0 :                         '" not found in ', raise=.false.)
     215             :                 call error_message('          Gauge location input file: ', &
     216           0 :                         trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(file_gaugeloc)))
     217             :               end if
     218             :             end do
     219             :           end if
     220             : 
     221         105 :           call append(L0_InflowGaugeLoc, pack(data_i4_2d, level0_iDomain%mask))
     222             : 
     223             :         end select
     224             :         !
     225             :         ! deallocate arrays
     226          84 :         if (allocated(data_i4_2d)) deallocate(data_i4_2d)
     227         105 :         if (allocated(mask_2d)) deallocate(mask_2d)
     228             :         !
     229             :       end do
     230             :     end do
     231             : 
     232          12 :   end subroutine mrm_read_L0_data
     233             :   ! ---------------------------------------------------------------------------
     234             : 
     235             :   !    NAME
     236             :   !        mrm_read_discharge
     237             : 
     238             :   !    PURPOSE
     239             :   !>       \brief Read discharge timeseries from file
     240             : 
     241             :   !>       \details Read Observed discharge at the outlet of a catchment
     242             :   !>       and at the inflow of a catchment. Allocate global runoff
     243             :   !>       variable that contains the simulated runoff after the simulation.
     244             : 
     245             :   !    HISTORY
     246             :   !>       \authors Matthias Zink & Stephan Thober
     247             : 
     248             :   !>       \date Aug 2015
     249             : 
     250             :   ! Modifications:
     251             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     252             : 
     253          13 :   subroutine mrm_read_discharge
     254             : 
     255          12 :     use mo_append, only : paste
     256             :     use mo_common_constants, only : nodata_dp
     257             :     use mo_common_mHM_mRM_variables, only : evalPer, nTstepDay, opti_function, optimize, simPer
     258             :     use mo_common_variables, only : domainMeta
     259             :     use mo_mrm_file, only : udischarge
     260             :     use mo_mrm_global_variables, only : InflowGauge, gauge, mRM_runoff, nGaugesLocal, &
     261             :                                         nInflowGaugesTotal, nMeasPerDay, &
     262             :                                         riv_temp_pcs
     263             :     use mo_read_timeseries, only : read_timeseries
     264             :     use mo_string_utils, only : num2str
     265             : 
     266             :     implicit none
     267             : 
     268             :     integer(i4) :: iGauge
     269             : 
     270             :     integer(i4) :: iDomain
     271             : 
     272             :     integer(i4) :: maxTimeSteps
     273             : 
     274             :     ! file name of file to read
     275             :     character(256) :: fName
     276             : 
     277             :     integer(i4), dimension(3) :: start_tmp, end_tmp
     278             : 
     279          13 :     real(dp), dimension(:), allocatable :: data_dp_1d
     280             : 
     281          13 :     logical, dimension(:), allocatable :: mask_1d
     282             : 
     283             : 
     284             :     !----------------------------------------------------------
     285             :     ! INITIALIZE RUNOFF
     286             :     !----------------------------------------------------------
     287          51 :     maxTimeSteps = maxval(simPer(1 : domainMeta%nDomains)%julEnd - simPer(1 : domainMeta%nDomains)%julStart + 1) * nTstepDay
     288          52 :     allocate(mRM_runoff(maxTimeSteps, nGaugesLocal))
     289      358584 :     mRM_runoff = nodata_dp
     290             : 
     291             :     ! READ GAUGE DATA
     292          35 :     do iGauge = 1, nGaugesLocal
     293             :       ! get domain id
     294          22 :       iDomain = gauge%domainId(iGauge)
     295             :       ! get start and end dates
     296          88 :       start_tmp = (/evalPer(iDomain)%yStart, evalPer(iDomain)%mStart, evalPer(iDomain)%dStart/)
     297          88 :       end_tmp = (/evalPer(iDomain)%yEnd, evalPer(iDomain)%mEnd, evalPer(iDomain)%dEnd  /)
     298             :       ! evaluation gauge
     299          22 :       fName = trim(adjustl(gauge%fname(iGauge)))
     300             :       call read_timeseries(trim(fName), udischarge, &
     301             :               start_tmp, end_tmp, optimize, opti_function, &
     302          22 :               data_dp_1d, mask = mask_1d, nMeasPerDay = nMeasPerDay)
     303        9540 :       data_dp_1d = merge(data_dp_1d, nodata_dp, mask_1d)
     304          22 :       call paste(gauge%Q, data_dp_1d, nodata_dp)
     305          35 :       deallocate (data_dp_1d)
     306             :       ! TODO-RIV-TEMP: read temperature at gauge
     307             :     end do
     308             :     !
     309             :     ! inflow gauge
     310             :     !
     311             :     ! in mhm call InflowGauge%Q has to be initialized -- dummy allocation with period of domain 1 and initialization
     312          13 :     if (nInflowGaugesTotal .EQ. 0) then
     313          58 :       allocate(data_dp_1d(maxval(simPer(:)%julEnd - simPer(:)%julStart + 1)))
     314        8043 :       data_dp_1d = nodata_dp
     315          12 :       call paste(InflowGauge%Q, data_dp_1d, nodata_dp)
     316             :     else
     317           3 :       do iGauge = 1, nInflowGaugesTotal
     318             :         ! get domain id
     319           2 :         iDomain = InflowGauge%domainId(iGauge)
     320             :         ! get start and end dates
     321           8 :         start_tmp = (/simPer(iDomain)%yStart, simPer(iDomain)%mStart, simPer(iDomain)%dStart/)
     322           8 :         end_tmp = (/simPer(iDomain)%yEnd, simPer(iDomain)%mEnd, simPer(iDomain)%dEnd  /)
     323             :         ! inflow gauge
     324           2 :         fName = trim(adjustl(InflowGauge%fname(iGauge)))
     325             :         call read_timeseries(trim(fName), udischarge, &
     326             :                 start_tmp, end_tmp, optimize, opti_function, &
     327           2 :                 data_dp_1d, mask = mask_1d, nMeasPerDay = nMeasPerDay)
     328        1094 :         if (.NOT. (all(mask_1d))) then
     329           0 :           call error_message('***ERROR: Nodata values in inflow gauge time series. File: ', trim(fName), raise=.false.)
     330           0 :           call error_message('          During simulation period from ', num2str(simPer(iDomain)%yStart) &
     331           0 :                   , ' to ', num2str(simPer(iDomain)%yEnd))
     332             :         end if
     333        1096 :         data_dp_1d = merge(data_dp_1d, nodata_dp, mask_1d)
     334           2 :         call paste(InflowGauge%Q, data_dp_1d, nodata_dp)
     335           3 :         deallocate (data_dp_1d)
     336             :       end do
     337             :     end if
     338             : 
     339          13 :   end subroutine mrm_read_discharge
     340             : 
     341             :   ! ---------------------------------------------------------------------------
     342             : 
     343             :   !    NAME
     344             :   !        mrm_read_total_runoff
     345             : 
     346             :   !    PURPOSE
     347             :   !>       \brief read simulated runoff that is to be routed
     348             : 
     349             :   !>       \details read spatio-temporal field of total runoff that has been
     350             :   !>       simulated by a hydrologic model or land surface model. This
     351             :   !>       total runoff will then be aggregated to the level 11 resolution
     352             :   !>       and then routed through the stream network.
     353             : 
     354             :   !    INTENT(IN)
     355             :   !>       \param[in] "integer(i4) :: iDomain" domain id
     356             : 
     357             :   !    HISTORY
     358             :   !>       \authors Stephan Thober
     359             : 
     360             :   !>       \date Sep 2015
     361             : 
     362             :   ! Modifications:
     363             :   ! Stephan Thober  Feb 2016 - refactored deallocate statements
     364             :   ! Stephan Thober  Sep 2016 - added ALMA convention
     365             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     366             : 
     367           0 :   subroutine mrm_read_total_runoff(iDomain)
     368             : 
     369          13 :     use mo_append, only : append
     370             :     use mo_constants, only : HourSecs
     371             :     use mo_common_constants, only : nodata_dp
     372             :     use mo_common_mHM_mRM_variables, only : simPer, timestep
     373             :     use mo_common_variables, only : ALMA_convention, level1
     374             :     use mo_mrm_global_variables, only : L1_total_runoff_in, dirTotalRunoff, filenameTotalRunoff, &
     375             :                                         varnameTotalRunoff
     376             :     use mo_read_nc, only : read_nc
     377             : 
     378             :     implicit none
     379             : 
     380             :     ! domain id
     381             :     integer(i4), intent(in) :: iDomain
     382             : 
     383             :     integer(i4) :: tt
     384             : 
     385             :     integer(i4) :: nTimeSteps
     386             : 
     387             :     ! tell nc file to read daily or hourly values
     388             :     integer(i4) :: nctimestep
     389             : 
     390             :     ! read data from file
     391           0 :     real(dp), dimension(:, :, :), allocatable :: L1_data
     392             : 
     393           0 :     real(dp), dimension(:, :), allocatable :: L1_data_packed
     394             : 
     395             : 
     396           0 :     if (timestep .eq. 1) nctimestep = -4 ! hourly input
     397           0 :     if (timestep .eq. 24) nctimestep = -1 ! daily input
     398           0 :     call read_nc(trim(dirTotalRunoff(iDomain)), level1(iDomain)%nrows, level1(iDomain)%ncols, &
     399           0 :             varnameTotalRunoff, level1(iDomain)%mask, L1_data, target_period = simPer(iDomain), &
     400           0 :             nctimestep = nctimestep, filename = filenameTotalRunoff)
     401             :     ! pack variables
     402           0 :     nTimeSteps = size(L1_data, 3)
     403           0 :     allocate(L1_data_packed(level1(iDomain)%nCells, nTimeSteps))
     404           0 :     do tt = 1, nTimeSteps
     405           0 :       L1_data_packed(:, tt) = pack(L1_data(:, :, tt), mask = level1(iDomain)%mask)
     406             :     end do
     407             :     ! free space immediately
     408           0 :     deallocate(L1_data)
     409             : 
     410             :     ! convert if ALMA conventions have been given
     411           0 :     if (ALMA_convention) then
     412             :       ! convert from kg m-2 s-1 to mm TS-1
     413             :       ! 1 kg m-2 -> 1 mm depth
     414             :       ! multiply with time to obtain per timestep
     415           0 :       L1_data_packed = L1_data_packed * timestep * HourSecs
     416             :     end if
     417             : 
     418             :     ! append
     419           0 :     call append(L1_total_runoff_in, L1_data_packed(:, :), nodata_dp)
     420             : 
     421             :     !free space
     422           0 :     deallocate(L1_data_packed)
     423             : 
     424           0 :   end subroutine mrm_read_total_runoff
     425             : 
     426           0 :   subroutine mrm_read_bankfull_runoff(iDomain)
     427             :   ! ---------------------------------------------------------------------------
     428             : 
     429             :   !      NAME
     430             :   !          mrm_read_bankfull_runoff
     431             : 
     432             :   !>         \brief reads the bankfull runoff for approximating the channel widths
     433             : 
     434             :   !>         \details reads the bankfull runoff, which can be calculated with
     435             :   !>         the script in mhm/post_proc/bankfull_discharge.py
     436             : 
     437             :   !     INTENT(IN)
     438             :   !>        \param[in] "integer(i4)               :: iDomain"  domain id
     439             : 
     440             :   !     INTENT(INOUT)
     441             :   !         None
     442             : 
     443             :   !     INTENT(OUT)
     444             :   !         None
     445             : 
     446             :   !     INTENT(IN), OPTIONAL
     447             :   !         None
     448             : 
     449             :   !     INTENT(INOUT), OPTIONAL
     450             :   !         None
     451             : 
     452             :   !     INTENT(OUT), OPTIONAL
     453             :   !         None
     454             : 
     455             :   !     RETURN
     456             :   !         None
     457             : 
     458             :   !     RESTRICTIONS
     459             :   !>        \note The file read in must contain a double precision float variable with the name
     460             :   !>        "Q_bkfl".
     461             : 
     462             :   !     EXAMPLE
     463             :   !         None
     464             : 
     465             :   !     LITERATURE
     466             :   !         None
     467             : 
     468             :   !     HISTORY
     469             :   !         \author Lennart Schueler
     470             :   !         \date    May 2018
     471             : 
     472           0 :     use mo_mrm_global_variables, only: level11
     473             :     use mo_read_nc, only: read_const_nc
     474             :     use mo_mrm_global_variables, only: &
     475             :          dirBankfullRunoff, &   ! directory of bankfull_runoff file for each domain
     476             :          L11_bankfull_runoff_in ! bankfull runoff at L1
     477             : 
     478             :     implicit none
     479             : 
     480             :     ! input variables
     481             :     integer(i4), intent(in) :: iDomain
     482             : 
     483             :     ! local variables
     484           0 :     real(dp), dimension(:,:), allocatable :: L11_data ! read data from file
     485           0 :     real(dp), dimension(:), allocatable :: L11_data_packed
     486             : 
     487           0 :     call read_const_nc(trim(dirBankfullRunoff(iDomain)), &
     488           0 :                                level11(iDomain)%nrows, &
     489             :                                level11(iDomain)%ncols, &
     490           0 :                                "Q_bkfl", L11_data)
     491             : 
     492           0 :     allocate(L11_data_packed(level11(iDomain)%nCells))
     493           0 :     L11_data_packed(:) = pack(L11_data(:,:), mask=level11(iDomain)%mask)
     494             : 
     495             :     ! append
     496           0 :     if (allocated(L11_bankfull_runoff_in)) then
     497           0 :         L11_bankfull_runoff_in = [L11_bankfull_runoff_in, L11_data_packed]
     498             :     else
     499           0 :         allocate(L11_bankfull_runoff_in(size(L11_data_packed)))
     500           0 :         L11_bankfull_runoff_in = L11_data_packed
     501             :     end if
     502             : 
     503           0 :     deallocate(L11_data)
     504           0 :     deallocate(L11_data_packed)
     505             : 
     506           0 :   end subroutine mrm_read_bankfull_runoff
     507             : 
     508             :   !    NAME
     509             :   !        rotate_fdir_variable
     510             : 
     511             :   !    PURPOSE
     512             :   !>       \brief TODO: add description
     513             : 
     514             :   !>       \details TODO: add description
     515             : 
     516             :   !    INTENT(INOUT)
     517             :   !>       \param[inout] "integer(i4), dimension(:, :) :: x"
     518             : 
     519             :   !    HISTORY
     520             :   !>       \authors L. Samaniego & R. Kumar
     521             : 
     522             :   !>       \date Jun 2018
     523             : 
     524             :   ! Modifications:
     525             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     526             : 
     527          21 :   subroutine rotate_fdir_variable(x)
     528             : 
     529           0 :     use mo_common_constants, only : nodata_i4
     530             : 
     531             :     implicit none
     532             : 
     533             :     integer(i4), dimension(:, :), intent(INOUT) :: x
     534             : 
     535             : 
     536             :     ! 0  -1   0 |
     537             :     integer(i4) :: i, j
     538             : 
     539             : 
     540             :     !-------------------------------------------------------------------
     541             :     ! NOTE:
     542             :     !
     543             :     !     Since the DEM was transposed from (lat,lon) to (lon,lat), i.e.
     544             :     !     new DEM = transpose(DEM_old), then
     545             :     !
     546             :     !     the flow direction X (which was read) for every i, j cell needs
     547             :     !     to be rotated as follows
     548             :     !
     549             :     !                 X(i,j) = [R] * {uVector}
     550             :     !
     551             :     !     with
     552             :     !     {uVector} denoting the fDir_old (read) in vector form, and
     553             :     !               e.g. dir 8 => {-1, -1, 0 }
     554             :     !     [R]       denting a full rotation matrix to transform the flow
     555             :     !               direction into the new coordinate system (lon,lat).
     556             :     !
     557             :     !     [R] = [rx][rz]
     558             :     !
     559             :     !     with
     560             :     !           |      1       0      0 |
     561             :     !     [rx] =|      0   cos T  sin T | = elemental rotation along x axis
     562             :     !           |      0  -sin T  cos T |
     563             :     !
     564             :     !           |  cos F   sin F      0 |
     565             :     !     [rz] =| -sin F   cos F      0 | = elemental rotation along z axis
     566             :     !           |      0       0      1 |
     567             :     !
     568             :     !     and T = pi, F = - pi/2
     569             :     !     thus
     570             :     !          |  0  -1   0 |
     571             :     !     [R] =| -1   0   0 |
     572             :     !          |  0   0  -1 |
     573             :     !     making all 8 directions the following transformation were
     574             :     !     obtained.
     575             :     !-------------------------------------------------------------------
     576        5973 :     do i = 1, size(x, 1)
     577     2485077 :       do j = 1, size(x, 2)
     578     2479104 :         if (x(i, j)  .eq. nodata_i4) cycle
     579        5952 :         select case (x(i, j))
     580             :         case(1)
     581      122693 :           x(i, j) = 4
     582             :         case(2)
     583       73066 :           x(i, j) = 2
     584             :         case(4)
     585      123704 :           x(i, j) = 1
     586             :         case(8)
     587       84531 :           x(i, j) = 128
     588             :         case(16)
     589      159247 :           x(i, j) = 64
     590             :         case(32)
     591      102046 :           x(i, j) = 32
     592             :         case(64)
     593      163544 :           x(i, j) = 16
     594             :         case(128)
     595      926405 :           x(i, j) = 8
     596             :         end select
     597             :       end do
     598             :     end do
     599             : 
     600          21 :   end subroutine rotate_fdir_variable
     601             : 
     602             : end module mo_mrm_read_data

Generated by: LCOV version 1.16