awc__depth_integrated.F90

Go to the documentation of this file.

!> @file
!>  Contains the module awc__depth_integrated, which populates
!! the available water content by reading in and depth-averaging
!! soil available water contents over multiple soil horizons.
 
 
!> Populate the available water content by reading in and depth-averaging
!! soil available water contents over multiple soil horizons.
 
module awc__depth_integrated
 
  use iso_c_binding, only : c_short, c_int, c_float, c_double
  use constants_and_conversions
  use data_catalog
  use data_catalog_entry
  use exceptions
  use parameters
  use fstring
  use fstring_list
 
  implicit none
 
  private
 
  public :: awc_depth_integrated_read, awc_depth_integrated_initialize, available_water_content
 
  type (data_catalog_entry_t), pointer :: psoils_code_grid
  real (c_float), allocatable     :: available_water_content(:,:)
 
contains
 
  subroutine awc_depth_integrated_read( fRooting_Depth )
 
    real (c_float), intent(inout)       :: frooting_depth(:,:)
 
    ! [ LOCALS ]
 
    integer (c_int), allocatable  :: ilanduse_code(:)
    integer (c_int), allocatable  :: isoils_table_code(:)
    integer (c_int), allocatable  :: isoils_components(:)
    integer (c_int), allocatable  :: isoils_horizons(:)    
    real (c_float), allocatable   :: fsoils_awc(:) 
    real (c_float), allocatable   :: fsoils_top_depth(:)
    real (c_float), allocatable   :: fsoils_bottom_depth(:)
    real (c_float), allocatable   :: fsoils_component_fraction(:)
    real (c_float), allocatable   :: fsoils_horizon_thickness(:)
                   
    integer (c_int)               :: inumberoflanduses
    integer (c_int)               :: inumberofsoils
    integer (c_int)               :: inumberofsoilscomponents
    integer (c_int)               :: inumberofsoilshorizons
    integer (c_int)               :: inumberofsoilsawc
    integer (c_int)               :: inumberofsoilstopdepths
    integer (c_int)               :: inumberofsoilsbottomdepths
    integer (c_int)               :: inumberofsoilscomponentfractions     
 
    real (c_float)                :: ftemp_awc               
    real (c_float)                :: fdepthofdeepesthorizon
    real (c_float)                :: ffinal_awc
    integer (c_int)               :: ideepestsoilhorizon
    logical (c_bool)              :: lfirst
 
    type (fstring_list_t)               :: sllist
    integer (c_int)               :: istat
    integer (c_int)               :: iindex, iindex2
    integer (c_int)               :: iindex_x, iindex_y
    real (c_float)                :: frooting_depth_inches
    real (c_float)                :: fsoil_thickness_total
 
   call sllist%append("LU_Code")
   call sllist%append("Landuse_Code")
   call sllist%append("Landuse_Lookup_Code")
 
   call params%get_parameters( slkeys=sllist, ivalues=ilanduse_code )
   inumberoflanduses = count( ilanduse_code >= 0 )
   call sllist%clear()
 
 
   call sllist%append("Soils_Code")
   call sllist%append("Soils_Lookup_Code")
   call sllist%append("Soil_Code")
   call sllist%append("Soil_Lookup_Code")
 
   call params%get_parameters( slkeys=sllist, ivalues=isoils_table_code )
   inumberofsoils = count( isoils_table_code >= 0 )
   call sllist%clear()
 
   call sllist%append("Soils_Component")
   call sllist%append("Soils_Component_Number")
   call sllist%append("Soil_Component")
   call sllist%append("Soil_Component_Number")
 
   call params%get_parameters( slkeys=sllist, ivalues=isoils_components )
   inumberofsoilscomponents = count( isoils_components >= 0 )
   call sllist%clear()
 
   call sllist%append("Soils_Top_Depth")
   call sllist%append("Soil_Top_Depth")   
   call sllist%append("Soils_Z_Top")
   call sllist%append("Soils_Top_of_Horizon")
 
   call params%get_parameters( slkeys=sllist, fvalues=fsoils_top_depth )
   inumberofsoilstopdepths = count( fsoils_top_depth >= 0 )
   call sllist%clear()
 
   call sllist%append("Soils_Bottom_Depth")
   call sllist%append("Soil_Bottom_Depth")   
   call sllist%append("Soils_Z_Bottom")
   call sllist%append("Soils_Bottom_of_Horizon")
 
   call params%get_parameters( slkeys=sllist, fvalues=fsoils_bottom_depth )
   inumberofsoilsbottomdepths = count( fsoils_bottom_depth >= 0 )
   call sllist%clear()
 
   call sllist%append("Soils_Horizon")
   call sllist%append("Soils_Horizon_Number")
   call sllist%append("Soil_Horizon")
   call sllist%append("Soil_Horizon_Number")
 
   call params%get_parameters( slkeys=sllist, ivalues=isoils_horizons )
   inumberofsoilshorizons = count( isoils_horizons >= 0 )
   call sllist%clear()
 
   call sllist%append("Soils_Component_Fraction")
   call sllist%append("Soil_Component_Fraction")
 
   call params%get_parameters( slkeys=sllist, fvalues=fsoils_component_fraction )
   inumberofsoilscomponentfractions = count( fsoils_component_fraction >= 0 )
   call sllist%clear()
 
 
   call sllist%append("Soils_Available_Water_Content")
   call sllist%append("Soils_AWC")
   call sllist%append("Soil_Available_Water_Content")
   call sllist%append("Soil_AWC")
   call sllist%append("Available_Water_Content") 
   call sllist%append("AWC")
 
   call params%get_parameters( slkeys=sllist, fvalues=fsoils_awc )
   inumberofsoilsawc = count( fsoils_awc >= 0 )
   call sllist%clear()
 
   !> @todo Need a *LOT* more data validation and guard code here to prevent execution in the case
   !! where only some data fields are found, or are only partially populated.
 
    ! locate the data structure associated with the gridded rainfall zone entries
    psoils_code_grid => dat%find("SOILS_CODE")
    call assert( associated(psoils_code_grid), &
      "A SOILS_CODE grid must be supplied in order to make use of this option.", __file__, __line__)
 
    call psoils_code_grid%getvalues( )
 
    allocate (fsoils_horizon_thickness( ubound( fsoils_bottom_depth, 1) ), stat=istat )
 
    allocate(available_water_content( ubound(frooting_depth, 1), ubound(frooting_depth, 2) ), stat=istat )
 
    ! this should be in units of inches
    fsoils_horizon_thickness = fsoils_bottom_depth - fsoils_top_depth
 
    do iindex_x=lbound( frooting_depth, 1 ), ubound( frooting_depth, 1 )
 
      do iindex_y=lbound( frooting_depth, 2 ), ubound( frooting_depth, 2 )
 
        ! rooting depth in feet * 12 in/ft = rooting depth in inches
        frooting_depth_inches = frooting_depth( iindex_x, iindex_y ) * 12.0_c_float
 
        ftemp_awc = 0.0_c_float
        lfirst = true
        ideepestsoilhorizon = 0_c_int
        fdepthofdeepesthorizon = 0.0_c_float
        ffinal_awc = 0.0_c_float
 
        do iindex=1, inumberofsoils
          if ( psoils_code_grid%pGrdBase%iData( iindex_x, iindex_y ) == isoils_table_code( iindex ) ) then
 
            ! look for current soil code; calculate the component-weighted mean AWC for the deepest horizon
            if ( lfirst ) then
 
              lfirst = false
 
              ! find deepest soil horizon for the current soils code
              do iindex2=iindex, inumberofsoils
                if ( psoils_code_grid%pGrdBase%iData( iindex_x, iindex_y ) == isoils_table_code( iindex2 ) )  &
                  ideepestsoilhorizon = max( ideepestsoilhorizon, isoils_horizons( iindex2 ) )
              enddo
 
              ! if we are in the deepest horizon of the current soils code, calculate the composite averaged AWC
              do iindex2=iindex, inumberofsoils
                if ( ( psoils_code_grid%pGrdBase%iData( iindex_x, iindex_y ) == isoils_table_code( iindex2 ) )  &
                  .and. ( ideepestsoilhorizon == isoils_horizons( iindex2) ) )   &
                  ffinal_awc = ffinal_awc + fsoils_awc( iindex2 ) * fsoils_component_fraction( iindex2 )
                  fdepthofdeepesthorizon = fsoils_bottom_depth( iindex2 )
              enddo
 
            endif  
 
            ! rooting depth never reaches to this horizon; ignore
            if ( frooting_depth_inches < fsoils_top_depth( iindex ) ) cycle
 
            ! rooting depth exceeds lower bound of current horizon; entire horizon
            ! depth used as weight factor
            if ( frooting_depth_inches > fsoils_bottom_depth( iindex ) ) then
 
              ftemp_awc = ftemp_awc + fsoils_awc( iindex ) * fsoils_component_fraction( iindex )     &
                                    * fsoils_horizon_thickness( iindex )
 
            ! rooting depth falls in between the upper and lower bounds of the current soil horizon
            elseif ( frooting_depth_inches >= fsoils_top_depth( iindex )                    &
                   .and. frooting_depth_inches <= fsoils_bottom_depth( iindex ) ) then
            
              ftemp_awc = ftemp_awc + fsoils_awc( iindex ) * fsoils_component_fraction( iindex )     &
                                    * ( frooting_depth_inches - fsoils_top_depth( iindex ) )
 
            endif                         
 
          endif
 
        enddo  
 
        ! if the soil data from the table does not extend to the rooting depth, extrapolate
        if (frooting_depth_inches > fdepthofdeepesthorizon )  &
 
              ftemp_awc = ftemp_awc + ffinal_awc * ( frooting_depth_inches - fdepthofdeepesthorizon )
 
 
        if ( frooting_depth_inches > 0.0_c_float ) then
          ftemp_awc = ftemp_awc / frooting_depth_inches  ! divide by total rooting depth to obtain weighted average
        else
          ftemp_awc = 0.0_c_float
        endif    
 
        ! table values are in/in; multiply by 12 in/ft to return AWC in inches/foot
        available_water_content(iindex_x, iindex_y) = ftemp_awc * 12.0_c_float 
 
      enddo
 
    enddo
 
 
