mo_soil_database.f90

Go to the documentation of this file.

!> \file mo_soil_database.f90
!> \brief \copybrief mo_soil_database
!> \details \copydetails mo_soil_database
 
!> \brief Generating soil database from input file.
!> \details This module provides the routines for generating the soil database for mHM from an ASCII input file.
!! One routine \e read_soil_LUT reads a soil LookUpTable, performs some consistency checks and returns an initial
!! soil database.
!! The second routine \e generate_soil_database calculates based on the initial one the proper soil database.
!> \authors Juliane Mai
!> \date Dec 2012
!> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
!! mHM is released under the LGPLv3+ license \license_note
!> \ingroup f_mpr
MODULE mo_soil_database
 
  ! This module to provide a soil database for mHM.
 
  ! Written  Juliane Mai, Dec 2012
 
  use mo_kind, only : i4, dp
  use mo_message, only: message, error_message
  use mo_string_utils, only : num2str
  use mo_os, only : check_path_isfile
 
  IMPLICIT NONE
 
  PRIVATE
 
  PUBLIC :: read_soil_lut             ! Reads the soil LUT file
  PUBLIC :: generate_soil_database    ! Generates the soil database
 
  ! ------------------------------------------------------------------
 
CONTAINS
 
  ! ------------------------------------------------------------------
 
  !    NAME
  !        read_soil_LUT
 
  !    PURPOSE
  !>       \brief Reads the soil LUT file.
 
  !>       \details Reads the soil LookUpTable file and checks for consistency.
 
  !    INTENT(IN)
  !>       \param[in] "character(len = *) :: filename" filename of the soil LUT
 
  !    HISTORY
  !>       \authors Juliane Mai
 
  !>       \date Dec 2012
 
  ! Modifications:
  ! Luis Samaniego   Nov 2013 - transform relation op. == -> .eq. etc
  ! Rohini Kumar     Mar 2016 - new variables for handling different soil databases
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  subroutine read_soil_lut(filename)
 
    use mo_common_constants, only : eps_dp, nodata_dp, nodata_i4
    use mo_mpr_constants, only : nlcover_class
    use mo_mpr_file, only : usoil_database
    use mo_mpr_global_variables, only : horizondepth_mhm, iflag_soildb, nsoilhorizons_mhm, nsoiltypes, soildb, &
                                        tillagedepth
 
    implicit none
 
    ! filename of the soil LUT
    character(len = *), intent(in) :: filename
 
    character(len = 10) :: dummy
 
    integer(i4) :: ios
 
    integer(i4) :: ii, jj, kk
 
    integer(i4) :: nr, nh, max_i
 
    real(dp) :: up, down
 
    real(dp) :: cly, snd, bd
 
    real(dp) :: dmin
 
    integer(i4) :: maxnumberhorizons
 
    integer(i4) :: minnumbertillhorizons, maxnumbertillhorizons
 
 
    SELECT CASE (iflag_soildb)
      ! classical mHM soil database
    CASE(0)
      ios = 0_i4
      !checking whether the file exists
      call check_path_isfile(path = filename, raise=.true.)
      open(usoil_database, file = filename, status = 'old', iostat = ios)
      read(usoil_database, *) dummy, nsoiltypes
      dummy = dummy // ''   ! only to avoid warning
 
      ! allocate space
      allocate(soildb%Id(nsoiltypes))
      allocate(soildb%nHorizons(nsoiltypes))
      allocate(soildb%nTillHorizons(nsoiltypes))
      allocate(soildb%RZdepth(nsoiltypes))
      ! initialize
      soildb%Id(:) = nodata_i4
      soildb%nHorizons(:) = nodata_i4
      soildb%nTillHorizons(:) = nodata_i4
      soildb%RZdepth(:) = nodata_dp
 
      ! check if max soil_type id is equal to nSoilTypes
      max_i = 0_i4
      ! initalise total rows to read
      nr = 0_i4
      read(usoil_database, *) dummy
      do while (.NOT. (ios .ne. 0))
        read(usoil_database, *, iostat = ios) ii, jj, up, down, cly, snd, bd
 
        ! Checks
        if(up .ge. down) then
          call error_message('read_soil_LUT: ERROR occurred: Mixed horizon depths in soil type', &
                  num2str(ii, '(I3)'), ' and horizon no.', num2str(jj, '(I3)'))
        end if
        if(cly .lt. 0.0_dp .OR. cly .gt. 100.0_dp .OR. &
                snd .lt. 0.0_dp .OR. snd .gt. 100.0_dp .OR. &
                bd  .lt. 0.0_dp .OR. bd  .gt.   5.0_dp) then
          call error_message('read_soil_LUT: ERROR occurred: Inappropriate soil properties in soil type', &
                  num2str(ii, '(I3)'), ' and horizon no.', num2str(jj, '(I3)'))
        end if
 
        ! initalise soil id
        if ( ii > size(soildb%Id) ) call error_message( &
          "ERROR: nSoil_Types (", num2str(size(soildb%Id)), &
          ") in soil_classdefinition.txt seems to be to low! Tried to read: ", num2str(ii) &
        )
        soildb%Id(ii) = ii
        soildb%nHorizons(ii) = jj
        if( anint(down, dp) .gt. soildb%RZdepth(ii) ) soildb%RZdepth(ii) = anint(down, dp)
 
        nr = nr + 1_i4
        max_i = max(max_i, ii)
 
      end do
 
      if ( max_i < nsoiltypes ) call error_message( &
        "ERROR: nSoil_Types (", num2str(size(soildb%Id)), &
        ") in soil_classdefinition.txt seems to be to high! Highest read value: ", num2str(max_i) &
      )
 
      ! initalise minimum root zone depth among all soil types
      dmin = minval(soildb%RZdepth(:), soildb%RZdepth(:) .gt. 0.0_dp)
 
      ! check the tillage depth...(if possible adjust it..)
      if(tillagedepth .gt. dmin) then
        call error_message('read_soil_LUT: ERROR occurred: ', raise=.false.)
        call error_message('    Tillage depth is greater than overall minimum total soil depth ', raise=.false.)
        call error_message('    So tillage depth should be at least', num2str(dmin, '(F7.2)'), raise=.false.)
        call error_message('    Please adjust!')
      end if
 
      ! insert a new tillage soil layer, only in those soil types, in which it is not present
      rewind(usoil_database)
      read(usoil_database, *) dummy
      read(usoil_database, *) dummy
 
      ! Last row is read twice so, read only upto (nR - 1)
      do ii = 1, nr - 1
        read(usoil_database, *) jj, nh, up, down, cly, snd, bd
        if(tillagedepth .gt. anint(up, dp) .and. tillagedepth .lt. anint(down, dp)) then
          soildb%nHorizons(jj) = soildb%nHorizons(jj) + 1_i4
        end if
 
        ! identify upto which soil horizon does the tillage depth goes in
        if(tillagedepth .gt. anint(up, dp) .and. tillagedepth .le. anint(down, dp)) then
          soildb%nTillHorizons(jj) = nh
        end if
      end do
 
      maxnumberhorizons = maxval(soildb%nHorizons(:))
 
      ! the variables of soilType are all allocated with maximal number of Horizons although not for each of the
      ! nSoilTypes the array will be used
      ! loops for array(i,j) should be: j=1, nHorizons(i), otherwise a nodata_dp will appear
      allocate(soildb%UD(nsoiltypes, maxnumberhorizons))
      allocate(soildb%LD(nsoiltypes, maxnumberhorizons))
      allocate(soildb%depth(nsoiltypes, maxnumberhorizons))
      allocate(soildb%clay(nsoiltypes, maxnumberhorizons))
      allocate(soildb%sand(nsoiltypes, maxnumberhorizons))
      allocate(soildb%dbM(nsoiltypes, maxnumberhorizons))
 
      soildb%UD(:, :) = nodata_dp
      soildb%LD(:, :) = nodata_dp
      soildb%depth(:, :) = nodata_dp
      soildb%clay(:, :) = nodata_dp
      soildb%sand(:, :) = nodata_dp
      soildb%dbM(:, :) = nodata_dp
 
      ! allocate space for other derived soil hydraulic properties ...
      minnumbertillhorizons = minval(soildb%nTillHorizons(:))
      maxnumbertillhorizons = maxval(soildb%nTillHorizons(:))
 
      allocate(soildb%thetaS_till(nsoiltypes, 1 : maxnumbertillhorizons, nlcover_class))
      allocate(soildb%thetaS(nsoiltypes, minnumbertillhorizons + 1 : maxnumberhorizons))
 
      allocate(soildb%thetaFC_till(nsoiltypes, 1 : maxnumbertillhorizons, nlcover_class))
      allocate(soildb%thetaFC(nsoiltypes, minnumbertillhorizons + 1 : maxnumberhorizons))
 
      allocate(soildb%thetaPW_till(nsoiltypes, 1 : maxnumbertillhorizons, nlcover_class))
      allocate(soildb%thetaPW(nsoiltypes, minnumbertillhorizons + 1 : maxnumberhorizons))
 
      allocate(soildb%Db(nsoiltypes, 1 : maxnumbertillhorizons, nlcover_class))
 
      allocate(soildb%Ks(nsoiltypes, 1 : maxnumberhorizons, nlcover_class))
 
      ! Initialize with default nodata_dp
      soildb%thetaS_till(:, :, :) = nodata_dp
      soildb%thetaS(:, :) = nodata_dp
 
      soildb%thetaFC_till(:, :, :) = nodata_dp
      soildb%thetaFC(:, :) = nodata_dp
 
      soildb%thetaPW_till(:, :, :) = nodata_dp
      soildb%thetaPW(:, :) = nodata_dp
 
      soildb%Db(:, :, :) = nodata_dp
      soildb%Ks(:, :, :) = nodata_dp
 
      ! Read again soil properties  from the data base
      rewind(usoil_database)
      read(usoil_database, *) dummy
      read(usoil_database, *) dummy
 
      ! Last row is read twice so, read only upto (nR - 1)
      kk = 0_i4
      do ii = 1, nr - 1_i4
        read(usoil_database, *) jj, nh, up, down, cly, snd, bd
 
        ! to avoid numerical errors in PTF
        if(cly .lt. 1.0_dp) cly = 1.0_dp
        if(snd .lt. 1.0_dp) snd = 1.0_dp
 
        ! Physical consistency
        if((cly + snd) .gt. 100.0_dp) then
          cly = cly / (cly + snd)
          snd = snd / (cly + snd)
        end if
 
        ! check for an extra tillage horizon (if not exists create a layer)
        if(tillagedepth .gt. anint(up, dp) .and. tillagedepth .lt. anint(down, dp)) then
 
          soildb%UD(jj, nh) = anint(up, dp)
          soildb%LD(jj, nh) = tillagedepth
          soildb%depth(jj, nh) = soildb%LD(jj, nh) - soildb%UD(jj, nh)
          soildb%clay(jj, nh) = cly
          soildb%sand(jj, nh) = snd
          soildb%dbM(jj, nh) = bd
 
          ! initalise the upper depth to new value... and the increment counter
          up = tillagedepth
          kk = 1_i4
        end if
 
        ! increment nH by one once it encounter the upper loop..
        if(kk .eq. 1_i4) nh = nh + 1_i4
 
        soildb%UD(jj, nh) = anint(up, dp)
        soildb%LD(jj, nh) = anint(down, dp)
        soildb%depth(jj, nh) = soildb%LD(jj, nh) - soildb%UD(jj, nh)
        soildb%clay(jj, nh) = cly
        soildb%sand(jj, nh) = snd
        soildb%dbM(jj, nh) = bd
 
        ! check for number of soil horizons...
        if(nh .gt. soildb%nHorizons(jj)) then
          call error_message('read_soil_LUT: ERROR occurred: ', raise=.false.)
          call error_message('    There is something wrong in allocating horizons in soil data base.', raise=.false.)
          call error_message('    Please check in code !')
        end if
 
        ! initalise the increment counter to zero
        if(nh .eq. soildb%nHorizons(jj)) kk = 0_i4
 
      end do
      close(usoil_database)
 
      ! soil database for the horizon specific case
    CASE(1)
 
      allocate(soildb%RZdepth(1))
      allocate(soildb%UD(1, 1))
      allocate(soildb%LD(1, 1))
      allocate(soildb%depth(1, 1))
      allocate(soildb%thetaS_till(1, 1, 1))
      allocate(soildb%thetaS(1, 1))
      allocate(soildb%thetaFC_till(1, 1, 1))
      allocate(soildb%thetaFC(1, 1))
      allocate(soildb%thetaPW_till(1, 1, 1))
      allocate(soildb%thetaPW(1, 1))
      allocate(soildb%Db(1, 1, 1))
      allocate(soildb%Ks(1, 1, 1))
 
 
 
      !checking whether the file exists
      call check_path_isfile(path = filename, raise=.true.)
      open(usoil_database, file = filename, status = 'old', action = 'read')
      read(usoil_database, *) dummy, nsoiltypes
      dummy = dummy // ''   ! only to avoid warning
      allocate(soildb%Id  (nsoiltypes))
      allocate(soildb%clay(nsoiltypes, 1))
      allocate(soildb%sand(nsoiltypes, 1))
      allocate(soildb%dbM (nsoiltypes, 1))
      soildb%clay(:, :) = nodata_dp
      soildb%sand(:, :) = nodata_dp
      soildb%dbM (:, :) = nodata_dp
      read(usoil_database, *) dummy
      do kk = 1, nsoiltypes
        read(usoil_database, *, iostat=ios) jj, cly, snd, bd
        if ( ios /= 0 ) call error_message( &
          "ERROR: nSoil_Types (", num2str(nsoiltypes), ") in soil_classdefinition_iFlag_soilDB_1.txt ", &
          "seems to be higher than available soil types!" &
        )
        ! to avoid numerical errors in PTF
        if(cly .lt. 1.0_dp) cly = 1.0_dp
        if(snd .lt. 1.0_dp) snd = 1.0_dp
        soildb%Id(kk) = jj
        soildb%clay(kk, 1) = cly
        soildb%sand(kk, 1) = snd
        soildb%dbM (kk, 1) = bd
      end do
      close(usoil_database)
 
      ! assign up to which horizon layer a soil is treated as tillage layer
      !   since our horizon information is well defined for modeling too
      !   this information is uniform across all soils/modeling cells
      ! for compatibility with iFlag_soilDB == 0, assign
      ! both nHorizons & tillage horizons
      allocate(soildb%nHorizons(1))
      allocate(soildb%nTillHorizons(1))
      soildb%nHorizons(:) = nsoilhorizons_mhm
      soildb%nTillHorizons(:) = -9
      do kk = 1, nsoilhorizons_mhm
        if(abs(horizondepth_mhm(kk) - tillagedepth) .lt. eps_dp) then
          soildb%nTillHorizons(1) = kk
        end if
      end do
 
      ! check
      if(soildb%nTillHorizons(1) .eq. -9) then
        ! rarely could happen *** since this is checked in reading of horizons depths only
        ! but is checked here for double confirmation
        call error_message('***ERROR: specification of tillage depths is not confirming', raise=.false.)
        call error_message('          with given depths of soil horizons to be modeled.')
      else
        call message()
        call message('Tillage layers: the tillage horizons are modelled ')
        call message('                upto mHM layers: ', trim(num2str(soildb%nTillHorizons(1))))
      end if
 
    CASE DEFAULT
      call error_message('***ERROR: iFlag_soilDB option given does not exist. Only 0 and 1 is taken at the moment.')
 
    END SELECT
 
  end subroutine read_soil_lut
 
  ! ------------------------------------------------------------------
 
  !    NAME
  !        generate_soil_database
 
  !    PURPOSE
  !>       \brief Generates soil database.
 
  !>       \details Calculates the proper soil database using the initialized soil database from read_soil_LUT.
 
  !    HISTORY
  !>       \authors Juliane Mai
 
  !>       \date Dec 2012
 
  ! Modifications:
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  subroutine generate_soil_database
 
    use mo_common_constants, only : nodata_dp, nodata_i4
    use mo_mpr_global_variables, only : horizondepth_mhm, iflag_soildb, nsoilhorizons_mhm, nsoiltypes, soildb
 
    implicit none
 
    integer(i4) :: ii, jj, kk
 
    real(dp) :: dmin
 
    real(dp) :: dpth_f, dpth_t
 
    integer(i4) :: layer_f, layer_t
 
    real(dp), parameter :: small = 0.000001_dp
 
    ! [mm]       soil depth accuracy
    real(dp), parameter :: soil_daccuracy = 0.5_dp
 
 
    SELECT CASE (iflag_soildb)
      ! classical mHM soil database
    CASE(0)
      ! initalise minimum root zone depth among all soil types
      dmin = minval(soildb%RZdepth(:), soildb%RZdepth(:) > 0.0_dp)
 
      ! check
      if (horizondepth_mhm(nsoilhorizons_mhm - 1) .ge. dmin) then
        call error_message('generate_soil_database: ERROR occurred: ', raise=.false.)
        call error_message('   The depth of soil Horizons provided for modelling is not appropriate', raise=.false.)
        call error_message('   The global minimum of total soil horizon depth among all soil type is ', num2str(dmin, '(F7.2)'), raise=.false.)
        call error_message('   Adjust your modeling soil horizon depth in this range', raise=.false.)
        call error_message('   OR Increase the soil depth in data base for only those soil types', raise=.false.)
        call error_message('   whose total depth is smaller than your given modeling depth.')
     end if
 
      ! allocate and initalise depth weight
      allocate(soildb%Wd(nsoiltypes, nsoilhorizons_mhm, maxval(soildb%nHorizons(:))))
      soildb%Wd(:, :, :) = 0.0_dp
 
      ! Process further to estimate weight of each horizons
      ! weightage according to soil depths
      do ii = 1, nsoiltypes
        soildb%Wd(ii, :, soildb%nHorizons(ii) + 1_i4 : maxval(soildb%nHorizons(:))) = nodata_dp
 
        ! initalise last horizon depth to model w.r.t to surface
        ! NOTE:: it is different for each soil
        horizondepth_mhm(nsoilhorizons_mhm) = soildb%RZdepth(ii)
 
        ! Estimate soil properties for each modeling layers
        do jj = 1, nsoilhorizons_mhm
 
          ! modeling depth ( **from --> to ** )
          ! take into account the depth accuracy [0.5mm, defined in module..]
          dpth_f = 0.0_dp
          if(jj .ne. 1_i4) dpth_f = horizondepth_mhm(jj - 1)
          dpth_t = horizondepth_mhm(jj) - soil_daccuracy
 
          ! identify to which layer of batabase this mHM  horizon lies
          layer_f = nodata_i4
          layer_t = nodata_i4
          do kk = 1, soildb%nHorizons(ii)
            if(dpth_f .ge. soildb%UD(ii, kk) .and. dpth_f .le. (soildb%LD(ii, kk) - soil_daccuracy)) layer_f = kk
            if(dpth_t .ge. soildb%UD(ii, kk) .and. dpth_t .le. (soildb%LD(ii, kk) - soil_daccuracy)) layer_t = kk
          end do
 
          ! Check
          if(layer_f .le. 0_i4 .or. layer_t .le. 0_i4) then
            call error_message('generate_soil_database: ERROR occurred: ', raise=.false.)
            call error_message('     Horizon depths to model do not lie in database for soil type', num2str(ii, '(I3)'), raise=.false.)
            call error_message('     Please check!')
          end if
          if(layer_f .gt. layer_t) then
            call error_message('generate_soil_database: ERROR occurred: ', raise=.false.)
            call error_message('     Something is wrong in assignment of modeling soil horizons or', raise=.false.)
            call error_message('     database of soil type ', num2str(ii, '(I3)'), raise=.false.)
            call error_message('     Please check!')
          end if
 
          ! iF modeling depth of a given horizon falls in a same soil layer
          if(layer_f .eq. layer_t) then
            soildb%Wd(ii, jj, layer_f) = 1.0_dp
 
            ! else estimate depth weightage...
          else
 
            ! for starting layer...
            soildb%Wd(ii, jj, layer_f) = soildb%LD(ii, layer_f) - dpth_f
 
            ! for ending layer
            soildb%Wd(ii, jj, layer_t) = (dpth_t + soil_daccuracy) - soildb%UD(ii, layer_t)
 
            ! other intermediate layers weight, if exit
            if(layer_t - layer_f .gt. 1_i4) then
              do kk = layer_f + 1, layer_t - 1
                soildb%Wd(ii, jj, kk) = soildb%LD(ii, kk) - soildb%UD(ii, kk)
              end do
            end if
 
            ! Estimate depth weightage
            if(jj .ne. 1_i4) then
              soildb%Wd(ii, jj, 1 : soildb%nHorizons(ii)) = soildb%Wd(ii, jj, 1 : soildb%nHorizons(ii)) / &
                      (horizondepth_mhm(jj) - horizondepth_mhm(jj - 1_i4))
            else
              soildb%Wd(ii, jj, 1 : soildb%nHorizons(ii)) = soildb%Wd(ii, jj, 1 : soildb%nHorizons(ii)) / &
                      horizondepth_mhm(jj)
            end if
 
            ! Check (small margin for numerical errors)
            if(sum(soildb%Wd(ii, jj, :), soildb%Wd(ii, jj, :) .gt. 0.0_dp) .le. 1.0_dp - small .or. &
                    sum(soildb%Wd(ii, jj, :), soildb%Wd(ii, jj, :) .gt. 0.0_dp) .ge. 1.0_dp + small) then
              call error_message('generate_soil_database: ERROR occurred: ', raise=.false.)
              call error_message('     Weight assigned for each soil horizons are not correct.', raise=.false.)
              call error_message('     Please check!')
            end if
 
          end if
 
        end do
      end do
      ! soil database for the horizon specific case
    CASE(1)
      ! right now nothing is done here
      ! *** reserved for future changes
      allocate(soildb%Wd(1,1,1))
 
    CASE DEFAULT
      call error_message('***ERROR: iFlag_soilDB option given does not exist. Only 0 and 1 is taken at the moment.')
    END SELECT
 
  end subroutine generate_soil_database
 
END MODULE mo_soil_database