Merge pull request ESCOMP#6 from billsacks/agsys

Rework phases and params derived types

src/agsys/ctsm_interface/AgSysClimateInterface.F90

@@ -9,7 +9,7 @@ module AgSysClimateInterface
   use shr_kind_mod     , only : r8 => shr_kind_r8
   use shr_infnan_mod   , only : nan => shr_infnan_nan, assignment(=)
   use decompMod        , only : bounds_type
-  use AgSysPhases      , only : agsys_phases_type
+  use AgSysRuntimeConstants, only : agsys_max_phases
   !
   implicit none
   private
@@ -70,34 +70,32 @@ end subroutine AgSysClimateTimeStep
   ! ========================================================================
 
   !-----------------------------------------------------------------------
-  subroutine Init(this, bounds, agsys_phases_inst)
+  subroutine Init(this, bounds)
     !
     ! !DESCRIPTION:
     ! Initialize this agsys_climate_type instance
     !
     ! !ARGUMENTS:
     class(agsys_climate_type), intent(inout) :: this
     type(bounds_type), intent(in) :: bounds
-    type(agsys_phases_type), intent(in) :: agsys_phases_inst
     !
     ! !LOCAL VARIABLES:
 
     character(len=*), parameter :: subname = 'Init'
     !-----------------------------------------------------------------------
 
-    call this%InitAllocate(bounds, agsys_phases_inst)
+    call this%InitAllocate(bounds)
   end subroutine Init
 
   !-----------------------------------------------------------------------
-  subroutine InitAllocate(this, bounds, agsys_phases_inst)
+  subroutine InitAllocate(this, bounds)
     !
     ! !DESCRIPTION:
     ! Allocate components of this agsys_climate_type instance
     !
     ! !ARGUMENTS:
     class(agsys_climate_type), intent(inout) :: this
     type(bounds_type), intent(in) :: bounds
-    type(agsys_phases_type), intent(in) :: agsys_phases_inst
     !
     ! !LOCAL VARIABLES:
 
@@ -112,7 +110,7 @@ subroutine InitAllocate(this, bounds, agsys_phases_inst)
     allocate(this%accumulated_thermal_time_patch(begp:endp))         ; this%accumulated_thermal_time_patch(:) = nan
     allocate(this%accumulated_emerged_thermal_time_patch(begp:endp)) ; this%accumulated_emerged_thermal_time_patch(:) = nan
 
-    allocate(this%accumulated_thermal_time_phases_patch(1:agsys_phases_inst%max_phases, begp:endp))
+    allocate(this%accumulated_thermal_time_phases_patch(1:agsys_max_phases, begp:endp))
     this%accumulated_thermal_time_phases_patch(:,:) = nan
 
     end associate

src/agsys/ctsm_interface/AgSysGeneralType.F90

@@ -26,12 +26,10 @@ module AgSysGeneralType
      ! that these may differ from the constants in pftconMod
      integer, pointer, public :: crop_type_patch(:)
 
-     ! Cultivar type. A given value implies a given crop type; for example, cultivar
-     ! values 1-3 may always be maize cultivars, and cultivar values 4-5 always soybean
-     ! cultivars, etc. - so you can index cultivar-specific parameters just with
-     ! cultivar_patch, without needing to also reference crop_type_patch. Each crop type
-     ! has at least one cultivar, and may have many. This is currently constant in time,
-     ! but eventually may be dynamic.
+     ! Cultivar type. For a given crop type, the cultivar type numbering starts at 1. So
+     ! maize may have cultivars 1-3, soybean 1-2, wheat 1-4, etc. Each crop type has at
+     ! least one cultivar, and may have many. This is currently constant in time, but
+     ! eventually may be dynamic.
      integer, pointer, public :: cultivar_patch(:)
 
      ! ------------------------------------------------------------------------

src/agsys/ctsm_interface/AgSysInterface.F90

@@ -10,10 +10,12 @@ module AgSys
 #include "shr_assert.h"
   use clm_time_manager, only : is_beg_curr_day
   use decompMod, only : bounds_type
+  use AgSysConstants, only : crop_type_maxval
   use AgSysGeneralType, only : agsys_general_type
-  use AgSysParams, only : agsys_params_type
+  use AgSysParams, only : agsys_crop_params_type, agsys_crop_cultivar_params_type
   use AgSysPhases, only : agsys_phases_type
   use AgSysParamReader, only : ReadParams, ReadPhases
+  use AgSysRuntimeConstants, only : InitRuntimeConstants
   use AgSysClimateInterface, only : agsys_climate_type
   use AgSysPhenologyInterface, only : agsys_phenology_type
   !
@@ -24,16 +26,37 @@ module AgSys
 
   type, public :: agsys_type
      private
