LCOV - code coverage report
Current view: top level - mHM - mo_mhm_interface_run.f90 (source / functions) Hit Total Coverage
Test: mHM coverage Lines: 206 338 60.9 %
Date: 2024-04-30 08:53:32 Functions: 9 9 100.0 %

          Line data    Source code
       1             : !> \file    mo_mhm_interface_run.f90
       2             : !> \brief   \copybrief mo_mhm_interface_run
       3             : !> \details \copydetails mo_mhm_interface_run
       4             : 
       5             : !> \brief   Module providing interfaces for running preconfigured mHM.
       6             : !> \details Interfaces to control the mHM run from outside (prepare domain, do timestep, ...).
       7             : !> \authors Sebastian Mueller, Matthias Kelbling
       8             : !> \version 0.1
       9             : !> \date    Jan 2022
      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_mhm
      13             : module mo_mhm_interface_run
      14             : 
      15             :   ! forces
      16             :   use mo_kind, only: i4, dp
      17             :   use mo_message, only: message, error_message
      18             :   use mo_string_utils, only : num2str
      19             :   ! mhm
      20             :   use mo_common_run_variables, only : run_cfg
      21             :   use mo_optimization_types, only : optidata_sim
      22             :   use mo_common_datetime_type, only : datetimeinfo
      23             :   use mo_common_mHM_mRM_variables, only : &
      24             :     resolutionRouting, &
      25             :     LCyearId, &
      26             :     mhmFileRestartIn, &
      27             :     mrmFileRestartIn, &
      28             :     nTstepDay, &
      29             :     optimize, &
      30             :     read_restart, &
      31             :     simPer, &
      32             :     timeStep, &
      33             :     warmingDays, &
      34             :     c2TSTu, &
      35             :     restart_reset_fluxes_states
      36             :   use mo_common_variables, only : &
      37             :     global_parameters, &
      38             :     level0, &
      39             :     level1, &
      40             :     domainMeta, &
      41             :     processMatrix
      42             :   use mo_global_variables, only : &
      43             :     meteo_handler, &
      44             :     L1_Throughfall, &
      45             :     L1_aETCanopy, &
      46             :     L1_aETSealed, &
      47             :     L1_aETSoil, &
      48             :     L1_baseflow, &
      49             :     L1_fastRunoff, &
      50             :     L1_infilSoil, &
      51             :     L1_inter, &
      52             :     L1_melt, &
      53             :     L1_neutrons, &
      54             :     L1_percol, &
      55             :     L1_pet_calc, &
      56             :     L1_temp_calc, &
      57             :     L1_prec_calc, &
      58             :     L1_preEffect, &
      59             :     L1_rain, &
      60             :     L1_runoffSeal, &
      61             :     L1_satSTW, &
      62             :     L1_sealSTW, &
      63             :     L1_slowRunoff, &
      64             :     L1_snow, &
      65             :     L1_snowPack, &
      66             :     L1_soilMoist, &
      67             :     L1_total_runoff, &
      68             :     L1_unsatSTW, &
      69             :     L1_twsaObs, &
      70             :     L1_etObs, &
      71             :     L1_smObs, &
      72             :     L1_neutronsObs, &
      73             :     evap_coeff, &
      74             :     nSoilHorizons_sm_input, &
      75             :     neutron_integral_AFast, &
      76             :     outputFlxState, &
      77             :     timeStep_model_outputs, &
      78             :     BFI_qBF_sum, &
      79             :     BFI_qT_sum
      80             :   use mo_init_states, only : variables_default_init, fluxes_states_default_init
      81             :   use mo_julian, only : caldat, julday
      82             :   use mo_string_utils, only : num2str
      83             :   use mo_mhm, only : mhm
      84             :   use mo_restart, only : read_restart_states
      85             :   use mo_write_fluxes_states, only : mHM_updateDataset, mHM_OutputDataset
      86             :   use mo_mrm_write_fluxes_states, only : mRM_updateDataset, mRM_OutputDataset, GW_OutputDataset, GW_updateDataset
      87             :   use mo_constants, only : HourSecs
      88             :   use mo_common_variables, only : resolutionHydrology
      89             :   use mo_mrm_global_variables, only : &
      90             :     InflowGauge, &
      91             :     L11_C1, &
      92             :     L11_C2, &
      93             :     L11_L1_Id, &
      94             :     L11_TSrout, &
      95             :     L11_fromN, &
      96             :     L11_length, &
      97             :     L11_nLinkFracFPimp, &
      98             :     L11_nOutlets, &
      99             :     L11_netPerm, &
     100             :     L11_qMod, &
     101             :     L11_qOUT, &
     102             :     L11_qTIN, &
     103             :     L11_qTR, &
     104             :     L11_slope, &
     105             :     L11_toN, &
     106             :     L1_L11_Id, &
     107             :     domain_mrm, &
     108             :     level11, &
     109             :     mRM_runoff, &
     110             :     outputFlxState_mrm, &
     111             :     timeStep_model_outputs_mrm, &
     112             :     gw_coupling, &
     113             :     L0_river_head_mon_sum, &
     114             :     riv_temp_pcs
     115             :   use mo_mpr_global_variables, only : &
     116             :     L1_HarSamCoeff, &
     117             :     L1_PrieTayAlpha, &
     118             :     L1_aeroResist, &
     119             :     L1_alpha, &
     120             :     L1_degDay, &
     121             :     L1_degDayInc, &
     122             :     L1_degDayMax, &
     123             :     L1_degDayNoPre, &
     124             :     L1_fAsp, &
     125             :     L1_fRoots, &
     126             :     L1_fSealed, &
     127             :     L1_jarvis_thresh_c1, &
     128             :     L1_kBaseFlow, &
     129             :     L1_kPerco, &
     130             :     L1_kSlowFlow, &
     131             :     L1_karstLoss, &
     132             :     L1_kfastFlow, &
     133             :     L1_maxInter, &
     134             :     L1_petLAIcorFactor, &
     135             :     L1_sealedThresh, &
     136             :     L1_soilMoistExp, &
     137             :     L1_soilMoistFC, &
     138             :     L1_soilMoistSat, &
     139             :     L1_surfResist, &
     140             :     L1_tempThresh, &
     141             :     L1_unsatThresh, &
     142             :     L1_wiltingPoint, &
     143             :     L1_No_Count, &
     144             :     L1_bulkDens, &
     145             :     L1_latticeWater, &
     146             :     L1_COSMICL3, &
     147             :     HorizonDepth_mHM, &
     148             :     nSoilHorizons_mHM
     149             :   use mo_mrm_init, only : variables_default_init_routing, fluxes_states_default_init_routing
     150             :   use mo_mrm_mpr, only : mrm_update_param
     151             :   use mo_mrm_restart, only : mrm_read_restart_states
     152             :   use mo_mrm_routing, only : mrm_routing
     153             :   use mo_utils, only : ge
     154             :   use mo_mrm_river_head, only: calc_river_head
     155             :   use mo_mpr_eval, only : mpr_eval
     156             : 
     157             : contains
     158             : 
     159             :   !> \brief prepare single run of mHM
     160         122 :   subroutine mhm_interface_run_prepare(parameterset, opti_domain_indices, runoff_present, BFI_present)
     161             :     implicit none
     162             :     !> a set of global parameter (gamma) to run mHM, DIMENSION [no. of global_Parameters]
     163             :     real(dp), dimension(:), optional, intent(in) :: parameterset
     164             :     !> selected domains for optimization
     165             :     integer(i4), dimension(:), optional, intent(in) :: opti_domain_indices
     166             :     !> whether runoff is present
     167             :     logical, optional, intent(in) :: runoff_present
     168             :     !> whether BFI is present
     169             :     logical, optional, intent(in) :: BFI_present
     170             : 
     171             :     integer(i4) :: i, iDomain, domainID
     172             : 
     173             :     ! run_cfg%output_runoff = .false. by default
     174          61 :     if (present(runoff_present)) run_cfg%output_runoff = runoff_present
     175             :     ! run_cfg%output_BFI = .false. by default
     176          61 :     if (present(BFI_present)) run_cfg%output_BFI = BFI_present
     177             : 
     178             :     ! store current parameter set
     179         183 :     allocate(run_cfg%parameterset(size(global_parameters, dim=1)))
     180          61 :     if (.not. present(parameterset) .and. optimize) then
     181           0 :       call error_message("mhm_interface_run_prepare: Can't optimize without parameter!")
     182          61 :     else if (.not. present(parameterset)) then
     183           0 :       run_cfg%parameterset = global_parameters(:, 3)
     184             :     else
     185        3283 :       run_cfg%parameterset = parameterset
     186             :     end if
     187             : 
     188             :     ! store domain indices for domain loop
     189          61 :     if (optimize .and. present(opti_domain_indices)) then
     190          28 :       run_cfg%nDomains = size(opti_domain_indices)
     191          84 :       allocate(run_cfg%domain_indices(run_cfg%nDomains))
     192          98 :       run_cfg%domain_indices = opti_domain_indices
     193             :     else
     194          33 :       run_cfg%nDomains = domainMeta%nDomains
     195          99 :       allocate(run_cfg%domain_indices(run_cfg%nDomains))
     196         186 :       run_cfg%domain_indices = [(i, i=1, run_cfg%nDomains)]
     197             :     end if
     198             : 
     199             :     !----------------------------------------------------------
     200             :     ! Check optionals and initialize
     201             :     !----------------------------------------------------------
     202          61 :     if (run_cfg%output_runoff) then
     203          50 :       do i = 1, run_cfg%nDomains
     204          25 :         iDomain = run_cfg%get_domain_index(i)
     205          25 :         domainID = domainMeta%indices(iDomain)
     206          50 :         if (.not. domainMeta%doRouting(iDomain)) then
     207             :           call error_message("***ERROR: runoff for domain", trim(num2str(domainID)),&
     208           0 :                         "can not be produced, since routing process is off in Process Matrix")
     209             :         end if
     210             :       end do
     211             :     end if
     212             : 
     213             :     ! prepare BFI calculation
     214          61 :     if (run_cfg%output_BFI) then
     215           0 :       allocate(BFI_qBF_sum(run_cfg%nDomains))
     216           0 :       allocate(BFI_qT_sum(run_cfg%nDomains))
     217           0 :       BFI_qBF_sum = 0.0_dp
     218           0 :       BFI_qT_sum = 0.0_dp
     219             :     end if
     220             : 
     221          61 :     if (read_restart) then
     222           2 :       do i = 1, run_cfg%nDomains
     223           1 :         iDomain = run_cfg%get_domain_index(i)
     224           1 :         domainID = domainMeta%indices(iDomain)
     225             :         ! this reads the eff. parameters and optionally the states and fluxes
     226           2 :         call read_restart_states(iDomain, domainID, mhmFileRestartIn(iDomain))
     227             :       end do
     228           1 :       if (restart_reset_fluxes_states) then
     229           0 :         call message('    Resetting mHM states and fluxes from restart files ...')
     230           0 :         call fluxes_states_default_init()
     231             :       end if
     232             :     else
     233          60 :       call variables_default_init()
     234          60 :       call mpr_eval(run_cfg%parameterset)
     235          60 :       if (processMatrix(8, 1) > 0) call variables_default_init_routing()
     236             :     end if
     237             : 
     238         305 :     allocate(run_cfg%L1_fNotSealed(size(L1_fSealed, 1), size(L1_fSealed, 2), size(L1_fSealed, 3)))
     239       15650 :     run_cfg%L1_fNotSealed = 1.0_dp - L1_fSealed
     240             : 
     241          61 :   end subroutine mhm_interface_run_prepare
     242             : 
     243             :   !> \brief get number of domains for looping
     244          61 :   subroutine mhm_interface_run_get_ndomains(ndomains)
     245             :     implicit none
     246             :     integer(i4), intent(inout) :: ndomains !< number of domains
     247          61 :     ndomains = run_cfg%nDomains
     248          61 :   end subroutine mhm_interface_run_get_ndomains
     249             : 
     250             :   !> \brief prepare single domain to run mHM on
     251         164 :   subroutine mhm_interface_run_prepare_domain(domain, etOptiSim, twsOptiSim, neutronsOptiSim, smOptiSim)
     252             :     implicit none
     253             :     !> domain loop counter
     254             :     integer(i4), intent(in), optional :: domain
     255             :     !> returns soil moisture time series for all grid cells (of multiple Domains concatenated),DIMENSION [nCells, nTimeSteps]
     256             :     type(optidata_sim), dimension(:), optional, intent(inout) :: smOptiSim
     257             :     !> dim1=ncells, dim2=time
     258             :     type(optidata_sim), dimension(:), optional, intent(inout) :: neutronsOptiSim
     259             :     !> returns evapotranspiration time series for all grid cells (of multiple Domains concatenated),DIMENSION [nCells, nTimeSteps]
     260             :     type(optidata_sim), dimension(:), optional, intent(inout) :: etOptiSim
     261             :     !> returns tws time series for all grid cells (of multiple Domains concatenated),DIMENSION [nCells, nTimeSteps]
     262             :     type(optidata_sim), dimension(:), optional, intent(inout) :: twsOptiSim
     263             : 
     264             :     integer(i4) :: iDomain, domainID
     265             : 
     266             :     ! get domain index
     267          82 :     run_cfg%selected_domain = 1_i4
     268          82 :     if (present(domain)) run_cfg%selected_domain = domain
     269          82 :     iDomain = run_cfg%get_domain_index(run_cfg%selected_domain)
     270          82 :     domainID = domainMeta%indices(iDomain)
     271             : 
     272             :     ! evapotranspiration optimization
     273          82 :     if (present(etOptiSim)) call etOptiSim(iDomain)%init(L1_etObs(iDomain))
     274             :     ! total water storage optimization
     275          82 :     if (present(twsOptiSim)) call twsOptiSim(iDomain)%init(L1_twsaObs(iDomain))
     276             :     ! neutrons optimization
     277          82 :     if (present(neutronsOptiSim)) call neutronsOptiSim(iDomain)%init(L1_neutronsObs(iDomain))
     278             :     ! sm optimization
     279          82 :     if (present(smOptiSim)) call smOptiSim(iDomain)%init(L1_smObs(iDomain))
     280             : 
     281             :     ! get Domain information
     282          82 :     run_cfg%nCells = level1(iDomain)%nCells
     283          82 :     run_cfg%mask1 => level1(iDomain)%mask
     284          82 :     run_cfg%s1 = level1(iDomain)%iStart
     285          82 :     run_cfg%e1 = level1(iDomain)%iEnd
     286             : 
     287          82 :     if (domainMeta%doRouting(iDomain)) then
     288             :       ! ----------------------------------------
     289             :       ! initialize factor between routing resolution and hydrologic model resolution
     290             :       ! ----------------------------------------
     291          46 :       run_cfg%tsRoutFactor = 1_i4
     292         138 :       allocate(run_cfg%InflowDischarge(size(InflowGauge%Q, dim = 2)))
     293          95 :       run_cfg%InflowDischarge = 0._dp
     294             : 
     295             :       ! read states from restart
     296          46 :       if (read_restart) then
     297           1 :         call mrm_read_restart_states(iDomain, domainID, mrmFileRestartIn(iDomain))
     298           1 :         if (restart_reset_fluxes_states) then
     299           0 :           call message('    Resetting mRM states and fluxes from restart files for domain ', num2str(iDomain), ' ...')
     300           0 :           call fluxes_states_default_init_routing(iDomain)
     301             :         end if
     302             :       end if
     303             :       ! get Domain information at L11 if routing is activated
     304          46 :       run_cfg%s11 = level11(iDomain)%iStart
     305          46 :       run_cfg%e11 = level11(iDomain)%iEnd
     306          46 :       run_cfg%mask11 => level11(iDomain)%mask
     307             : 
     308             :       ! initialize routing parameters (has to be called for routing option 2)
     309          46 :       if ((processMatrix(8, 1) .eq. 2) .or. (processMatrix(8, 1) .eq. 3)) &
     310             :           call mrm_update_param( &
     311          10 :             iDomain, run_cfg%parameterset(processMatrix(8, 3) - processMatrix(8, 2) + 1 : processMatrix(8, 3)))
     312             :       ! initialize variable for runoff for routing
     313         138 :       allocate(run_cfg%RunToRout(run_cfg%e1 - run_cfg%s1 + 1))
     314        2180 :       run_cfg%RunToRout = 0._dp
     315             : 
     316          46 :       if ( riv_temp_pcs%active ) then
     317             :         ! set indices for current L11 domain
     318           1 :         riv_temp_pcs%s11 = run_cfg%s11
     319           1 :         riv_temp_pcs%e11 = run_cfg%e11
     320             :         ! allocate current L1 lateral components
     321           1 :         call riv_temp_pcs%alloc_lateral(run_cfg%nCells)
     322             :       end if
     323             :     end if
     324             : 
     325             :     ! init datetime variable
     326          82 :     call run_cfg%domainDateTime%init(iDomain)
     327             : 
     328             :     ! init time-loop
     329          82 :     run_cfg%time_step = 0_i4
     330             : 
     331         143 :   end subroutine mhm_interface_run_prepare_domain
     332             : 
     333             :   !> \brief check if current time loop is finished
     334     2020272 :   subroutine mhm_interface_run_finished(time_loop_finished)
     335             :     implicit none
     336             :     logical, intent(inout) :: time_loop_finished !< flag to indicate end of timeloop
     337     2020272 :     time_loop_finished = run_cfg%time_step == run_cfg%domainDateTime%nTimeSteps
     338          82 :   end subroutine mhm_interface_run_finished
     339             : 
     340             :   !> \brief do one time-step on current domain
     341     2020272 :   subroutine mhm_interface_run_do_time_step()
     342             :     implicit none
     343             : 
     344             :     integer(i4) :: iDomain, domainID, tt, jj, s1, e1
     345             : 
     346             :     ! increment time step count (first input is 0)
     347     2020272 :     run_cfg%time_step = run_cfg%time_step + 1_i4
     348             :     ! current time counter and domain indices
     349     2020272 :     tt = run_cfg%time_step
     350     2020272 :     s1 = run_cfg%s1
     351     2020272 :     e1 = run_cfg%e1
     352     2020272 :     iDomain = run_cfg%get_domain_index(run_cfg%selected_domain)
     353     2020272 :     domainID = domainMeta%indices(iDomain)
     354             : 
     355     2020272 :     call run_cfg%domainDateTime%update_LAI_timestep()
     356             : 
     357             :     ! update the meteo-handler
     358             :     call meteo_handler%update_timestep( &
     359             :       tt=tt, &
     360             :       time=run_cfg%domainDateTime%newTime - 0.5_dp, &
     361             :       iDomain=iDomain, &
     362             :       level1=level1, &
     363     2020272 :       simPer=simPer)
     364             : 
     365             :     ! get corrected pet
     366           0 :     call meteo_handler%get_corrected_pet(pet_calc=L1_pet_calc(s1 : e1), &
     367             :       ! pet calculation dependencies
     368           0 :       petLAIcorFactorL1 = L1_petLAIcorFactor(s1 : e1, run_cfg%domainDateTime%iLAI, run_cfg%domainDateTime%yId), &
     369           0 :       fAsp              = L1_fAsp(s1 : e1, 1, 1), &
     370           0 :       HarSamCoeff       = L1_HarSamCoeff(s1 : e1, 1, 1), &
     371           0 :       latitude          = pack(level1(iDomain)%y, level1(iDomain)%mask), &
     372           0 :       PrieTayAlpha      = L1_PrieTayAlpha(s1 : e1, run_cfg%domainDateTime%iLAI, 1), &
     373           0 :       aeroResist        = L1_aeroResist(s1 : e1, run_cfg%domainDateTime%iLAI, run_cfg%domainDateTime%yId), &
     374     2020272 :       surfResist        = L1_surfResist(s1 : e1, run_cfg%domainDateTime%iLAI, 1))
     375             : 
     376             :     ! get temperature and precipitation
     377     2020272 :     call meteo_handler%get_temp(temp_calc=L1_temp_calc(s1 : e1))
     378     2020272 :     call meteo_handler%get_prec(prec_calc=L1_prec_calc(s1 : e1))
     379             : 
     380             :     ! -------------------------------------------------------------------------
     381             :     ! ARGUMENT LIST KEY FOR mHM
     382             :     ! -------------------------------------------------------------------------
     383             :     !  C    CONFIGURATION
     384             :     !  F    FORCING DATA L2
     385             :     !  Q    INFLOW FROM UPSTREAM AREAS
     386             :     !  L0   MORPHOLOGIC DATA L0
     387             :     !  L1   MORPHOLOGIC DATA L1
     388             :     !  L11  MORPHOLOGIC DATA L11
     389             :     !  P    GLOBAL PARAMETERS
     390             :     !  E    EFFECTIVE PARAMETER FIELDS (L1, L11 levels)
     391             :     !  S    STATE VARIABLES L1
     392             :     !  X    FLUXES (L1, L11 levels)
     393             :     ! --------------------------------------------------------------------------
     394             :     call mhm( &
     395             :       read_states            = read_restart, & ! IN C
     396             :       tt                     = tt, &
     397             :       time                   = run_cfg%domainDateTime%newTime - 0.5_dp, &
     398             :       processMatrix          = processMatrix, &
     399             :       horizon_depth          = HorizonDepth_mHM, & ! IN C
     400             :       nCells1                = run_cfg%nCells, &
     401             :       nHorizons_mHM          = nSoilHorizons_mHM, &
     402             :       c2TSTu                 = c2TSTu,  & ! IN C
     403             :       neutron_integral_AFast = neutron_integral_AFast, & ! IN C
     404             :       evap_coeff             = evap_coeff, &
     405           0 :       fSealed1               = L1_fSealed(s1 : e1, 1, run_cfg%domainDateTime%yId), & ! INOUT L1
     406           0 :       interc                 = L1_inter(s1 : e1), &
     407           0 :       snowpack               = L1_snowPack(s1 : e1), &
     408           0 :       sealedStorage          = L1_sealSTW(s1 : e1), & ! INOUT S
     409           0 :       soilMoisture           = L1_soilMoist(s1 : e1, :), &
     410           0 :       unsatStorage           = L1_unsatSTW(s1 : e1), &
     411           0 :       satStorage             = L1_satSTW(s1 : e1), & ! INOUT S
     412           0 :       neutrons               = L1_neutrons(s1 : e1), & ! INOUT S
     413           0 :       pet_calc               = L1_pet_calc(s1 : e1), & ! INOUT X
     414           0 :       temp_calc              = L1_temp_calc(s1 : e1), & ! INOUT X
     415           0 :       prec_calc              = L1_prec_calc(s1 : e1), & ! INOUT X
     416           0 :       aet_soil               = L1_aETSoil(s1 : e1, :), &
     417           0 :       aet_canopy             = L1_aETCanopy(s1 : e1), &
     418           0 :       aet_sealed             = L1_aETSealed(s1 : e1), & ! INOUT X
     419           0 :       baseflow               = L1_baseflow(s1 : e1), &
     420           0 :       infiltration           = L1_infilSoil(s1 : e1, :), &
     421           0 :       fast_interflow         = L1_fastRunoff(s1 : e1), & ! INOUT X
     422           0 :       melt                   = L1_melt(s1 : e1), &
     423           0 :       perc                   = L1_percol(s1 : e1), &
     424           0 :       prec_effect            = L1_preEffect(s1 : e1), &
     425           0 :       rain                   = L1_rain(s1 : e1), & ! INOUT X
     426           0 :       runoff_sealed          = L1_runoffSeal(s1 : e1), &
     427           0 :       slow_interflow         = L1_slowRunoff(s1 : e1), &
     428           0 :       snow                   = L1_snow(s1 : e1), & ! INOUT X
     429           0 :       throughfall            = L1_Throughfall(s1 : e1), &
     430           0 :       total_runoff           = L1_total_runoff(s1 : e1), & ! INOUT X
     431             :       ! MPR
     432           0 :       alpha                  = L1_alpha(s1 : e1, 1, run_cfg%domainDateTime%yId), &
     433           0 :       deg_day_incr           = L1_degDayInc(s1 : e1, 1, run_cfg%domainDateTime%yId), &
     434           0 :       deg_day_max            = L1_degDayMax(s1 : e1, 1, run_cfg%domainDateTime%yId), & ! INOUT E1
     435           0 :       deg_day_noprec         = L1_degDayNoPre(s1 : e1, 1, run_cfg%domainDateTime%yId), &
     436           0 :       deg_day                = L1_degDay(s1 : e1, 1, 1), & ! INOUT E1
     437           0 :       frac_roots             = L1_fRoots(s1 : e1, :, run_cfg%domainDateTime%yId), & ! INOUT E1
     438           0 :       interc_max             = L1_maxInter(s1 : e1, run_cfg%domainDateTime%iLAI, 1), &
     439           0 :       karst_loss             = L1_karstLoss(s1 : e1, 1, 1), & ! INOUT E1
     440           0 :       k0                     = L1_kFastFlow(s1 : e1, 1, run_cfg%domainDateTime%yId), &
     441           0 :       k1                     = L1_kSlowFlow(s1 : e1, 1, run_cfg%domainDateTime%yId), & ! INOUT E1
     442           0 :       k2                     = L1_kBaseFlow(s1 : e1, 1, run_cfg%domainDateTime%yId), &
     443           0 :       kp                     = L1_kPerco(s1 : e1, 1, run_cfg%domainDateTime%yId), & ! INOUT E1
     444           0 :       soil_moist_FC          = L1_soilMoistFC(s1 : e1, :, run_cfg%domainDateTime%yId), & ! INOUT E1
     445           0 :       soil_moist_sat         = L1_soilMoistSat(s1 : e1, :, run_cfg%domainDateTime%yId), & ! INOUT E1
     446           0 :       soil_moist_exponen     = L1_soilMoistExp(s1 : e1, :, run_cfg%domainDateTime%yId), &
     447           0 :       jarvis_thresh_c1       = L1_jarvis_thresh_c1(s1 : e1, 1, 1), & ! INOUT E1
     448           0 :       temp_thresh            = L1_tempThresh(s1 : e1, 1, run_cfg%domainDateTime%yId), &
     449           0 :       unsat_thresh           = L1_unsatThresh(s1 : e1, 1, 1), & ! INOUT E1
     450           0 :       water_thresh_sealed    = L1_sealedThresh(s1 : e1, 1, 1), & ! INOUT E1
     451           0 :       wilting_point          = L1_wiltingPoint(s1 : e1, :, run_cfg%domainDateTime%yId), & ! INOUT E1
     452             :       ! >> neutron count
     453           0 :       No_count               = L1_No_Count(s1:e1, 1, 1),  &                     ! INOUT E1
     454           0 :       bulkDens               = L1_bulkDens(s1:e1,     :, run_cfg%domainDateTime%yId), & ! INOUT E1
     455           0 :       latticeWater           = L1_latticeWater(s1:e1, :, run_cfg%domainDateTime%yId), & ! INOUT E1
     456           0 :       COSMICL3               = L1_COSMICL3(s1:e1,     :, run_cfg%domainDateTime%yId)  & ! INOUT E1
     457     2020272 :     )
     458             : 
     459             :     ! call mRM routing
     460     2020272 :     run_cfg%doRoute = .false.
     461     2020272 :     if (domainMeta%doRouting(iDomain)) then
     462             :       ! set discharge timestep
     463      779928 :       run_cfg%iDischargeTS = ceiling(real(tt, dp) / real(nTstepDay, dp))
     464             :       ! set input variables for routing
     465      779928 :       if (processMatrix(8, 1) .eq. 1) then
     466             :         ! >>>
     467             :         ! >>> original Muskingum routing, executed every time
     468             :         ! >>>
     469      495096 :         run_cfg%doRoute = .True.
     470      495096 :         run_cfg%tsRoutFactorIn = 1._dp
     471      495096 :         run_cfg%timestep_rout = timestep
     472    25781256 :         run_cfg%RunToRout = L1_total_runoff(s1 : e1) ! runoff [mm TS-1] mm per timestep
     473     1524600 :         run_cfg%InflowDischarge = InflowGauge%Q(run_cfg%iDischargeTS, :) ! inflow discharge in [m3 s-1]
     474             :         !
     475      284832 :       else if ((processMatrix(8, 1) .eq. 2) .or. &
     476             :                (processMatrix(8, 1) .eq. 3)) then
     477             :         ! >>>
     478             :         ! >>> adaptive timestep
     479             :         ! >>>
     480             :         run_cfg%doRoute = .False.
     481             :         ! calculate factor
     482      284832 :         run_cfg%tsRoutFactor = L11_tsRout(iDomain) / (timestep * HourSecs)
     483             :         ! print *, 'routing factor: ', tsRoutFactor
     484             :         ! prepare routing call
     485      284832 :         if (run_cfg%tsRoutFactor .lt. 1._dp) then
     486             :           ! ----------------------------------------------------------------
     487             :           ! routing timesteps are shorter than hydrologic time steps
     488             :           ! ----------------------------------------------------------------
     489             :           ! set all input variables
     490           0 :           run_cfg%tsRoutFactorIn = run_cfg%tsRoutFactor
     491           0 :           run_cfg%RunToRout = L1_total_runoff(s1 : e1) ! runoff [mm TS-1] mm per timestep
     492           0 :           run_cfg%InflowDischarge = InflowGauge%Q(run_cfg%iDischargeTS, :) ! inflow discharge in [m3 s-1]
     493           0 :           run_cfg%timestep_rout = timestep
     494           0 :           run_cfg%doRoute = .True.
     495             :         else
     496             :           ! ----------------------------------------------------------------
     497             :           ! routing timesteps are longer than hydrologic time steps
     498             :           ! ----------------------------------------------------------------
     499             :           ! set all input variables
     500      284832 :           run_cfg%tsRoutFactorIn = run_cfg%tsRoutFactor
     501             :           ! Runoff is accumulated in [mm]
     502    10253952 :           run_cfg%RunToRout = run_cfg%RunToRout + L1_total_runoff(s1 : e1)
     503      854496 :           run_cfg%InflowDischarge = run_cfg%InflowDischarge + InflowGauge%Q(run_cfg%iDischargeTS, :)
     504             :           ! reset tsRoutFactorIn if last period did not cover full period
     505      284832 :           if ((tt == run_cfg%domainDateTime%nTimeSteps) .and. (mod(tt, nint(run_cfg%tsRoutFactorIn)) /= 0_i4)) &
     506           0 :             run_cfg%tsRoutFactorIn = mod(tt, nint(run_cfg%tsRoutFactorIn))
     507      284832 :           if ((mod(tt, nint(run_cfg%tsRoutFactorIn)) .eq. 0_i4) .or. (tt .eq. run_cfg%domainDateTime%nTimeSteps)) then
     508             :             ! Inflow discharge is given as flow-rate and has to be converted to [m3]
     509      278280 :             run_cfg%InflowDischarge = run_cfg%InflowDischarge / run_cfg%tsRoutFactorIn
     510      139140 :             run_cfg%timestep_rout = timestep * nint(run_cfg%tsRoutFactorIn, i4)
     511      139140 :             run_cfg%doRoute = .True.
     512             :           end if
     513             :         end if
     514             :       end if
     515             :       ! prepare temperature routing
     516      779928 :       if ( riv_temp_pcs%active ) then
     517             :         ! init riv-temp from current air temp
     518       13104 :         if ( tt .eq. 1_i4 ) call riv_temp_pcs%init_riv_temp( &
     519           0 :           temp_air     = L1_temp_calc(s1 : e1), &
     520           0 :           efecarea     = level1(iDomain)%CellArea * 1.E-6_dp, &
     521           0 :           L1_L11_Id    = L1_L11_Id(s1 : e1), &
     522           0 :           L11_areacell = level11(iDomain)%CellArea * 1.E-6_dp, &
     523           0 :           L11_L1_Id    = L11_L1_Id(run_cfg%s11 : run_cfg%e11), &
     524           0 :           map_flag     = ge(resolutionRouting(iDomain), resolutionHydrology(iDomain)) &
     525          69 :         )
     526             :         ! get riv-temp specific meteo arrays
     527       13104 :         call meteo_handler%get_ssrd(riv_temp_pcs%L1_ssrd_calc)
     528       13104 :         call meteo_handler%get_strd(riv_temp_pcs%L1_strd_calc)
     529       13104 :         call meteo_handler%get_tann(riv_temp_pcs%L1_tann_calc)
     530             :         ! accumulate source Energy at L1 level
     531             :         call riv_temp_pcs%acc_source_E( &
     532           0 :           fSealed_area_fraction = L1_fSealed(s1 : e1, 1, run_cfg%domainDateTime%yId), &
     533           0 :           fast_interflow        = L1_fastRunoff(s1 : e1), &
     534           0 :           slow_interflow        = L1_slowRunoff(s1 : e1), &
     535           0 :           baseflow              = L1_baseflow(s1 : e1), &
     536           0 :           direct_runoff         = L1_runoffSeal(s1 : e1), &
     537           0 :           temp_air              = L1_temp_calc(s1 : e1) &
     538       13104 :         )
     539             :         ! if routing should be performed, scale source energy to L11 level
     540       13104 :         if ( run_cfg%doRoute ) call riv_temp_pcs%finalize_source_E( &
     541           0 :           efecarea     = level1(iDomain)%CellArea * 1.E-6_dp, &
     542           0 :           L1_L11_Id    = L1_L11_Id(s1 : e1), &
     543           0 :           L11_areacell = level11(iDomain)%CellArea * 1.E-6_dp, &
     544           0 :           L11_L1_Id    = L11_L1_Id(run_cfg%s11 : run_cfg%e11), &
     545             :           timestep     = run_cfg%timestep_rout, &
     546           0 :           map_flag     = ge(resolutionRouting(iDomain), resolutionHydrology(iDomain)) &
     547      452088 :         )
     548             :       end if
     549             :       ! -------------------------------------------------------------------
     550             :       ! execute routing
     551             :       ! -------------------------------------------------------------------
     552      779928 :       if (run_cfg%doRoute) call mRM_routing(&
     553             :         ! general INPUT variables
     554             :         read_restart, &
     555             :         processMatrix(8, 1), & ! parse process Case to be used
     556           0 :         run_cfg%parameterset(processMatrix(8, 3) - processMatrix(8, 2) + 1 : processMatrix(8, 3)), & ! routing par.
     557             :         run_cfg%RunToRout, & ! runoff [mm TS-1] mm per timestep old: L1_total_runoff_in(run_cfg%s1:run_cfg%e1, tt), &
     558           0 :         level1(iDomain)%CellArea * 1.E-6_dp, &
     559           0 :         L1_L11_Id(s1 : e1), &
     560           0 :         level11(iDomain)%CellArea * 1.E-6_dp, &
     561           0 :         L11_L1_Id(run_cfg%s11 : run_cfg%e11), &
     562           0 :         L11_netPerm(run_cfg%s11 : run_cfg%e11), & ! routing order at L11
     563           0 :         L11_fromN(run_cfg%s11 : run_cfg%e11), & ! link source at L11
     564           0 :         L11_toN(run_cfg%s11 : run_cfg%e11), & ! link target at L11
     565           0 :         L11_nOutlets(iDomain), & ! number of outlets
     566             :         run_cfg%timestep_rout, & ! timestep of runoff to rout [h]
     567             :         run_cfg%tsRoutFactorIn, & ! simulate timestep in [h]
     568           0 :         level11(iDomain)%nCells, & ! number of Nodes
     569           0 :         domain_mrm(iDomain)%nInflowGauges, &
     570             :         domain_mrm(iDomain)%InflowGaugeIndexList(:), &
     571             :         domain_mrm(iDomain)%InflowGaugeHeadwater(:), &
     572             :         domain_mrm(iDomain)%InflowGaugeNodeList(:), &
     573             :         run_cfg%InflowDischarge, &
     574             :         domain_mrm(iDomain)%nGauges, &
     575             :         domain_mrm(iDomain)%gaugeIndexList(:), &
     576             :         domain_mrm(iDomain)%gaugeNodeList(:), &
     577           0 :         ge(resolutionRouting(iDomain), resolutionHydrology(iDomain)), &
     578             :         ! original routing specific input variables
     579           0 :         L11_length(run_cfg%s11 : run_cfg%e11 - 1), & ! link length
     580           0 :         L11_slope(run_cfg%s11 : run_cfg%e11 - 1), &
     581           0 :         L11_nLinkFracFPimp(run_cfg%s11 : run_cfg%e11, run_cfg%domainDateTime%yId), & ! fraction of impervious layer at L11 scale
     582             :         ! general INPUT/OUTPUT variables
     583           0 :         L11_C1(run_cfg%s11 : run_cfg%e11), & ! first muskingum parameter
     584           0 :         L11_C2(run_cfg%s11 : run_cfg%e11), & ! second muskigum parameter
     585           0 :         L11_qOUT(run_cfg%s11 : run_cfg%e11), & ! routed runoff flowing out of L11 cell
     586           0 :         L11_qTIN(run_cfg%s11 : run_cfg%e11, :), & ! inflow water into the reach at L11
     587           0 :         L11_qTR(run_cfg%s11 : run_cfg%e11, :), & !
     588           0 :         L11_qMod(run_cfg%s11 : run_cfg%e11), &
     589             :         mRM_runoff(tt, :) &
     590    58038084 :       )
     591             :       ! -------------------------------------------------------------------
     592             :       ! groundwater coupling
     593             :       ! -------------------------------------------------------------------
     594      779928 :       if (gw_coupling) then
     595             :         call calc_river_head(iDomain, L11_Qmod, L0_river_head_mon_sum)
     596             :       end if
     597             :       ! -------------------------------------------------------------------
     598             :       ! reset variables
     599             :       ! -------------------------------------------------------------------
     600      779928 :       if (processMatrix(8, 1) .eq. 1) then
     601             :         ! reset Input variables
     602     1029504 :         run_cfg%InflowDischarge = 0._dp
     603    25286160 :         run_cfg%RunToRout = 0._dp
     604      284832 :       else if ((processMatrix(8, 1) .eq. 2) .or. (processMatrix(8, 1) .eq. 3)) then
     605      284832 :         if ((.not. (run_cfg%tsRoutFactorIn .lt. 1._dp)) .and. run_cfg%doRoute) then
     606      423972 :           do jj = 1, nint(run_cfg%tsRoutFactorIn) ! BUG: this should start at 2
     607      708804 :             mRM_runoff(tt - jj + 1, :) = mRM_runoff(tt, :)
     608             :           end do
     609             :           ! reset Input variables
     610      278280 :           run_cfg%InflowDischarge = 0._dp
     611     4869900 :           run_cfg%RunToRout = 0._dp
     612             :           ! reset lateral fluxes and time-step counter if routing was done
     613      139140 :           if ( riv_temp_pcs%active ) call riv_temp_pcs%reset_timestep()
     614             :         end if
     615             :         ! if routing is done every time-step, reset river-temp time step
     616      284832 :         if (run_cfg%tsRoutFactor .lt. 1._dp .and. riv_temp_pcs%active ) call riv_temp_pcs%reset_timestep()
     617             :       end if
     618             :     end if
     619             : 
     620             :     ! output only for evaluation period
     621     2020272 :     run_cfg%domainDateTime%tIndex_out = (tt - warmingDays(iDomain) * nTstepDay) ! tt if write out of warming period
     622             : 
     623     2020272 :     call run_cfg%domainDateTime%increment()
     624             : 
     625             :     ! update the year-dependent domainDateTime%yId (land cover id)
     626     2020272 :     if (run_cfg%domainDateTime%is_new_year .and. tt < run_cfg%domainDateTime%nTimeSteps) then
     627         165 :       run_cfg%domainDateTime%yId = LCyearId(run_cfg%domainDateTime%year, iDomain)
     628             :     end if
     629             : 
     630             :     ! calculate BFI releated after warming days if wanted
     631     2020272 :     if ( run_cfg%output_BFI .and. (run_cfg%domainDateTime%tIndex_out > 0_i4) ) then
     632           0 :       BFI_qBF_sum(iDomain) = BFI_qBF_sum(iDomain) &
     633           0 :         + sum(L1_baseflow(s1 : e1) * level1(iDomain)%CellArea) / level1(iDomain)%nCells
     634           0 :       BFI_qT_sum(iDomain) = BFI_qT_sum(iDomain) &
     635           0 :         + sum(L1_total_runoff(s1 : e1) * level1(iDomain)%CellArea) / level1(iDomain)%nCells
     636             :     end if
     637             : 
     638     4040544 :   end subroutine mhm_interface_run_do_time_step
     639             : 
     640             :   !> \brief write output after current time-step
     641     2020272 :   subroutine mhm_interface_run_write_output()
     642             :     implicit none
     643             :     integer(i4) :: iDomain, tt, s0, e0
     644             : 
     645             :     ! get time step
     646     2020272 :     tt = run_cfg%time_step
     647             : 
     648             :     ! get domain index
     649     2020272 :     iDomain = run_cfg%get_domain_index(run_cfg%selected_domain)
     650             : 
     651     2020272 :     if ( (.not. optimize) .and. (run_cfg%domainDateTime%tIndex_out > 0_i4)) then
     652      350760 :       if (any(outputFlxState_mrm) .AND. (domainMeta%doRouting(iDomain))) then
     653             : 
     654       35040 :         if (run_cfg%domainDateTime%tIndex_out == 1) then
     655           4 :           run_cfg%nc_mrm = mRM_OutputDataset(iDomain, run_cfg%mask11)
     656             :         end if
     657             : 
     658             :         ! update Dataset (riv-temp as optional input)
     659       35040 :         if ( riv_temp_pcs%active ) then
     660             :           call mRM_updateDataset(run_cfg%nc_mrm, &
     661        8760 :             L11_qmod(run_cfg%s11 : run_cfg%e11), riv_temp_pcs%river_temp(riv_temp_pcs%s11 : riv_temp_pcs%e11))
     662             :         else
     663       26280 :           call mRM_updateDataset(run_cfg%nc_mrm, L11_qmod(run_cfg%s11 : run_cfg%e11))
     664             :         end if
     665             : 
     666             :         ! write data
     667       35040 :         if (run_cfg%domainDateTime%writeout(timeStep_model_outputs_mrm, tt)) then
     668         732 :           call run_cfg%nc_mrm%writeTimestep(run_cfg%domainDateTime%tIndex_out * timestep)
     669             :         end if
     670             : 
     671       35040 :         if(tt == run_cfg%domainDateTime%nTimeSteps) then
     672           4 :           call run_cfg%nc_mrm%close()
     673             :         end if
     674             : 
     675       35040 :         if ( gw_coupling ) then
     676             :           ! create
     677           0 :           if (run_cfg%domainDateTime%tIndex_out == 1) then
     678           0 :             run_cfg%nc_gw = GW_OutputDataset(iDomain, level0(iDomain)%mask)
     679             :           end if
     680             :           ! add data
     681           0 :           s0 = level0(iDomain)%iStart
     682           0 :           e0 = level0(iDomain)%iEnd
     683           0 :           call GW_updateDataset(run_cfg%nc_gw, L0_river_head_mon_sum(s0 : e0))
     684             :           ! write
     685           0 :           if (run_cfg%domainDateTime%writeout(-2, tt)) then ! -2 for monthly
     686           0 :             call run_cfg%nc_gw%writeTimestep(run_cfg%domainDateTime%tIndex_out * timestep)
     687             :           end if
     688             :           ! close
     689           0 :           if(tt == run_cfg%domainDateTime%nTimeSteps) then
     690           0 :             call run_cfg%nc_gw%close()
     691             :           end if
     692             :         end if
     693             : 
     694             :       end if
     695             : 
     696      587280 :       if (any(outputFlxState)) then
     697             : 
     698      131520 :         if (run_cfg%domainDateTime%tIndex_out == 1) then
     699          38 :           run_cfg%nc_mhm = mHM_OutputDataset(iDomain, run_cfg%mask1)
     700             :         end if
     701             : 
     702             :         call mHM_updateDataset(&
     703             :           run_cfg%nc_mhm, &
     704             :           run_cfg%s1, run_cfg%e1, &
     705             :           L1_fSealed(:, 1, run_cfg%domainDateTime%yId), &
     706             :           run_cfg%L1_fNotSealed(:, 1, run_cfg%domainDateTime%yId), &
     707             :           L1_inter, &
     708             :           L1_snowPack, &
     709             :           L1_soilMoist, &
     710             :           L1_soilMoistSat(:, :, run_cfg%domainDateTime%yId), &
     711             :           L1_sealSTW, &
     712             :           L1_unsatSTW, &
     713             :           L1_satSTW, &
     714             :           L1_neutrons, &
     715             :           L1_pet_calc, &
     716             :           L1_aETSoil, &
     717             :           L1_aETCanopy, &
     718             :           L1_aETSealed, &
     719             :           L1_total_runoff, &
     720             :           L1_runoffSeal, &
     721             :           L1_fastRunoff, &
     722             :           L1_slowRunoff, &
     723             :           L1_baseflow, &
     724             :           L1_percol, &
     725             :           L1_infilSoil, &
     726             :           L1_preEffect, &
     727             :           L1_melt &
     728      131520 :         )
     729             : 
     730             :         ! write data
     731      131520 :         if (run_cfg%domainDateTime%writeout(timeStep_model_outputs, tt)) then
     732          93 :           call run_cfg%nc_mhm%writeTimestep(run_cfg%domainDateTime%tIndex_out * timestep)
     733             :         end if
     734             : 
     735      131520 :         if(tt == run_cfg%domainDateTime%nTimeSteps) then
     736          15 :           call run_cfg%nc_mhm%close()
     737             :         end if
     738             : 
     739             :       end if
     740             :     end if ! <-- if (.not. optimize)
     741             : 
     742     2020272 :   end subroutine mhm_interface_run_write_output
     743             : 
     744             :   !> \brief add simulation data to optimization data types
     745     4040544 :   subroutine mhm_interface_run_update_optisim(etOptiSim, twsOptiSim, neutronsOptiSim, smOptiSim)
     746             :     implicit none
     747             :     !> returns soil moisture time series for all grid cells (of multiple Domains concatenated),DIMENSION [nCells, nTimeSteps]
     748             :     type(optidata_sim), dimension(:), optional, intent(inout) :: smOptiSim
     749             :     !> dim1=ncells, dim2=time
     750             :     type(optidata_sim), dimension(:), optional, intent(inout) :: neutronsOptiSim
     751             :     !> returns evapotranspiration time series for all grid cells (of multiple Domains concatenated),DIMENSION [nCells, nTimeSteps]
     752             :     type(optidata_sim), dimension(:), optional, intent(inout) :: etOptiSim
     753             :     !> returns tws time series for all grid cells (of multiple Domains concatenated),DIMENSION [nCells, nTimeSteps]
     754             :     type(optidata_sim), dimension(:), optional, intent(inout) :: twsOptiSim
     755             : 
     756             :     integer(i4) :: gg, s1, e1, lcId
     757             : 
     758             :     integer(i4) :: iDomain, tt
     759             : 
     760             :     ! get time step
     761     2020272 :     tt = run_cfg%time_step
     762             : 
     763             :     ! get domain index
     764     2020272 :     iDomain = run_cfg%get_domain_index(run_cfg%selected_domain)
     765             : 
     766             :     ! slice settings
     767     2020272 :     s1 = run_cfg%s1
     768     2020272 :     e1 = run_cfg%e1
     769     2020272 :     lcId = run_cfg%domainDateTime%yId
     770             : 
     771             :     !----------------------------------------------------------------------
     772             :     ! FOR SOIL MOISTURE
     773             :     ! NOTE:: modeled soil moisture is averaged according to input time step
     774             :     !        soil moisture (timeStep_sm_input)
     775             :     !----------------------------------------------------------------------
     776     2020272 :     if (present(smOptiSim)) then
     777             :       ! only for evaluation period - ignore warming days
     778       78624 :       if ((tt - warmingDays(iDomain) * nTstepDay) .GT. 0) then
     779             :         ! decide for daily, monthly or yearly aggregation
     780       52560 :         call smOptiSim(iDomain)%average_per_timestep(L1_smObs(iDomain)%timeStepInput, &
     781      105120 :                      run_cfg%domainDateTime%is_new_day, run_cfg%domainDateTime%is_new_month, run_cfg%domainDateTime%is_new_year)
     782             :         ! last timestep is already done - write_counter exceeds size(smOptiSim(iDomain)%dataSim, dim=2)
     783       52560 :         if (tt /= run_cfg%domainDateTime%nTimeSteps) then
     784             :           ! aggregate soil moisture to needed time step for optimization
     785       52554 :           call smOptiSim(iDomain)%average_add(sum(L1_soilMoist(:, 1 : nSoilHorizons_sm_input), dim = 2) / &
     786     5465616 :                           sum(L1_soilMoistSat(:, 1 : nSoilHorizons_sm_input, lcId), dim = 2))
     787             :         end if
     788             :       end if
     789             :     end if
     790             : 
     791             :     !----------------------------------------------------------------------
     792             :     ! FOR NEUTRONS
     793             :     ! NOTE:: modeled neutrons are averaged daily
     794             :     !----------------------------------------------------------------------
     795     2020272 :     if (present(neutronsOptiSim)) then
     796             :       ! only for evaluation period - ignore warming days
     797       78624 :       if ((tt - warmingDays(iDomain) * nTstepDay) .GT. 0) then
     798             :         ! decide for daily, monthly or yearly aggregation
     799             :         ! daily
     800       52560 :         if (run_cfg%domainDateTime%is_new_day)   then
     801        2190 :           call neutronsOptiSim(iDomain)%average()
     802             :         end if
     803             : 
     804             :         ! last timestep is already done - write_counter exceeds size(sm_opti, dim=2)
     805       52560 :         if (tt /= run_cfg%domainDateTime%nTimeSteps) then
     806             :           ! aggregate neutrons to needed time step for optimization
     807       52554 :           call neutronsOptiSim(iDomain)%average_add(L1_neutrons(s1 : e1))
     808             :         end if
     809             :       end if
     810             :     end if
     811             : 
     812             :     !----------------------------------------------------------------------
     813             :     ! FOR EVAPOTRANSPIRATION
     814             :     ! NOTE:: modeled evapotranspiration is averaged according to input time step
     815             :     !        evapotranspiration (timeStep_et_input)
     816             :     !----------------------------------------------------------------------
     817     2020272 :     if (present(etOptiSim)) then
     818             :       ! only for evaluation period - ignore warming days
     819      736344 :       if ((tt - warmingDays(iDomain) * nTstepDay) .GT. 0) then
     820             :         ! decide for daily, monthly or yearly aggregation
     821      736344 :         call etOptiSim(iDomain)%increment_counter(L1_etObs(iDomain)%timeStepInput, &
     822     1472688 :           run_cfg%domainDateTime%is_new_day, run_cfg%domainDateTime%is_new_month, run_cfg%domainDateTime%is_new_year)
     823             : 
     824             :         ! last timestep is already done - write_counter exceeds size(etOptiSim(iDomain)%dataSim, dim=2)
     825      736344 :         if (tt /= run_cfg%domainDateTime%nTimeSteps) then
     826             :           ! aggregate evapotranspiration to needed time step for optimization
     827           0 :           call etOptiSim(iDomain)%add(sum(L1_aETSoil(s1 : e1, :), dim = 2) * &
     828           0 :             run_cfg%L1_fNotSealed(s1 : e1, 1, lcId) + &
     829           0 :             L1_aETCanopy(s1 : e1) + &
     830    75841269 :             L1_aETSealed(s1 : e1) * L1_fSealed(s1 : e1, 1, lcId))
     831             :         end if
     832             :       end if
     833             :     end if
     834             : 
     835             :     !----------------------------------------------------------------------
     836             :     ! FOR TWS
     837             :     ! NOTE:: modeled tws is averaged according to input time step
     838             :     !        (timeStepInput)
     839             :     !----------------------------------------------------------------------
     840     2020272 :     if (present(twsOptiSim)) then
     841             :       ! only for evaluation period - ignore warming days
     842      841464 :       if ((tt - warmingDays(iDomain) * nTstepDay) > 0) then
     843             :         ! decide for daily, monthly or yearly aggregation
     844      841464 :         call twsOptiSim(iDomain)%average_per_timestep(L1_twsaObs(iDomain)%timeStepInput, &
     845     1682928 :           run_cfg%domainDateTime%is_new_day, run_cfg%domainDateTime%is_new_month, run_cfg%domainDateTime%is_new_year)
     846             : 
     847             :         ! last timestep is already done - write_counter exceeds size(twsOptiSim(iDomain)%dataSim, dim=2)
     848      841464 :         if (tt /= run_cfg%domainDateTime%nTimeSteps) then
     849             :           ! aggregate evapotranspiration to needed time step for optimization
     850           0 :           call twsOptiSim(iDomain)%average_add(L1_inter(s1 : e1) + L1_snowPack(s1 : e1) + L1_sealSTW(s1 : e1) + &
     851    29450295 :             L1_unsatSTW(s1 : e1) + L1_satSTW(s1 : e1))
     852     2524311 :           do gg = 1, nSoilHorizons_mHM
     853     2524311 :             call twsOptiSim(iDomain)%add(L1_soilMoist (s1 : e1, gg))
     854             :           end do
     855             :         end if
     856             :       end if
     857             :     end if
     858             : 
     859             :     ! TODO-RIV-TEMP: add OptiSim for river temperature
     860             : 
     861     2020272 :   end subroutine mhm_interface_run_update_optisim
     862             : 
     863             :   !> \brief finalize current domain after running
     864          82 :   subroutine mhm_interface_run_finalize_domain()
     865             :     implicit none
     866             : 
     867             :     integer(i4) :: iDomain
     868             : 
     869             :     ! get domain index
     870          82 :     iDomain = run_cfg%get_domain_index(run_cfg%selected_domain)
     871             : 
     872          46 :     if (allocated(run_cfg%InflowDischarge)) deallocate(run_cfg%InflowDischarge)
     873          82 :       if (domainMeta%doRouting(iDomain)) then
     874             :       ! clean runoff variable
     875          46 :       deallocate(run_cfg%RunToRout)
     876          46 :       if ( riv_temp_pcs%active ) call riv_temp_pcs%dealloc_lateral()
     877             :     end if
     878             : 
     879     2020272 :   end subroutine mhm_interface_run_finalize_domain
     880             : 
     881             :   !> \brief finalize run
     882          61 :   subroutine mhm_interface_run_finalize(runoff, BFI)
     883             :     implicit none
     884             :     !> returns runoff time series, DIMENSION [nTimeSteps, nGaugesTotal]
     885             :     real(dp), dimension(:, :), allocatable, optional, intent(out) :: runoff
     886             :     !> baseflow index, dim1=domainID
     887             :     real(dp), dimension(:), allocatable, optional, intent(out) :: BFI
     888             : 
     889      508027 :     if (present(runoff) .and. (processMatrix(8, 1) > 0)) runoff = mRM_runoff
     890          61 :     if (present(BFI)) then
     891           0 :       BFI = BFI_qBF_sum / BFI_qT_sum
     892           0 :       deallocate(BFI_qBF_sum)
     893           0 :       deallocate(BFI_qT_sum)
     894             :     end if
     895             : 
     896          61 :     if (allocated(run_cfg%parameterset)) deallocate(run_cfg%parameterset)
     897          61 :     if (allocated(run_cfg%L1_fNotSealed)) deallocate(run_cfg%L1_fNotSealed)
     898          61 :     if (allocated(run_cfg%domain_indices)) deallocate(run_cfg%domain_indices)
     899             : 
     900          82 :   end subroutine mhm_interface_run_finalize
     901             : 
     902             : end module mo_mhm_interface_run

Generated by: LCOV version 1.16