! Original HWB code...
 
! c........compute area-weighted soil-moisture storage capacity and account for
! c        depth-varying AWC
!          smca(i,k)=0.
!          do 420 j=1,nseq(i)
!          do 400 l=1,nlay(i,j)
!             if(z.ge.zt(i,j,l).and.z.le.zb(i,j,l))then
! c...........bottom of root zone is within current layer
!                smca(i,k)=smca(i,k)+pct(i,j,l)*awc(i,j,l)*(z-zt(i,j,l))
!             elseif(z.gt.zb(i,j,l).and.l.eq.nlay(i,j))then
! c...........bottom of root zone is below current layer, which is the
! c           bottom layer; extrapolate bottom layer down to root depth
!                smca(i,k)=smca(i,k)+pct(i,j,l)*awc(i,j,l)*(z-zt(i,j,l))
!             elseif(z.gt.zb(i,j,l).and.l.lt.nlay(i,j))then
! c...........bottom of root zone is below current layer, which is not the
! c           bottom layer; include entire current layer
!                smca(i,k)=smca(i,k)+pct(i,j,l)*awc(i,j,l)*
!      1         (zb(i,j,l)-zt(i,j,l))
! c...........exit if bottom of root zone is above current layer
!             elseif(z.lt.zt(i,j,l))then
!                goto 440
!             endif
!  400     continue
!  420     continue
!  440     continue
 
 
 
 
  end subroutine awc_depth_integrated_read
 
!--------------------------------------------------------------------------------------------------
 
  subroutine awc_depth_integrated_initialize( lActive, fAWC, iSoils_Code )
 
    logical (c_bool), intent(in)     :: lactive(:,:)
    real (c_float), intent(inout)    :: fawc(:)
    integer (c_int), intent(inout)   :: isoils_code(:)
 
    isoils_code = pack( psoils_code_grid%pGrdBase%iData, lactive )
 
    fawc = pack( available_water_content, lactive )
 
  end subroutine awc_depth_integrated_initialize
 
end module awc__depth_integrated