+     ! Parameters that vary by crop type
+     type(agsys_crop_params_type) :: crop_params(crop_type_maxval)
+
+     ! Parameters that vary by cultivar; these are first indexed by crop type, then
+     ! indexed by the specific cultivar for that crop. For example:
+     !
+     ! do p = ...  <- loop over patches (typically a filter loop)
+     !    crop_type = agsys_general_inst%crop_type_patch(p)
+     !    cultivar  = agsys_general_inst%cultivar_patch(p)
+     !    call SomeAgsysRoutine( &
+     !         crop_params     = agsys_inst%crop_params(crop_type), &
+     !         cultivar_params = agsys_inst%crop_cultivar_params(crop_type)%cultivar_params(cultivar), &
+     !         ...)
+     ! end do
+     !
+     ! (Note that this way, the science code - such as SomeAgsysRoutine - doesn't need to
+     ! index the params by crop type or cultivar: a given call for a single patch only
+     ! has access to the parameters for the appropriate crop type and cultivar for that
+     ! patch.)
+     type(agsys_crop_cultivar_params_type) :: crop_cultivar_params(crop_type_maxval)
+
+     ! Information about the phases for each crop type
+     type(agsys_phases_type) :: crop_phases(crop_type_maxval)
+
      type(agsys_general_type)   :: agsys_general_inst
-     type(agsys_params_type)    :: agsys_params_inst
-     type(agsys_phases_type)    :: agsys_phases_inst
      type(agsys_climate_type)   :: agsys_climate_inst
      type(agsys_phenology_type) :: agsys_phenology_inst
 
    contains
      procedure, public :: AgSysDriver
      procedure, public :: Init
-
   end type agsys_type
 
   character(len=*), parameter, private :: sourcefile = &
@@ -89,10 +112,11 @@ subroutine Init(this, bounds)
     character(len=*), parameter :: subname = 'Init'
     !-----------------------------------------------------------------------
 
+    call ReadParams(this%crop_params, this%crop_cultivar_params)
+    call ReadPhases(this%crop_phases)
+    call InitRuntimeConstants(this%crop_phases)
     call this%agsys_general_inst%Init(bounds)
-    call ReadParams(this%agsys_params_inst)
-    call ReadPhases(this%agsys_phases_inst)
-    call this%agsys_climate_inst%Init(bounds, this%agsys_phases_inst)
+    call this%agsys_climate_inst%Init(bounds)
     call this%agsys_phenology_inst%Init(bounds)
 
   end subroutine Init

src/agsys/ctsm_interface/AgSysParamReader.F90

@@ -8,8 +8,9 @@ module AgSysParamReader
 #include "shr_assert.h"
   use shr_kind_mod     , only : r8 => shr_kind_r8
   use shr_infnan_mod   , only : nan => shr_infnan_nan, assignment(=)
-  use AgSysParams      , only : agsys_params_type
+  use AgSysParams      , only : agsys_crop_params_type, agsys_crop_cultivar_params_type
   use AgSysPhases      , only : agsys_phases_type
+  use AgSysConstants   , only : crop_type_maxval
   !
   implicit none
   private
@@ -24,35 +25,41 @@ module AgSysParamReader
 contains
 
   !-----------------------------------------------------------------------
-  subroutine ReadParams(params)
+  subroutine ReadParams(crop_params, crop_cultivar_params)
     !
     ! !DESCRIPTION:
     ! Read parameters
     !
     ! !ARGUMENTS:
-    type(agsys_params_type), intent(inout) :: params
+    type(agsys_crop_params_type), intent(inout) :: crop_params(:)
+    type(agsys_crop_cultivar_params_type), intent(inout) :: crop_cultivar_params(:)
     !
     ! !LOCAL VARIABLES:
 
     character(len=*), parameter :: subname = 'ReadParams'
     !-----------------------------------------------------------------------
 
+    SHR_ASSERT_FL((size(crop_params) == crop_type_maxval), sourcefile, __LINE__)
+    SHR_ASSERT_FL((size(crop_cultivar_params) == crop_type_maxval), sourcefile, __LINE__)
+
   end subroutine ReadParams
 
   !-----------------------------------------------------------------------
-  subroutine ReadPhases(phases)
+  subroutine ReadPhases(crop_phases)
     !
     ! !DESCRIPTION:
     ! Read phase descriptions for each crop
     !
     ! !ARGUMENTS:
-    type(agsys_phases_type), intent(inout) :: phases
+    type(agsys_phases_type), intent(inout) :: crop_phases(:)
     !
     ! !LOCAL VARIABLES:
 
     character(len=*), parameter :: subname = 'ReadPhases'
     !-----------------------------------------------------------------------
 
