crop_coefficients__fao56.F90

Go to the documentation of this file.

!> @file
!>  Contains a single module, \ref crop_coefficients__fao56, which
!>  provides support for modifying reference ET through the use of
!>  crop coefficients
 
!> Update crop coefficients for crop types in simulation.
 
module crop_coefficients__fao56
 
  use iso_c_binding, only             : c_bool, c_short, c_int, c_float, c_double
  use constants_and_conversions, only : m_per_foot, true, false, ftinyval,       &
                                        itinyval, asint, asfloat, fzero, in_to_mm, &
                                        true, false, clip, near_zero
  use data_catalog, only              : dat
  use data_catalog_entry, only        : data_catalog_entry_t
  use datetime
  use logfiles, only                  : logs, log_all
  use exceptions, only                : assert, warn, die
  use parameters, only                : params
  use simulation_datetime, only       : sim_dt
  use fstring, only                   : ascharacter, squote, operator(.contains.)
  use fstring_list
  implicit none
 
  private
 
  public :: crop_coefficients_fao56_initialize, crop_coefficients_fao56_calculate
  public :: crop_coefficients_fao56_update_growth_stage_dates
  public :: crop_coefficients_fao56_update_growing_season
  public :: crop_coefficients_fao56_calculate_kcb_max
  public :: update_crop_coefficient_date_as_threshold, update_crop_coefficient_gdd_as_threshold
  public :: growth_stage_date, planting_date, growth_stage_length_in_days
  public :: kcb_min, kcb_ini, kcb_mid, kcb_end
  public :: kcb_l, jan, dec, kcb_method, kcb_method_gdd, kcb_method_fao56
  public :: kcb_method_monthly_values
 
  enum, bind(c)
    enumerator :: l_doy_ini=1, l_doy_dev, l_doy_mid, l_doy_late, l_doy_fallow
  end enum
 
  enum, bind(c)
    enumerator :: gdd_plant=1, gdd_ini, gdd_dev, gdd_mid, gdd_late
  end enum
 
  enum, bind(c)
    enumerator :: planting_date=1, enddate_ini, enddate_dev, enddate_mid, enddate_late, &
                    enddate_fallow
  end enum
 
  enum, bind(c)
    enumerator :: kcb_ini=13, kcb_mid, kcb_end, kcb_min
  end enum
 
  enum, bind(c)
    enumerator :: jan = 1, feb, mar, apr, may, jun, jul, aug, sep, oct, nov, dec
  end enum
 
   enum, bind(c)
     enumerator :: kcb_method_gdd = 1, kcb_method_monthly_values, kcb_method_fao56
   end enum
 
  ! Private, module level variables
  ! kept at a landuse code level (i.e. same value applies to all cells with same LU codes)
  integer (c_int), allocatable  :: landuse_code(:)
!  real (c_float), allocatable   :: REW(:,:)
!  real (c_float), allocatable   :: TEW(:,:)
  real (c_float), allocatable   :: kcb_l(:,:)
  integer (c_int), allocatable  :: kcb_method(:)
  real (c_float), allocatable   :: growth_stage_length_in_days(:,:)
  real (c_float), allocatable   :: growth_stage_gdd(:,:)
  type (datetime_t), allocatable     :: growth_stage_date(:,:)
  type (fstring_list_t)         :: sl_planting_date
 
