growing_degree_day.F90

Go to the documentation of this file.

module growing_degree_day
 
  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_calculate, growing_degree_day_initialize
  public :: modified_growing_degree_day_calculate
 
  real (c_float), allocatable  :: gdd_base(:)
  real (c_float), allocatable  :: gdd_max(:)
  integer (c_int), allocatable :: gdd_reset_date(:)
  type (t_netcdf4_file), pointer    :: pncfile
 
contains
 
  subroutine growing_degree_day_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)
        gdd_max( indx ) = gdd_max_l( landuse_index( indx ) )
      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)
        gdd_base( indx ) = gdd_base_l( landuse_index( indx ) )
      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_initialize
 
!--------------------------------------------------------------------------------------------------
 
  impure elemental subroutine growing_degree_day_calculate( gdd, tmean, order )
 
    ! [ ARGUMENTS ]
    real (c_float), intent(inout)       :: gdd
    real (c_float), intent(in)          :: tmean
    integer (c_int), intent(in)         :: order
 
    associate( doy_to_reset_gdd => gdd_reset_date( order ),         &
               gdd_max => gdd_max( order ),                         &
               gdd_base => gdd_base( order ) )
 
      if ( sim_dt%iDOY == doy_to_reset_gdd )  gdd = 0.0_c_float
      
      gdd = gdd + max(tmean, gdd_base) - gdd_base
 
    end associate
 
  end subroutine growing_degree_day_calculate
 
!--------------------------------------------------------------------------------------------------
 
  impure elemental subroutine modified_growing_degree_day_calculate( gdd, tmin, tmax, order )
 
  ! [ ARGUMENTS ]
  real (c_float), intent(inout)       :: gdd
  real (c_float), intent(in)          :: tmin
  real (c_float), intent(in)          :: tmax
  integer (c_int), intent(in)         :: order
 
  associate( doy_to_reset_gdd => gdd_reset_date( order ),         &
             gdd_max => gdd_max( order ),                         &
             gdd_base => gdd_base( order ) )
 
    if ( sim_dt%iDOY == doy_to_reset_gdd )  gdd = 0.0_c_float
    
    gdd = gdd + max((max(tmin, gdd_base) + min(tmax,gdd_max))/2.0_c_float, gdd_base) - gdd_base
 
  end associate
 
end subroutine modified_growing_degree_day_calculate
 
end module growing_degree_day