+    SHR_ASSERT_FL((size(crop_phases) == crop_type_maxval), sourcefile, __LINE__)
+
   end subroutine ReadPhases
 
 end module AgSysParamReader

src/agsys/ctsm_interface/AgSysRuntimeConstants.F90

@@ -0,0 +1,54 @@
+module AgSysRuntimeConstants
+
+  !-----------------------------------------------------------------------
+  ! !DESCRIPTION:
+  ! Runtime constants needed throughout AgSys. These are NOT crop-specific parameters
+  ! (which are stored elsewhere); rather, these are more general constants.
+  !
+  ! Constants should only be added to this module sparingly: this should be reserved for
+  ! constants whose use is widespread, inside and/or outside of AgSys (e.g., to define
+  ! dimension sizes on restart and/or history files).
+  !
+  ! !USES:
+  use AgSysPhases, only : agsys_phases_type
+  !
+  implicit none
+  private
+  save
+
+  ! !PUBLIC DATA:
+  integer, public, protected :: agsys_max_phases ! maximum number of phases used by any crop
+
+  ! !PUBLIC ROUTINES:
+  public :: InitRuntimeConstants ! Initialize runtime constants in this module
+
+  character(len=*), parameter, private :: sourcefile = &
+       __FILE__
+
+contains
+
+  !-----------------------------------------------------------------------
+  subroutine InitRuntimeConstants(crop_phases)
+    !
+    ! !DESCRIPTION:
+    ! Initialize runtime constants in this module
+    !
+    ! !ARGUMENTS:
+    type(agsys_phases_type), intent(in) :: crop_phases(:) ! phases for each crop
+    !
+    ! !LOCAL VARIABLES:
+    integer :: crop
+
+    character(len=*), parameter :: subname = 'InitRuntimeConstants'
+    !-----------------------------------------------------------------------
+
+    agsys_max_phases = 0
+    do crop = 1, ubound(crop_phases,1)
+       if (crop_phases(crop)%num_phases > agsys_max_phases) then
+          agsys_max_phases = crop_phases(crop)%num_phases
+       end if
+    end do
+
+  end subroutine InitRuntimeConstants
+
+end module AgSysRuntimeConstants

src/agsys/science/AgSysConstants.F90

@@ -16,7 +16,7 @@ module AgSysConstants
   integer, parameter, public :: crop_type_maize = 1
   integer, parameter, public :: crop_type_soybean = 2
   integer, parameter, public :: crop_type_wheat = 3
-  integer, parameter, public :: crop_type_soyghum = 4
+  integer, parameter, public :: crop_type_sorghum = 4
   integer, parameter, public :: crop_type_maxval = 4
 
 end module AgSysConstants

src/agsys/science/AgSysParams.F90

@@ -14,21 +14,41 @@ module AgSysParams
   implicit none
   private
 
+  ! !PUBLIC TYPES:
+
   type, public :: response_curve_type
-     integer :: num_pts
-     real(r8), allocatable :: x(:)
-     real(r8), allocatable :: y(:)
+     private
+
+     ! Public data members
+     integer, public :: num_pts
+     real(r8), allocatable, public :: x(:)
+     real(r8), allocatable, public :: y(:)
   end type response_curve_type
 
-  ! !PUBLIC TYPES:
-  type, public :: agsys_params_type
+  ! Parameters that vary by crop
+  type, public :: agsys_crop_params_type
      private
 
      ! Public data members
-     real(r8), allocatable, public :: shoot_lag_cultivar(:)
-     real(r8), allocatable, public :: shoot_rate_cultivar(:)
+     real(r8), public :: shoot_lag
+     real(r8), public :: shoot_rate
+  end type agsys_crop_params_type
 
-     type(response_curve_type), allocatable, public :: target_tt_from_photoperiod_end_of_juvenile_cultivar(:)
-  end type agsys_params_type
+  ! Parameters that vary by cultivar
+  type, public :: agsys_cultivar_params_type
+     private
+
+     ! Public data members
+     type(response_curve_type), public :: target_tt_from_photoperiod_end_of_juvenile
+  end type agsys_cultivar_params_type
+
+  ! Each crop has its own vector of cultivar-specific parameters. There is one instance
+  ! of this derived type for each crop.
+  type, public :: agsys_crop_cultivar_params_type
+     private
+
+     ! Public data members
+     type(agsys_cultivar_params_type), allocatable :: cultivar_params(:)
+  end type agsys_crop_cultivar_params_type
 
 end module AgSysParams

src/agsys/science/AgSysPhases.F90

@@ -25,27 +25,70 @@ module AgSysPhases
   integer, parameter, public :: phase_type_end = 7
   integer, parameter, public :: phase_type_maxval = 7
 