contains
 
 
  function crop_coefficients_calculate_planting_date(sPlantingDate, iYear)   result(dtNewPlantingDate)
    character(len=*), intent(in)      :: splantingdate
    integer (c_int), intent(in)       :: iyear
    type (datetime_t)                 :: dtnewplantingdate
 
    ! [ LOCAL ]
    character (len=10)                :: smmddyyyy
 
    if ( splantingdate .contains. "/" ) then
      smmddyyyy = trim(splantingdate)//"/"//ascharacter( iyear )
      call dtnewplantingdate%parsedate( smmddyyyy, __file__, __line__ )
 
    else
      ! assume the value is a day-of-year value
      smmddyyyy = "01/01/"//ascharacter( iyear )
      call dtnewplantingdate%parsedate( smmddyyyy, __file__, __line__ )
      dtnewplantingdate = dtnewplantingdate + ( asfloat( splantingdate ) - 1.0_c_float )
      call dtnewplantingdate%calcJulianDay()
    endif
 
  end function crop_coefficients_calculate_planting_date
 
 
  subroutine crop_coefficients_fao56_initialize()
 
    ! [ LOCALS ]
    ! type (FSTRING_LIST_T)            :: slREW, slTEW
    type (fstring_list_t)             :: sllist
    ! integer (c_int), allocatable :: iTEWSeqNums(:)
    ! integer (c_int), allocatable :: iREWSeqNums(:)
    integer (c_int)                   :: inumberoflanduses
    integer (c_int)                   :: iindex
    integer (c_int)                   :: istat
 
 
    !type (FSTRING_LIST_T)             :: slPlantingDate
    character (len=:), allocatable    :: plantingdate_str
 
    real (c_float), allocatable       :: l_ini_l(:)
    real (c_float), allocatable       :: l_dev_l(:)
    real (c_float), allocatable       :: l_mid_l(:)
    real (c_float), allocatable       :: l_late_l(:)
    real (c_float), allocatable       :: l_fallow_l(:)
 
    real (c_float), allocatable       :: gdd_plant_l(:)
    real (c_float), allocatable       :: gdd_ini_l(:)
    real (c_float), allocatable       :: gdd_dev_l(:)
    real (c_float), allocatable       :: gdd_mid_l(:)
    real (c_float), allocatable       :: gdd_late_l(:)
 
 
    real (c_float), allocatable       :: kcb_ini_l(:)
    real (c_float), allocatable       :: kcb_mid_l(:)
    real (c_float), allocatable       :: kcb_end_l(:)
    real (c_float), allocatable       :: kcb_min_l(:)
 
    real (c_float), allocatable       :: kcb_jan(:)
    real (c_float), allocatable       :: kcb_feb(:)
    real (c_float), allocatable       :: kcb_mar(:)
    real (c_float), allocatable       :: kcb_apr(:)
    real (c_float), allocatable       :: kcb_may(:)
    real (c_float), allocatable       :: kcb_jun(:)
    real (c_float), allocatable       :: kcb_jul(:)
    real (c_float), allocatable       :: kcb_aug(:)
    real (c_float), allocatable       :: kcb_sep(:)
    real (c_float), allocatable       :: kcb_oct(:)
    real (c_float), allocatable       :: kcb_nov(:)
    real (c_float), allocatable       :: kcb_dec(:)
 
 
 
   !> create string list that allows for alternate heading identifiers for the landuse code
   sllist = create_list("LU_Code, Landuse_Code, Landuse_Lookup_Code")
 
   !> Determine how many landuse codes are present
   call params%get_parameters( slkeys=sllist, ivalues=landuse_code )
   inumberoflanduses = count( landuse_code >= 0 )
   !> @todo Implement thorough input error checking:
   !! are all soils in grid included in table values?
   !> is soil suffix vector continuous?
 
   ! Retrieve and populate the Readily Evaporable Water (REW) table values
  !  CALL PARAMS%get_parameters( fValues=REW, sPrefix="REW_", iNumRows=iNumberOfLanduses )
 
   ! Retrieve and populate the Total Evaporable Water (TEW) table values
  !  CALL PARAMS%get_parameters( fValues=TEW, sPrefix="TEW_", iNumRows=iNumberOfLanduses )
 
   !> @TODO What should happen if the TEW / REW header entries do *not* fall in a
   !!       logical sequence of values? In other words, if the user has columns named
   !!       REW_1, REW_3, REW_5, only the values associated with "REW_1" would be retrieved.
   !!       Needless to say, this would be catastrophic.
 
   call params%get_parameters( skey="Planting_date", slvalues=sl_planting_date )
 
   call params%get_parameters( skey="L_ini", fvalues=l_ini_l)
   call params%get_parameters( skey="L_dev", fvalues=l_dev_l)
   call params%get_parameters( skey="L_mid", fvalues=l_mid_l)
   call params%get_parameters( skey="L_late", fvalues=l_late_l)
   call params%get_parameters( skey="L_fallow", fvalues=l_fallow_l)
 
   call params%get_parameters( skey="GDD_plant", fvalues=gdd_plant_l)
   call params%get_parameters( skey="GDD_ini", fvalues=gdd_ini_l)
   call params%get_parameters( skey="GDD_dev", fvalues=gdd_dev_l)
   call params%get_parameters( skey="GDD_mid", fvalues=gdd_mid_l)
   call params%get_parameters( skey="GDD_late", fvalues=gdd_late_l)
 
   call params%get_parameters( skey="Kcb_ini", fvalues=kcb_ini_l)
   call params%get_parameters( skey="Kcb_mid", fvalues=kcb_mid_l)
   call params%get_parameters( skey="Kcb_end", fvalues=kcb_end_l)
   call params%get_parameters( skey="Kcb_min", fvalues=kcb_min_l)
 
   call params%get_parameters( skey="Kcb_Jan", fvalues=kcb_jan )
   call params%get_parameters( skey="Kcb_Feb", fvalues=kcb_feb )
   call params%get_parameters( skey="Kcb_Mar", fvalues=kcb_mar )
   call params%get_parameters( skey="Kcb_Apr", fvalues=kcb_apr )
   call params%get_parameters( skey="Kcb_May", fvalues=kcb_may )
   call params%get_parameters( skey="Kcb_Jun", fvalues=kcb_jun )
   call params%get_parameters( skey="Kcb_Jul", fvalues=kcb_jul )
   call params%get_parameters( skey="Kcb_Aug", fvalues=kcb_aug )
   call params%get_parameters( skey="Kcb_Sep", fvalues=kcb_sep )
   call params%get_parameters( skey="Kcb_Oct", fvalues=kcb_oct )
   call params%get_parameters( skey="Kcb_Nov", fvalues=kcb_nov )
   call params%get_parameters( skey="Kcb_Dec", fvalues=kcb_dec )
 
 
    allocate( growth_stage_length_in_days( 5, inumberoflanduses ), stat=istat )
    call assert( istat==0, "Failed to allocate memory for GROWTH_STAGE_LENGTH_IN_DAYS array", &
      __file__, __line__ )
 
    allocate( growth_stage_gdd( 5, inumberoflanduses ), stat=istat )
    call assert( istat==0, "Failed to allocate memory for GROWTH_STAGE_GDD array", &
      __file__, __line__ )
 
    allocate( growth_stage_date( 6, inumberoflanduses ), stat=istat )
    call assert( istat==0, "Failed to allocate memory for GROWTH_STAGE_DATE array", &
      __file__, __line__ )
 
    allocate( kcb_l( 16, inumberoflanduses ), stat=istat )
    call assert( istat==0, "Failed to allocate memory for KCB_l array", &
      __file__, __line__ )
 
    allocate( kcb_method( inumberoflanduses ), stat=istat )
    call assert( istat==0, "Failed to allocate memory for KCB_METHOD vector", &
      __file__, __line__ )
 
    kcb_method = -9999
    kcb_l = -9999.
    growth_stage_gdd = -9999.
    growth_stage_length_in_days = 0.
 
    if ( ubound(l_ini_l,1) == inumberoflanduses )           &
      growth_stage_length_in_days( l_doy_ini,  : ) = l_ini_l
 
    if ( ubound(l_dev_l,1) == inumberoflanduses )           &
      growth_stage_length_in_days( l_doy_dev,  : ) = l_dev_l
 
    if ( ubound(l_mid_l,1) == inumberoflanduses )           &
      growth_stage_length_in_days( l_doy_mid,  : ) = l_mid_l
 
    if ( ubound(l_late_l,1) == inumberoflanduses )          &
      growth_stage_length_in_days( l_doy_late, : ) = l_late_l
 
    if ( ubound(l_fallow_l,1) == inumberoflanduses )        &
      growth_stage_length_in_days( l_doy_fallow, : ) = l_fallow_l
 
    call logs%write(" ## Crop Kcb Curve Summary ##", ilinesafter=1)
    call logs%write(" _only meaningful for landuses where the Kcb curve is defined " &
      //"in terms of days _", ilinesafter=1)
    call logs%write("Landuse Code | Planting Date | End of 'ini' | End of 'dev' " &
      //"| End of 'mid' | End of 'late' | End of 'fallow' ")
    call logs%write("-------------|---------------|--------------|--------------" &
      //"|--------------|---------------|-----------------")
 
    if ( sl_planting_date%count == inumberoflanduses .and. sl_planting_date%count > 0 ) then
 
      do iindex=1, sl_planting_date%count
 
        plantingdate_str = sl_planting_date%get( iindex )
 
        ! if there is no planting date entry, assume user is specifying
        ! a planting GDD for this landuse code
        if ( len_trim(plantingdate_str) == 0 ) cycle
 
        ! if ( PlantingDate_str .contains. "/" ) then
 
        !   ! append current year to the end of the user-entered planting date in mm/dd
        !   sMMDDYYYY = trim(PlantingDate_str)//"/"//asCharacter( SIM_DT%start%iYear )
        !   call GROWTH_STAGE_DATE( PLANTING_DATE, iIndex)%parsedate( sMMDDYYYY, __FILE__, __LINE__ )
 
        !   GROWTH_STAGE_DATE( PLANTING_DATE, iIndex) = GROWTH_STAGE_DATE( PLANTING_DATE, iIndex)
 
        ! else
        !   ! assume the value is a day-of-year value
        !   dtPlantingDate = SIM_DT%start + asFloat( PlantingDate_str)
        !   call dtPlantingDate%calcJulianDay()
        !   GROWTH_STAGE_DATE( PLANTING_DATE, iIndex) = dtPlantingDate
        ! endif
 
        growth_stage_date( planting_date, iindex) =                                     &
            crop_coefficients_calculate_planting_date(splantingdate=plantingdate_str,   &
                                                      iyear=sim_dt%start%iYear)
 
        ! march forward through time calculating the various dates on the Kcb curve
        ! GROWTH_STAGE_DATE( ENDDATE_INI, iIndex ) = GROWTH_STAGE_DATE( PLANTING_DATE, iIndex ) + L_ini_l( iIndex )
        ! GROWTH_STAGE_DATE( ENDDATE_DEV, iIndex ) = GROWTH_STAGE_DATE( ENDDATE_INI, iIndex ) + L_dev_l( iIndex )
        ! GROWTH_STAGE_DATE( ENDDATE_MID, iIndex ) = GROWTH_STAGE_DATE( ENDDATE_DEV, iIndex ) + L_mid_l( iIndex )
        ! GROWTH_STAGE_DATE( ENDDATE_LATE, iIndex ) = GROWTH_STAGE_DATE( ENDDATE_MID, iIndex ) + L_late_l( iIndex )
        ! GROWTH_STAGE_DATE( ENDDATE_FALLOW, iIndex ) = GROWTH_STAGE_DATE( ENDDATE_LATE, iIndex ) + L_fallow_l( iIndex )
 
        ! if any of the L_* length values is missing, a value of zero will be used, resulting in a wierd looking Kcb curve
        growth_stage_date( enddate_ini, iindex ) = growth_stage_date( planting_date, iindex )                        &
                                                   + (growth_stage_length_in_days( l_doy_ini, iindex ) )
        growth_stage_date( enddate_dev, iindex ) = growth_stage_date( enddate_ini, iindex )                          &
                                                   + (growth_stage_length_in_days( l_doy_dev, iindex ) )
        growth_stage_date( enddate_mid, iindex ) = growth_stage_date( enddate_dev, iindex )                          &
                                                   + (growth_stage_length_in_days( l_doy_mid, iindex ) )
        growth_stage_date( enddate_late, iindex ) = growth_stage_date( enddate_mid, iindex )                         &
                                                   + (growth_stage_length_in_days( l_doy_late, iindex ) )
        growth_stage_date( enddate_fallow, iindex ) = growth_stage_date( enddate_late, iindex )                      &
                                                     + (growth_stage_length_in_days( l_doy_fallow, iindex ) )
 
        call logs%write( "| "//ascharacter( landuse_code( iindex ))//" | "                                    &
           //trim( growth_stage_date( planting_date, iindex )%prettydate() )                                  &
             //" (doy:"//ascharacter( growth_stage_date( planting_date, iindex )%getDayOfYear() )//") | "     &
           //trim( growth_stage_date( enddate_ini, iindex )%prettydate() )//" | "                             &
             //" (doy:"//ascharacter( growth_stage_date( enddate_ini, iindex )%getDayOfYear() )//") | "       &
           //trim( growth_stage_date( enddate_dev, iindex )%prettydate() )//" | "                             &
             //" (doy:"//ascharacter( growth_stage_date( enddate_dev, iindex )%getDayOfYear() )//") | "       &
           //trim( growth_stage_date( enddate_mid, iindex )%prettydate() )//" | "                             &
             //" (doy:"//ascharacter( growth_stage_date( enddate_mid, iindex )%getDayOfYear() )//") | "       &
           //trim( growth_stage_date( enddate_late, iindex )%prettydate() )//" | "                            &
             //" (doy:"//ascharacter( growth_stage_date( enddate_late, iindex )%getDayOfYear() )//") | "      &
           //trim( growth_stage_date( enddate_fallow, iindex )%prettydate() )                                 &
             //" (doy:"//ascharacter( growth_stage_date( enddate_fallow, iindex )%getDayOfYear() )//") | ")
      enddo
 
    endif
 
    if (ubound(gdd_plant_l,1) == inumberoflanduses)  growth_stage_gdd( gdd_plant,  : ) = gdd_plant_l
    if (ubound(gdd_ini_l,1) == inumberoflanduses)    growth_stage_gdd( gdd_ini,  : ) = gdd_ini_l
    if (ubound(gdd_dev_l,1) == inumberoflanduses)    growth_stage_gdd( gdd_dev,  : ) = gdd_dev_l
    if (ubound(gdd_mid_l,1) == inumberoflanduses)    growth_stage_gdd( gdd_mid,  : ) = gdd_mid_l
    if (ubound(gdd_late_l,1) == inumberoflanduses)   growth_stage_gdd( gdd_late, : ) = gdd_late_l
 
    if (ubound(kcb_ini_l,1) == inumberoflanduses)  kcb_l( kcb_ini, :) = kcb_ini_l
    if (ubound(kcb_mid_l,1) == inumberoflanduses)  kcb_l( kcb_mid, :) = kcb_mid_l
    if (ubound(kcb_end_l,1) == inumberoflanduses)  kcb_l( kcb_end, :) = kcb_end_l
    if (ubound(kcb_min_l,1) == inumberoflanduses)  kcb_l( kcb_min, :) = kcb_min_l
 
    if (ubound(kcb_jan,1) == inumberoflanduses)   kcb_l( jan, :) = kcb_jan
    if (ubound(kcb_feb,1) == inumberoflanduses)   kcb_l( feb, :) = kcb_feb
    if (ubound(kcb_mar,1) == inumberoflanduses)   kcb_l( mar, :) = kcb_mar
    if (ubound(kcb_apr,1) == inumberoflanduses)   kcb_l( apr, :) = kcb_apr
    if (ubound(kcb_may,1) == inumberoflanduses)   kcb_l( may, :) = kcb_may
    if (ubound(kcb_jun,1) == inumberoflanduses)   kcb_l( jun, :) = kcb_jun
    if (ubound(kcb_jul,1) == inumberoflanduses)   kcb_l( jul, :) = kcb_jul
    if (ubound(kcb_aug,1) == inumberoflanduses)   kcb_l( aug, :) = kcb_aug
    if (ubound(kcb_sep,1) == inumberoflanduses)   kcb_l( sep, :) = kcb_sep
    if (ubound(kcb_oct,1) == inumberoflanduses)   kcb_l( oct, :) = kcb_oct
    if (ubound(kcb_nov,1) == inumberoflanduses)   kcb_l( nov, :) = kcb_nov
    if (ubound(kcb_dec,1) == inumberoflanduses)   kcb_l( dec, :) = kcb_dec
 
    ! go through the table values and try to figure out how Kcb curves should be constructed:
    ! Monthly Kcb, GDD-based, or DOY-based
    do iindex = lbound( kcb_method, 1), ubound( kcb_method, 1)
 
      if ( all( kcb_l( jan:dec, iindex ) > 0.0_c_float ) ) then
        kcb_method( iindex ) = kcb_method_monthly_values
        kcb_l( kcb_min, iindex ) = minval( kcb_l(jan:dec, iindex) )
        kcb_l( kcb_mid, iindex) = minval( kcb_l(jan:dec, iindex) )
 
      elseif ( all( growth_stage_gdd( :, iindex ) >= 0.0_c_float )              &
         .and. all( kcb_l( kcb_ini:kcb_min, iindex ) > 0.0_c_float ) ) then
        kcb_method( iindex ) = kcb_method_gdd
 
      elseif ( all( growth_stage_length_in_days( planting_date:, iindex ) >= 0.0_c_float )              &
         .and. all( kcb_l( kcb_ini:kcb_min, iindex ) > 0.0_c_float ) ) then
        kcb_method( iindex ) = kcb_method_fao56
      endif
 
      if ( kcb_method( iindex ) < 0 ) then
        call warn("There are missing day-of-year (L_ini, L_dev, L_mid, L_late, L_fallow), " &
          //"growing degree-day ~(GDD_plant, GDD_ini, GDD_dev, GDD_mid, GDD_late)," &
          //" or monthly crop ~coefficients (Kcb_jan...Kcb_dec) for" &
          //" landuse "//ascharacter( landuse_code( iindex ) ), lfatal=true )
      endif
 
    enddo
 
!    do iIndex = lbound( fSoilStorage, 1 ), ubound( fSoilStorage,1 )
 
!      fKcb_initial = update_crop_coefficient_date_as_threshold( iLanduseIndex( iIndex ) )
 
      ! call calc_effective_root_depth( fRz_i=fRz_initial, iLanduseIndex=iLanduseIndex( iIndex ),    &
      !                                 fZr_max=fMax_Rooting_Depths( iLanduseIndex( iIndex ),        &
      !                                 iSoilGroup( iIndex ) ),                                      &
      !                                 Kcb=fKcb_initial )
      !
      ! fSoilStorage( iIndex ) = INITIAL_PERCENT_SOIL_MOISTURE( iIndex ) / 100.0_c_float             &
      !                          * fRz_initial * fAvailable_Water_Content( iIndex )
 
!    enddo
 
  !> @TODO Add more logic here to perform checks on the validity of this data.
 
  !> @TODO Need to handle missing values. WHat do we do if an entire column of values
  !!       is missing?
 
  end subroutine crop_coefficients_fao56_initialize
 
!------------------------------------------------------------------------------
 
 !> Update the current basal crop coefficient (Kcb) for
 !! a SINGLE irrigation table entry
 !!
 !! @param[inout] pIRRIGATION pointer to a single line of information in the irrigation file.
 !! @param[in] iThreshold either the current day of year or the number of growing degree days.
 !! @retval rKcb Basal crop coefficient given the irrigation table entries and the
 !!         current threshold values.
 
 impure elemental function update_crop_coefficient_date_as_threshold( iLanduseIndex )           &
                                                                          result(kcb)
 
  integer (c_int), intent(in)   :: ilanduseindex
  real (c_float)                :: kcb
 
  ! [ LOCALS ]
  real (c_double) :: ffrac
 
  if ( kcb_method( ilanduseindex ) == kcb_method_monthly_values ) then
 
    kcb = kcb_l( sim_dt%curr%iMonth, ilanduseindex )
 
  else
 
    ! define shorthand variable names for remainder of function
    associate( date_ini => growth_stage_date( enddate_ini, ilanduseindex ),           &
                date_dev => growth_stage_date( enddate_dev, ilanduseindex ),           &
                date_mid => growth_stage_date( enddate_mid, ilanduseindex ),           &
                date_late => growth_stage_date( enddate_late, ilanduseindex ),         &
                date_fallow => growth_stage_date( enddate_fallow, ilanduseindex ),     &
                kcb_ini => kcb_l(kcb_ini, ilanduseindex),                              &
                kcb_mid => kcb_l(kcb_mid, ilanduseindex),                              &
                kcb_min => kcb_l(kcb_min, ilanduseindex),                              &
                plantingdate => growth_stage_date( planting_date, ilanduseindex),      &
                kcb_end => kcb_l(kcb_end, ilanduseindex),                              &
                current_date => sim_dt%curr )
 
      ! now calculate Kcb for the given landuse
 
      if( current_date >= date_late ) then
 
        kcb = kcb_min
 
      elseif ( current_date >= date_mid ) then
 
        ffrac = ( current_date - date_mid ) / ( date_late - date_mid )
 
        kcb =  real(kcb_mid * (1.0_c_double - ffrac) + kcb_end * ffrac, c_float)
 
      elseif ( current_date >= date_dev ) then
 
        kcb = kcb_mid
 
      elseif ( current_date >= date_ini ) then
 
        ffrac = ( current_date - date_ini ) / ( date_dev - date_ini )
 
        kcb = real(kcb_ini * (1.0_c_double - ffrac) + kcb_mid * ffrac, c_float)
 
      elseif ( current_date >= plantingdate ) then
 
        kcb = kcb_ini
 
      else
 
        kcb = kcb_min
 
      endif
 
    end associate
 
  end if
 
end function update_crop_coefficient_date_as_threshold
 
!------------------------------------------------------------------------------
 
pure elemental function crop_coefficients_fao56_calculate_kcb_max(wind_speed_meters_per_sec,   &
                                          relative_humidity_min_pct,   &
                                          Kcb,                         & 
                                          plant_height_meters)                       result(kcb_max)
 
  real (c_float), intent(in) :: wind_speed_meters_per_sec
  real (c_float), intent(in) :: relative_humidity_min_pct
  real (c_float), intent(in) :: kcb
  real (c_float), intent(in) :: plant_height_meters
 
  real (c_float)  :: kcb_max
  real (c_double) :: u2
  real (c_double) :: rhmin
  real (c_double) :: plant_height
 
  ! Limits are as suggested on page 123 of FAO-56 with respect to
  ! modifying mid-season KCB_mid values 
  rhmin = clip( relative_humidity_min_pct, minval=20., maxval=80. )
  u2 = clip(wind_speed_meters_per_sec, minval=1., maxval=6.)
  plant_height = clip(plant_height_meters, minval=1., maxval=10.)
 
  ! equation 72, FAO-56, p 199
  kcb_max = real(max(  1.2_c_double + ( (0.04_c_double * (u2 - 2._c_double)               &
                                  - 0.004_c_double * (rhmin - 45._c_double) ) )      &
                                  * (plant_height_meters/3._c_double)**0.3_c_double, &
                  kcb + 0.05_c_double ), c_float)
 
end function crop_coefficients_fao56_calculate_kcb_max
 
!------------------------------------------------------------------------------
 
 !> Update the current basal crop coefficient (Kcb), with GDD as the threhold
 !!
 !! @param[in] fGDD current growing degree day value associated with the cell.
 !! @retval fKcb Basal crop coefficient given the irrigation table entries and the
 !!         current threshold values.
 
 impure elemental function update_crop_coefficient_gdd_as_threshold( iLanduseIndex, fGDD )   &
                                                                         result(fkcb)
 
  integer (c_int), intent(in)   :: ilanduseindex
  real (c_float), intent(in)    :: fgdd
  real (c_float)                :: fkcb
 
  ! [ LOCALS ]
  real (c_double) :: ffrac
 
  ! define shorthand variable names for remainder of function
  associate( gdd_ini_l => growth_stage_gdd( gdd_ini, ilanduseindex ),          &
              gdd_dev_l => growth_stage_gdd( gdd_dev, ilanduseindex ),          &
              gdd_mid_l => growth_stage_gdd( gdd_mid, ilanduseindex ),          &
              gdd_late_l => growth_stage_gdd( gdd_late, ilanduseindex ),        &
              kcb_ini => kcb_l(kcb_ini, ilanduseindex),                         &
              kcb_mid => kcb_l(kcb_mid, ilanduseindex),                         &
              kcb_min => kcb_l(kcb_min, ilanduseindex),                         &
              plantingdoy => growth_stage_date( planting_date,                  &
                                                ilanduseindex)%getDayOfYear(),  &
              gdd_plant_l => growth_stage_gdd( gdd_plant, ilanduseindex),       &
              kcb_end => kcb_l(kcb_end, ilanduseindex),                         &
              current_doy => sim_dt%curr%getDayOfYear() )
 
    ! now calculate Kcb for the given landuse
    if( fgdd > gdd_late_l ) then
 
      fkcb = kcb_min
 
    elseif ( fgdd > gdd_mid_l ) then
 
      ffrac = ( fgdd - gdd_mid_l ) / ( gdd_late_l - gdd_mid_l )
 
      fkcb =  real(kcb_mid * (1.0_c_double - ffrac) + kcb_end * ffrac, c_float)
 
    elseif ( fgdd > gdd_dev_l ) then
 
      fkcb = kcb_mid
 
    elseif ( fgdd > gdd_ini_l ) then
 
      ffrac = ( fgdd - gdd_ini_l ) / ( gdd_dev_l - gdd_ini_l )
 
      fkcb = real(kcb_ini * (1_c_double - ffrac) + kcb_mid * ffrac, c_float)
 
    elseif ( (plantingdoy > 0) .and. (current_doy >= plantingdoy) ) then
 
      ! if there is a valid value for the Planting Date, use it
      fkcb = kcb_ini
 
    elseif ( (gdd_plant_l >= 0.) .and. (fgdd >= gdd_plant_l) ) then
 
        fkcb = kcb_ini
 
    else
 
      fkcb = kcb_min
 
    endif
 
  end associate
 
end function update_crop_coefficient_gdd_as_threshold
 
!------------------------------------------------------------------------------
 
  subroutine crop_coefficients_fao56_update_growth_stage_dates( )
 
    ! [ LOCALS ]
    integer (c_int) :: iindex
 
    do iindex=lbound(growth_stage_date,2), ubound(growth_stage_date,2)
 
  !     print *, SIM_DT%curr%prettydate(), " | ", GROWTH_STAGE_DATE( ENDDATE_LATE, iIndex )%prettydate()
 
      ! if we have not yet reached the enddate associated with the fallow period, skip
      if ( sim_dt%curr < growth_stage_date( enddate_fallow, iindex ) ) cycle
 
      if ( kcb_method( iindex ) /= kcb_method_fao56 ) cycle
 
      ! current date is beyond the enddate associated with fallow period;
      ! update Kcb curve and dates
 
!      print *, "a) ",SIM_DT%curr%prettydate(), " | ", GROWTH_STAGE_DATE( ENDDATE_LATE, iIndex )%prettydate()
 
      ! it's possible that the planting date might be later in the current calendar year
      !call GROWTH_STAGE_DATE( PLANTING_DATE, iIndex )%setYear( SIM_DT%curr%iYear )
 
!      print *, "b) ",SIM_DT%curr%prettydate(), " | ", GROWTH_STAGE_DATE( PLANTING_DATE, iIndex )%prettydate()
 
      ! if we are already past the planting doy, planting date must be
      ! next calendar year
      if ( sim_dt%iDOY > growth_stage_date( planting_date, iindex )%getDayOfYear() )  then
 
        growth_stage_date( planting_date, iindex ) =                                                &
            crop_coefficients_calculate_planting_date(splantingdate = sl_planting_date%get(iindex), &
                                                      iyear=sim_dt%curr%iYear + 1)
 
      else
 
        ! allow for planting date to be later in the current year
        growth_stage_date( planting_date, iindex ) =                                                &
            crop_coefficients_calculate_planting_date(splantingdate = sl_planting_date%get(iindex), &
                                                      iyear=sim_dt%curr%iYear )
        
      endif 
 
      ! now calculate dates associated with the rest of the Kcb curve
      growth_stage_date( enddate_ini, iindex ) = growth_stage_date( planting_date, iindex ) &
                                                + (growth_stage_length_in_days( l_doy_ini, iindex ) )
      growth_stage_date( enddate_dev, iindex ) = growth_stage_date( enddate_ini, iindex )&
                                                + (growth_stage_length_in_days( l_doy_dev, iindex ) )
      growth_stage_date( enddate_mid, iindex ) = growth_stage_date( enddate_dev, iindex )&
                                                + (growth_stage_length_in_days( l_doy_mid, iindex ) )
      growth_stage_date( enddate_late, iindex ) = growth_stage_date( enddate_mid, iindex )&
                                                + (growth_stage_length_in_days( l_doy_late, iindex ) )
      growth_stage_date( enddate_fallow, iindex ) = growth_stage_date( enddate_late, iindex )&
                                                + (growth_stage_length_in_days( l_doy_fallow, iindex ) )
 
      call logs%write("## Updated Crop Kcb Curve Summary ##", ilinesafter=1, lecho=false )
      call logs%write("Landuse Code | Planting Date | End of 'ini' | End of 'dev' " &
        //"| End of 'mid' | End of 'late' | End of 'fallow' ", lecho=false )
      call logs%write("-------------|---------------|--------------|--------------" &
        //"|--------------|---------------|-----------------", lecho=false )
 
      ! call LOGS%write( asCharacter( LANDUSE_CODE( iIndex ))//" | "                &
      !    //trim( GROWTH_STAGE_DATE( PLANTING_DATE, iIndex )%prettydate() )//" | " &
      !    //trim( GROWTH_STAGE_DATE( ENDDATE_INI, iIndex )%prettydate() )//" | "   &
      !    //trim( GROWTH_STAGE_DATE( ENDDATE_DEV, iIndex )%prettydate() )//" | "   &
      !    //trim( GROWTH_STAGE_DATE( ENDDATE_MID, iIndex )%prettydate() )//" | "   &
      !    //trim( GROWTH_STAGE_DATE( ENDDATE_LATE, iIndex )%prettydate() )//" | "  &
      !    //trim( GROWTH_STAGE_DATE( ENDDATE_FALLOW, iIndex )%prettydate() ),      &
      !    lEcho=FALSE, iLogLevel=LOG_ALL )
 
 
        call logs%write( "| "//ascharacter( landuse_code( iindex ))//" | "                                    &
           //trim( growth_stage_date( planting_date, iindex )%prettydate() )                                  &
             //" (doy:"//ascharacter( growth_stage_date( planting_date, iindex )%getDayOfYear() )//") | "     &
           //trim( growth_stage_date( enddate_ini, iindex )%prettydate() )//" | "                             &
             //" (doy:"//ascharacter( growth_stage_date( enddate_ini, iindex )%getDayOfYear() )//") | "       &
           //trim( growth_stage_date( enddate_dev, iindex )%prettydate() )//" | "                             &
             //" (doy:"//ascharacter( growth_stage_date( enddate_dev, iindex )%getDayOfYear() )//") | "       &
           //trim( growth_stage_date( enddate_mid, iindex )%prettydate() )//" | "                             &
             //" (doy:"//ascharacter( growth_stage_date( enddate_mid, iindex )%getDayOfYear() )//") | "       &
           //trim( growth_stage_date( enddate_late, iindex )%prettydate() )//" | "                            &
             //" (doy:"//ascharacter( growth_stage_date( enddate_late, iindex )%getDayOfYear() )//") | "      &
           //trim( growth_stage_date( enddate_fallow, iindex )%prettydate() )                                 &
             //" (doy:"//ascharacter( growth_stage_date( enddate_fallow, iindex )%getDayOfYear() )//") | ",   &
           lecho=false, iloglevel=log_all)
 
    enddo
 
 
  end subroutine crop_coefficients_fao56_update_growth_stage_dates
 
!--------------------------------------------------------------------------------------------------
 
  impure elemental subroutine crop_coefficients_fao56_calculate( Kcb, landuse_index, GDD )
 
    real (c_float), intent(inout)          :: kcb
    integer (c_int), intent(in)            :: landuse_index
    real (c_float), intent(in), optional   :: gdd
 
 
    if ( kcb_method( landuse_index )  == kcb_method_fao56  &
      .or. kcb_method( landuse_index ) == kcb_method_monthly_values ) then
 
      kcb = update_crop_coefficient_date_as_threshold( landuse_index )
 
    else
 
      kcb = update_crop_coefficient_gdd_as_threshold( landuse_index, gdd )
 
    endif
 
  end subroutine crop_coefficients_fao56_calculate
 
!--------------------------------------------------------------------------------------------------
 
  impure elemental subroutine crop_coefficients_fao56_update_growing_season(   &
                                           landuse_index,                      &
                                           Kcb,                                &
                                           it_is_growing_season)
 
    real (c_float), intent(in)             :: kcb
    integer (c_int), intent(in)            :: landuse_index
    logical (c_bool), intent(out)          :: it_is_growing_season
 
    if ( kcb > kcb_l( kcb_min, landuse_index) ) then
      it_is_growing_season = true
    else
      it_is_growing_season = false
    endif
 
  end subroutine crop_coefficients_fao56_update_growing_season
 
end module crop_coefficients__fao56