mo_write_fluxes_states.f90

Go to the documentation of this file.

!> \file mo_write_fluxes_states.f90
!> \brief \copybrief mo_write_fluxes_states
!> \details \copydetails mo_write_fluxes_states
 
!> \brief Creates NetCDF output for different fluxes and state variables of mHM.
!> \details NetCDF is first initialized and later on variables are put to the NetCDF.
!> \changelog
!! - David Schaefer       Aug 2015
!!   - major rewrite
!! - O. Rakovec, R. Kumar Nov 2017
!!   - added project description for the netcdf outputs
!> \authors Matthias Zink
!> \date Apr 2013
!> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
!! mHM is released under the LGPLv3+ license \license_note
!> \ingroup f_mhm
module mo_write_fluxes_states
 
  use mo_nc_output, only: outputdataset, outputvariable, set_attributes, data_dims, data_dtype
  use mo_kind, only : i4, dp
  use mo_string_utils, only : num2str
  use mo_common_constants, only : nodata_dp
  use mo_netcdf, only : ncdataset
  use mo_global_variables, only : &
    output_deflate_level, output_double_precision, output_time_reference, timestep_model_outputs, outputflxstate
  use mo_common_mhm_mrm_variables, only: timestep
  use mo_mpr_global_variables, only : nsoilhorizons_mhm
  use mo_common_variables, only : iflag_cordinate_sys, level1
  use mo_file, only : file_mhm_output
 
  implicit none
 