+  integer, parameter, public :: composite_phase_type_vernalization = 1
+  integer, parameter, public :: composite_phase_type_emerge_to_end_of_juvenile = 2
+  integer, parameter, public :: composite_phase_type_maxval = 2
+
   integer, parameter, public :: max_phase_name_len = 64  ! maximum number of characters in names of phases / stages
 
   ! !PUBLIC TYPES:
+
+  ! Information for a single composite phase for a given crop. A composite phase is
+  ! something like vernalization, which is defined as a period spanning multiple real
+  ! phases. The real phases included in a composite phase may be contiguous or
+  ! non-contiguous.
+  type, public :: composite_phase_type
+     private
+
+     ! Public data members
+     character(len=max_phase_name_len), public              :: name                ! name of this composite phase
+     integer, public                                        :: num_child_phases    ! number of child phases (arrays are dimensioned to be this large)
+     integer, allocatable, public                           :: child_phase_id(:)   ! index of each child phase (index into arrays in agsys_phases_type)
+
+     ! TODO(wjs, 2019-11-07) child_phase_name is redundant with information we can already
+     ! get from child_phase_id combined with the phase_name array in agsys_phases_type; we
+     ! may want to keep this for the sake of consistency checking, but we may want to
+     ! remove it.
+     character(len=max_phase_name_len), allocatable, public :: child_phase_name(:) ! name of each child phase
+  end type composite_phase_type
+
+  ! Information on the phases for a given crop
   type, public :: agsys_phases_type
      private
 
      ! Public data members
-     integer, public :: max_phases ! maximum number of phases used by any crop (arrays are dimensioned to be this large)
-     integer, allocatable, public :: num_phases_for_crop(:) ! [crop_type] number of phases used by each crop
-
-     ! Each of these arrays are dimensioned [stage, crop_type]. Stage is defined as a
-     ! point in time (e.g., the time of sowing, or the time of germination). Phase is
-     ! defined as the span of time between two stages (e.g., germinating is the phase
-     ! between sowing and germination). 'phase' variables give the phase whose start point
-     ! is the given stage. For example, if crop type 7 has the first two stages 'sowing'
-     ! and 'germination', then phase_name(1,7) = 'germinating'. Conversely, if
-     ! phase_name(3,7) = 'juvenile', then the name of the stage that triggers the start
-     ! of the juvenile phase is given by stage_name(3,7) (e.g., 'emergence').
-     character(len=max_phase_name_len), allocatable, public :: stage_name(:,:) ! [stage, crop_type]; name of the given stage
-     character(len=max_phase_name_len), allocatable, public :: phase_name(:,:) ! [stage, crop_type]; name of the phase whose start point is the given stage
-     integer                          , allocatable, public :: phase_type(:,:) ! [stage, crop_type]; each value is one of the above phase_type_* constants
+     integer, allocatable, public :: num_phases ! number of phases for this crop (arrays are dimensioned to be this large)
+
+     ! Each of these arrays are dimensioned by phase. Stage is defined as a point in time
+     ! (e.g., the time of sowing, or the time of germination). Phase is defined as the
+     ! span of time between two stages (e.g., germinating is the phase between sowing and
+     ! germination). 'stage' variables give the stage that occurs at the start of the
+     ! given phase. For example, if phase_name(3) = 'juvenile', then the name of the
+     ! stage that triggers the start of the juvenile phase is given by stage_name(3)
+     ! (e.g., 'emergence').
+     character(len=max_phase_name_len), allocatable, public :: stage_name(:) ! name of the stage that occurs at the start of the given phase
+     character(len=max_phase_name_len), allocatable, public :: phase_name(:) ! name of the given phase
+     integer                          , allocatable, public :: phase_type(:) ! each value is one of the above phase_type_* constants
+
+     ! Composite phases for this crop
+     !
+     ! 'composite_phases' has one element for each composite phase defined for this crop.
+     ! 'composite_phase_index_from_type' lets you find the appropriate composite phase
+     ! structure for a given composite phase type (e.g., vernalization).
+     !
+     ! For example, to find the appropriate composite phase structure for vernalization,
+     ! you would do:
+     !
+     !    vernalization_index = phases%composite_phase_index_from_type(composite_phase_type_vernalization)
+     !    phases%composite_phases(vernalization_index)
+     !
+     ! If the given crop doesn't have a vernalization composite phase, then
+     ! phases%composite_phase_index_from_type(composite_phase_type_vernalization) will be
+     ! 0. (This error condition doesn't need to be checked explicitly in the code, because
+     ! it will lead to a runtime error when we test in debug mode.)
+     integer, public :: composite_phase_index_from_type(composite_phase_type_maxval)
+     type(composite_phase_type), allocatable, public :: composite_phases(:)
+
   end type agsys_phases_type
 
 end module AgSysPhases