growing_degree_day_baskerville_emin.F90

Go to the documentation of this file.

module growing_degree_day_baskerville_emin
 
  use iso_c_binding
 
  use constants_and_conversions
  use datetime                   , only : mmdd2doy
  use exceptions
  use netcdf4_support
  use parameters
  use simulation_datetime
  use fstring_list
  implicit none
 
  private
 
  public :: gdd_base, gdd_max, gdd_reset_date
  public :: growing_degree_day_be_calculate, growing_degree_day_be_initialize
 
  real (c_float), allocatable  :: gdd_base(:)
  real (c_float), allocatable  :: gdd_max(:)
  integer (c_int), allocatable :: gdd_reset_date(:)
 
contains
 
  subroutine growing_degree_day_be_initialize( is_cell_active, landuse_index )
 
    logical (c_bool), intent(in)     :: is_cell_active(:,:)
    integer (c_int), intent(in)      :: landuse_index(:)
 
    ! [ LOCALS ]
    integer (c_int)              :: status
    integer (c_int)              :: indx
    type (fstring_list_t)              :: parameter_list
    type (fstring_list_t)              :: gdd_reset_val_list
    character (len=32)                :: sbuf
    real (c_float), allocatable  :: gdd_base_l(:)
    real (c_float), allocatable  :: gdd_max_l(:)
    integer (c_int)              :: number_of_landuse_codes
    integer (c_int), allocatable :: landuse_code(:)
 
    allocate( gdd_base( count( is_cell_active ) ), stat=status )
    call assert( status == 0, "Problem allocating memory", __file__, __line__ )
 
    allocate( gdd_max( count( is_cell_active ) ), stat=status )
    call assert( status == 0, "Problem allocating memory", __file__, __line__ )
 
    !> create string list that allows for alternate heading identifiers for the landuse code
    call parameter_list%append("LU_Code")
    call parameter_list%append("Landuse_Code")
    call parameter_list%append("Landuse_Lookup_Code")
 
    !> Determine how many landuse codes are present
    call params%get_parameters( slkeys=parameter_list, ivalues=landuse_code )
    number_of_landuse_codes = count( landuse_code >= 0 )
    call parameter_list%clear()
 
    call parameter_list%append("GDD_Base_Temp")
    call parameter_list%append("GDD_Base_Temperature")
    call parameter_list%append("GDD_Base")
 
    call params%get_parameters( slkeys=parameter_list, fvalues=gdd_base_l )
    call parameter_list%clear()
 
    call parameter_list%append("GDD_Max_Temp")
    call parameter_list%append("GDD_Maximum_Temperature")
    call parameter_list%append("GDD_Maximum_Temp")
    call parameter_list%append("GDD_Max")
 
    call params%get_parameters( slkeys=parameter_list, fvalues=gdd_max_l )
    call parameter_list%clear()
 
    call parameter_list%append("GDD_Reset_Date")
    call parameter_list%append("GDD_Reset")
 
    call params%get_parameters( slkeys=parameter_list, slvalues=gdd_reset_val_list )
    call parameter_list%clear()
 
    allocate( gdd_reset_date( count( is_cell_active ) ), stat=status )
    call assert( status==0, "Problem allocating memory.", __file__, __line__ )
 
    if ( gdd_reset_val_list%count == number_of_landuse_codes      &
         .and. gdd_reset_val_list%count_matching("<NA>") == 0 ) then
 
      ! retrieve gdd reset values; convert mm/dd to DOY
      do indx=1, gdd_reset_val_list%count
        sbuf = gdd_reset_val_list%get( indx )
 
        where ( landuse_index == indx )
          gdd_reset_date = mmdd2doy( sbuf, "GDD_RESET_DATE" )
        end where
 
      enddo
 
    else
 
      ! if no GDD_RESET_DATE found in parameter tables, assign the default: reset GDD at end of calendar year
      gdd_reset_date = 365_c_int
 
    endif
 
 
    if ( ubound( gdd_max_l, 1 ) == number_of_landuse_codes        &
        .and. gdd_max_l(1) > rtinyval ) then
 
      do indx=1, ubound( landuse_index, 1)
        if( landuse_index( indx ) >= lbound( gdd_max, 1) .and. landuse_index( indx ) <= ubound( gdd_max, 1) ) then
          gdd_max( indx ) = gdd_max_l( landuse_index( indx ) )
        endif
      enddo
 
    else
 
      ! if no GDD_MAX found in parameter tables, assign the default FAO-56 value
      gdd_max = 86.0_c_float
 
    endif
 
    if ( ubound( gdd_base_l, 1 ) == number_of_landuse_codes        &
        .and. gdd_base_l(1) > rtinyval  ) then
 
      do indx=1, ubound( landuse_index, 1)
        if( landuse_index( indx ) >= lbound( gdd_base, 1) .and. landuse_index( indx ) <= ubound( gdd_base, 1) ) then
          gdd_base( indx ) = gdd_base_l( landuse_index( indx ) )
        endif
      enddo
 
    else
 
      ! if no GDD_Base found in parameter tables, assign the default FAO-56 value
      gdd_base = 50.0_c_float
 
    endif
 
  end subroutine growing_degree_day_be_initialize
 
!--------------------------------------------------------------------------------------------------
 
  subroutine growing_degree_day_be_calculate( gdd, tmean, tmin, tmax, order )
 
    ! [ ARGUMENTS ]
    real (c_float), intent(inout)       :: gdd(:)
    real (c_float), intent(in)          :: tmean(:)
    real (c_float), intent(in)          :: tmin(:)
    real (c_float), intent(in)          :: tmax(:)
    integer (c_int), intent(in)         :: order(:)
 
    ! [ LOCALS ]
    real (c_float)    :: tmax_l
    real (c_float)    :: dd
    real (c_float)    :: w
    real (c_float)    :: at
    real (c_float)    :: a
    integer (c_int)   :: indx
 
    do indx=lbound(order,1),ubound(order,1)
 
      if ( sim_dt%iDOY == gdd_reset_date( order(indx) ) )  gdd(indx) = 0.0_c_float
 
      tmax_l = min( tmax(indx), gdd_max( order(indx) ) )
 
      if ( tmax_l <= gdd_base( order(indx) ) ) then
 
        dd = 0.0_c_float
 
      elseif ( tmin(indx) >= gdd_base( order(indx) ) ) then
 
        dd = min(( tmean(indx) - gdd_base( order(indx) )), (gdd_max( order(indx) ) - gdd_base( order(indx) ) ) )
 
      else
 
        w = ( tmax_l -  tmin(indx)) / 2.0_c_float
 
        at = ( gdd_base( order(indx) ) - tmean(indx) ) / w
 
        if ( at > 1 ) at = 1.0_c_float
        if ( at < -1 ) at = -1._c_float
 
        a = asin( at )
 
        dd = real(( ( w * cos( a ) ) - ( ( gdd_base( order(indx) ) - tmean(indx) )                    &
               * ( real( pi / 2._c_double, c_float ) - a ) ) ) / pi, c_float)
 
      endif
 
      gdd(indx) = gdd(indx) + dd
 
    enddo  
 
  end subroutine growing_degree_day_be_calculate
 
!--------------------------------------------------------------------------------------------------
 
end module growing_degree_day_baskerville_emin