contains
 
  !> \brief Initialize mHM OutputDataset
  !> \details Create and initialize the output file. If new a new output
  !!       variable needs to be written, this is the first of two
  !!       procedures to change (second: updateDataset)
  !!
  !> \changelog
  !! - Robert Schweppe Jun 2018
  !!   - refactoring and reformatting
  !! - Pallav Shrestha Mar 2020
  !!   - iFlag_cordinate_sys based dimensions (dims)
  !!
  !> \return type(OutputDataset)
  !> \authors Matthias Zink
  !> \date Apr 2013
  function mhm_outputdataset(iDomain, mask1) result(out)
    implicit none
 
    !> domain id
    integer(i4), intent(in) :: idomain
    !> L1 mask to reconstruct the data
    logical, pointer, intent(in), dimension(:, :) :: mask1
 
    type(outputdataset) :: out
 
    integer(i4) :: ii, nn, ncells
    character(3) :: dtype
    character(16), dimension(3) :: dims
    character(16) :: unit
    type(outputvariable), dimension(size(outputFlxState) * nSoilHorizons_mHM) :: tmpvars
 
    ncells = level1(idomain)%nCells
 
    out = outputdataset( &
      idomain=idomain, &
      level=level1, &
      file_name=file_mhm_output, &
      double_precision=output_double_precision, &
      outputs_frequence=timestep_model_outputs, &
      time_reference=output_time_reference &
    )
 
    dtype = data_dtype(output_double_precision)
    unit = fluxesunit(idomain)
    dims = data_dims()
 
    ii = 0
 
    if (outputflxstate(1)) then
      ii = ii + 1
      tmpvars(ii) = outputvariable(&
              out%nc, "interception", dtype, dims, ncells, mask1, output_deflate_level, .true.)
      call set_attributes(&
              tmpvars(ii)%nc, "canopy interception storage", "mm", output_double_precision)
    end if
 
    if (outputflxstate(2)) then
      ii = ii + 1
      tmpvars(ii) = outputvariable(&
              out%nc, "snowpack", dtype, dims, ncells, mask1, output_deflate_level, .true.)
      call set_attributes(tmpvars(ii)%nc, "depth of snowpack", "mm", output_double_precision)
    end if
 
    if (outputflxstate(3)) then
      do nn = 1, nsoilhorizons_mhm
        ii = ii + 1
        tmpvars(ii) = outputvariable(out%nc, "SWC_L" // trim(num2str(nn, '(i2.2)')), &
                dtype, dims, ncells, mask1, output_deflate_level, .true.)
        call set_attributes(tmpvars(ii)%nc, &
                'soil water content of soil layer' // trim(num2str(nn)), "mm", output_double_precision)
      end do
    end if
 
    if (outputflxstate(4)) then
      do nn = 1, nsoilhorizons_mhm
        ii = ii + 1
        tmpvars(ii) = outputvariable(out%nc, "SM_L" // trim(num2str(nn, '(i2.2)')), &
                dtype, dims, ncells, mask1, output_deflate_level, .true.)
        call set_attributes(tmpvars(ii)%nc, &
                'volumetric soil moisture of soil layer' // trim(num2str(nn)), "mm mm-1", output_double_precision)
      end do
    end if
 
    if (outputflxstate(5)) then
      ii = ii + 1
      tmpvars(ii) = outputvariable(&
              out%nc, "SM_Lall", dtype, dims, ncells, mask1, output_deflate_level, .true.)
      call set_attributes(&
              tmpvars(ii)%nc, "average soil moisture over all layers", "mm mm-1", output_double_precision)
    end if
 
    if (outputflxstate(6)) then
      ii = ii + 1
      tmpvars(ii) = outputvariable(&
              out%nc, "sealedSTW", dtype, dims, ncells, mask1, output_deflate_level, .true.)
      call set_attributes(&
              tmpvars(ii)%nc, "reservoir of sealed areas (sealedSTW)", "mm", output_double_precision)
    end if
 
    if (outputflxstate(7)) then
      ii = ii + 1
      tmpvars(ii) = outputvariable(&
              out%nc, "unsatSTW", dtype, dims, ncells, mask1, output_deflate_level, .true.)
      call set_attributes(&
              tmpvars(ii)%nc, "reservoir of unsaturated zone", "mm", output_double_precision)
    end if
 
    if (outputflxstate(8)) then
      ii = ii + 1
      tmpvars(ii) = outputvariable(&
              out%nc, "satSTW", dtype, dims, ncells, mask1, output_deflate_level, .true.)
      call set_attributes(&
              tmpvars(ii)%nc, "water level in groundwater reservoir", "mm", output_double_precision)
    end if
 
    if (outputflxstate(18)) then
      ii = ii + 1
      tmpvars(ii) = outputvariable(&
              out%nc, "neutrons", dtype, dims, ncells, mask1, output_deflate_level, .true.)
      call set_attributes(&
              tmpvars(ii)%nc, "ground albedo neutrons", "cph", output_double_precision)
    end if
 
    if (outputflxstate(9)) then
      ii = ii + 1
      tmpvars(ii) = outputvariable(&
              out%nc, "PET", dtype, dims, ncells, mask1, output_deflate_level)
      call set_attributes(&
              tmpvars(ii)%nc, "potential Evapotranspiration", trim(unit), output_double_precision)
    end if
 
    if (outputflxstate(10)) then
      ii = ii + 1
      tmpvars(ii) = outputvariable(&
              out%nc, "aET", dtype, dims, ncells, mask1, output_deflate_level)
      call set_attributes(&
              tmpvars(ii)%nc, "actual Evapotranspiration", trim(unit), output_double_precision)
    end if
 
    if (outputflxstate(11)) then
      ii = ii + 1
      tmpvars(ii) = outputvariable(&
              out%nc, "Q", dtype, dims, ncells, mask1, output_deflate_level)
      call set_attributes(&
              tmpvars(ii)%nc, "total runoff generated by every cell", trim(unit), output_double_precision)
    end if
 
    if (outputflxstate(12)) then
      ii = ii + 1
      tmpvars(ii) = outputvariable(&
              out%nc, "QD", dtype, dims, ncells, mask1, output_deflate_level)
      call set_attributes(tmpvars(ii)%nc, &
              "direct runoff generated by every cell (runoffSeal)", trim(unit), output_double_precision)
    end if
 
    if (outputflxstate(13)) then
      ii = ii + 1
      tmpvars(ii) = outputvariable(&
              out%nc, "QIf", dtype, dims, ncells, mask1, output_deflate_level)
      call set_attributes(tmpvars(ii)%nc, &
              "fast interflow generated by every cell (fastRunoff)", trim(unit), output_double_precision)
    end if
 
    if (outputflxstate(14)) then
      ii = ii + 1
      tmpvars(ii) = outputvariable(&
              out%nc, "QIs", dtype, dims, ncells, mask1, output_deflate_level)
      call set_attributes(tmpvars(ii)%nc, &
              "slow interflow generated by every cell (slowRunoff)", trim(unit), output_double_precision)
    end if
 
    if (outputflxstate(15)) then
      ii = ii + 1
      tmpvars(ii) = outputvariable(&
              out%nc, "QB", dtype, dims, ncells, mask1, output_deflate_level)
      call set_attributes(&
              tmpvars(ii)%nc, "baseflow generated by every cell", trim(unit), output_double_precision)
    end if
 
    if (outputflxstate(16)) then
      ii = ii + 1
      tmpvars(ii) = outputvariable(&
              out%nc, "recharge", dtype, dims, ncells, mask1, output_deflate_level)
      call set_attributes(&
              tmpvars(ii)%nc, "groundwater recharge", trim(unit), output_double_precision)
    end if
 
    if (outputflxstate(17)) then
      do nn = 1, nsoilhorizons_mhm
        ii = ii + 1
        tmpvars(ii) = outputvariable(&
                out%nc, "soil_infil_L" // trim(num2str(nn, '(i2.2)')), &
                dtype, dims, ncells, mask1, output_deflate_level)
        call set_attributes(tmpvars(ii)%nc, &
                "infiltration flux from soil layer" // trim(num2str(nn)), unit, output_double_precision)
      end do
    end if
 
    if (outputflxstate(19)) then
      do nn = 1, nsoilhorizons_mhm
        ii = ii + 1
        tmpvars(ii) = outputvariable(&
                out%nc, "aET_L" // trim(num2str(nn, '(i2.2)')), &
                dtype, dims, ncells, mask1, output_deflate_level)
        call set_attributes(tmpvars(ii)%nc, &
                'actual Evapotranspiration from soil layer' // trim(num2str(nn)), &
                "mm " // trim(unit), output_double_precision)
      end do
    end if
 
    if (outputflxstate(20)) then
      ii = ii + 1
      tmpvars(ii) = outputvariable(&
              out%nc, "preEffect", dtype, dims, ncells, mask1, output_deflate_level)
      call set_attributes(&
              tmpvars(ii)%nc, "effective precipitation", trim(unit), output_double_precision)
    end if
 
    if (outputflxstate(21)) then
      ii = ii + 1
      tmpvars(ii) = outputvariable(&
              out%nc, "Qsm", dtype, dims, ncells, mask1, output_deflate_level)
      call set_attributes(&
              tmpvars(ii)%nc, "Average liquid water generated from solid to liquid phase change in the snow", &
              trim(unit), output_double_precision)
    end if
 
    allocate(out%vars(ii))
    out%vars = tmpvars(1 : ii)
 
  end function mhm_outputdataset
 
  !> \brief Update all variables.
  !> \details Call the type bound procedure updateVariable for
  !> all output variables. If a new output
  !> variable needs to be written, this is the second
  !> of two procedures to change (first: newOutputDataset)
  !!
  !> \changelog
  !! - L. Samaniego et al. Dec  2013
  !!   - nullify pointer Matthias Zink,        Feb. 2014
  !!   - added aditional output: pet V. Prykhodk, J. Mai,  Nov. 2014
  !!   - adding new variable infilSoil
  !!   - case 16 David Schaefer      , Jun. 2015
  !!   - major rewrite
  !! - Robert Schweppe Jun 2018
  !!   - refactoring and reformatting
  !> \authors Matthias Zink
  !> \date Apr 2013
  subroutine mhm_updatedataset(nc_mhm, sidx, eidx, L1_fSealed, L1_fNotSealed, L1_inter, L1_snowPack, L1_soilMoist, &
                          L1_soilMoistSat, L1_sealSTW, L1_unsatSTW, L1_satSTW, L1_neutrons, L1_pet, L1_aETSoil, &
                          L1_aETCanopy, L1_aETSealed, L1_total_runoff, L1_runoffSeal, L1_fastRunoff, L1_slowRunoff, &
                          L1_baseflow, L1_percol, L1_infilSoil, L1_preEffect, L1_melt)
 
    implicit none
 
    class(outputdataset), intent(inout), target :: nc_mhm
 
    integer(i4), intent(in) :: sidx !< start index
    integer(i4), intent(in) :: eidx !< end index
    real(dp), intent(in), dimension(:) :: L1_fSealed
    real(dp), intent(in), dimension(:) :: L1_fNotSealed
    real(dp), intent(in), dimension(:) :: L1_inter
    real(dp), intent(in), dimension(:) :: L1_snowPack
    real(dp), intent(in), dimension(:, :) :: L1_soilMoist
    real(dp), intent(in), dimension(:, :) :: L1_soilMoistSat
    real(dp), intent(in), dimension(:) :: L1_sealSTW
    real(dp), intent(in), dimension(:) :: L1_unsatSTW
    real(dp), intent(in), dimension(:) :: L1_satSTW
    real(dp), intent(in), dimension(:) :: L1_neutrons
    real(dp), intent(in), dimension(:) :: L1_pet
    real(dp), intent(in), dimension(:, :) :: L1_aETSoil
    real(dp), intent(in), dimension(:) :: L1_aETCanopy
    real(dp), intent(in), dimension(:) :: L1_aETSealed
    real(dp), intent(in), dimension(:) :: L1_total_runoff
    real(dp), intent(in), dimension(:) :: L1_runoffSeal
    real(dp), intent(in), dimension(:) :: L1_fastRunoff
    real(dp), intent(in), dimension(:) :: L1_slowRunoff
    real(dp), intent(in), dimension(:) :: L1_baseflow
    real(dp), intent(in), dimension(:) :: L1_percol
    real(dp), intent(in), dimension(:, :) :: L1_infilSoil
    real(dp), intent(in), dimension(:) :: L1_preEffect
    real(dp), intent(in), dimension(:) :: L1_melt
 
    type(outputvariable), pointer, dimension(:) :: vars
    integer(i4) :: ii, nn
 
    ii = 0
    vars => nc_mhm%vars
 
    if (outputflxstate(1)) then
      ii = ii + 1
      call vars(ii)%updateVariable(l1_inter(sidx : eidx))
    end if
 
    if (outputflxstate(2)) then
      ii = ii + 1
      call vars(ii)%updateVariable(l1_snowpack(sidx : eidx))
    end if
 
    if (outputflxstate(3)) then
      do nn = 1, nsoilhorizons_mhm
        ii = ii + 1
        call vars(ii)%updateVariable(l1_soilmoist(sidx : eidx, nn))
      end do
    end if
 
    if (outputflxstate(4)) then
      do nn = 1, nsoilhorizons_mhm
        ii = ii + 1
        call vars(ii)%updateVariable(l1_soilmoist(sidx : eidx, nn) / l1_soilmoistsat(sidx : eidx, nn))
      end do
    end if
 
    if (outputflxstate(5)) then
      ii = ii + 1
      call vars(ii)%updateVariable(sum(l1_soilmoist(sidx : eidx, :), dim = 2) / sum(l1_soilmoistsat(sidx : eidx, :), dim = 2))
    end if
 
    if (outputflxstate(6)) then
      ii = ii + 1
      call vars(ii)%updateVariable(l1_sealstw(sidx : eidx))
    end if
 
    if (outputflxstate(7)) then
      ii = ii + 1
      call vars(ii)%updateVariable(l1_unsatstw(sidx : eidx))
    end if
 
    if (outputflxstate(8)) then
      ii = ii + 1
      call vars(ii)%updateVariable(l1_satstw(sidx : eidx))
    end if
 
    if (outputflxstate(18)) then
      ii = ii + 1
      call vars(ii)%updateVariable(l1_neutrons(sidx : eidx))
    end if
 
    if (outputflxstate(9)) then
      ii = ii + 1
      call vars(ii)%updateVariable(l1_pet(sidx : eidx))
    end if
 
    if (outputflxstate(10)) then
      ii = ii + 1
      call vars(ii)%updateVariable(sum(l1_aetsoil(sidx : eidx, :), dim = 2) * l1_fnotsealed(sidx : eidx) &
              + l1_aetcanopy(sidx : eidx) + l1_aetsealed(sidx : eidx) * l1_fsealed(sidx : eidx))
    end if
 
    if (outputflxstate(11)) then
      ii = ii + 1
      call vars(ii)%updateVariable(l1_total_runoff(sidx : eidx))
    end if
 
    if (outputflxstate(12)) then
      ii = ii + 1
      call vars(ii)%updateVariable(l1_runoffseal(sidx : eidx) * l1_fsealed(sidx : eidx))
    end if
 
    if (outputflxstate(13)) then
      ii = ii + 1
      call vars(ii)%updateVariable(l1_fastrunoff(sidx : eidx) * l1_fnotsealed(sidx : eidx))
    end if
 
    if (outputflxstate(14)) then
      ii = ii + 1
      call vars(ii)%updateVariable(l1_slowrunoff(sidx : eidx) * l1_fnotsealed(sidx : eidx))
    end if
 
    if (outputflxstate(15)) then
      ii = ii + 1
      call vars(ii)%updateVariable(l1_baseflow(sidx : eidx) * l1_fnotsealed(sidx : eidx))
    end if
 
    if (outputflxstate(16)) then
      ii = ii + 1
      call vars(ii)%updateVariable(l1_percol(sidx : eidx) * l1_fnotsealed(sidx : eidx))
    end if
 
    if (outputflxstate(17)) then
      do nn = 1, nsoilhorizons_mhm
        ii = ii + 1
        call vars(ii)%updateVariable(l1_infilsoil(sidx : eidx, nn) * l1_fnotsealed(sidx : eidx))
      end do
    end if
 
    if (outputflxstate(19)) then
      do nn = 1, nsoilhorizons_mhm
        ii = ii + 1
        call vars(ii)%updateVariable(l1_aetsoil(sidx : eidx, nn) * l1_fnotsealed(sidx : eidx))
      end do
    end if
 
    if (outputflxstate(20)) then
      ii = ii + 1
      call vars(ii)%updateVariable(l1_preeffect(sidx : eidx))
    end if
 
    if (outputflxstate(21)) then
      ii = ii + 1
      call vars(ii)%updateVariable(l1_melt(sidx : eidx))
    end if
 
  end subroutine mhm_updatedataset
 
  !> \brief Generate a unit string
  !> \details Generate the unit string for the output variable netcdf attribute based on modeling timestep
  !!
  !> \changelog
  !! - Robert Schweppe Jun 2018
  !!   - refactoring and reformatting
  !!
  !> \return character(16)
  !> \authors David Schaefer
  !> \date June 2015
  function fluxesunit(iDomain)
 
    use mo_common_mhm_mrm_variables, only : ntstepday, simper
 
    implicit none
 
    integer(i4), intent(in) :: idomain !< selected domain
 
    character(16) :: fluxesunit
    real(dp) :: ntsteps
 
    if (timestep * timestep_model_outputs .eq. 1) then
      fluxesunit = 'mm h-1'
    else if (timestep_model_outputs > 1) then
      fluxesunit = 'mm ' // trim(adjustl(num2str(timestep))) // 'h-1'
    else if (timestep_model_outputs .eq. 0) then
      ntsteps = (simper(idomain)%julEnd - simper(idomain)%julStart + 1) * ntstepday
      fluxesunit = 'mm ' // trim(adjustl(num2str(nint(ntsteps)))) // 'h-1'
    else if (timestep_model_outputs .eq. -1) then
      fluxesunit = 'mm d-1'
    else if (timestep_model_outputs .eq. -2) then
      fluxesunit = 'mm month-1'
    else if (timestep_model_outputs .eq. -3) then
      fluxesunit = 'mm a-1'
    else
      fluxesunit = ''
    end if
 
  end function fluxesunit
 
end module mo_write_fluxes_states