Merge pull request ESCOMP#11 from billsacks/agsys

Latest agsys changes

src/agsys/ctsm_interface/AgSysInterface.F90

@@ -9,7 +9,7 @@ module AgSysInterface
   ! !USES:
 #include "shr_assert.h"
   use shr_kind_mod    , only : r8 => shr_kind_r8
-  use clm_time_manager, only : is_beg_curr_day
+  use clm_time_manager, only : is_beg_curr_day, get_curr_calday
   use decompMod, only : bounds_type
   use clm_varpar, only : nlevsoi
   use PatchType, only : patch_type
@@ -20,7 +20,6 @@ module AgSysInterface
   use AgSysPhases, only : agsys_phases_type
   use AgSysParamReader, only : ReadParams
   use AgSysRuntimeConstants, only : InitRuntimeConstants
-  use AgSysPlaceholder, only : DoTimeStep_Phenology_Placeholder
   use AgSysCropTypeGeneric, only : agsys_cultivars_of_crop_type
   use AgSysRoot, only : get_soil_condition
   use AgSysPhenology, only : AgSysRunPhenology
@@ -102,6 +101,9 @@ subroutine AgSysDriver(this, num_pcropp, filter_pcropp, &
     ! temperature accumulators, because the current routine is called before the various
     ! UpdateAccVars calls in the driver loop.)
     if (is_beg_curr_day()) then
+       call this%agsys_inst%agsys_environmental_inputs%SetSpatiallyConstantValues( &
+            calday = floor(get_curr_calday()))
+
        do fp = 1, num_pcropp
           p = filter_pcropp(fp)
 
@@ -111,38 +113,42 @@ subroutine AgSysDriver(this, num_pcropp, filter_pcropp, &
              c = patch%column(p)
              cultivar_type = this%agsys_inst%cultivar_patch(p)
 
-             call this%agsys_inst%agsys_environmental_inputs%SetValues( &
-                  photoperiod    = grc%dayl(g), &
-                  tair_max       = temperature_inst%t_ref2m_max_patch(p), &
-                  tair_min       = temperature_inst%t_ref2m_min_patch(p), &
-                  tc_24hr        = temperature_inst%t_veg24_patch(p), &
-                  h2osoi_liq_24hr = this%agsys_inst%h2osoi_liq_24hr_col(c, 1:nlevsoi))
-
-             call get_soil_condition( &
-                  crop      = this%crops(crop_type)%cultivars(cultivar_type), &
-                  env       = this%agsys_inst%agsys_environmental_inputs, &
-                  soil_prop = this%agsys_inst%agsys_soil_properties, &
-                  root      = this%agsys_inst%agsys_root_properties, &
-                  soil_cond  = this%agsys_inst%agsys_soil_condition)
-
-             !call DoTimeStep_Phenology_Placeholder( &
-             call AgSysRunPhenology ( &
-                  ! Inputs, time-constant
-                  crop      = this%crops(crop_type)%cultivars(cultivar_type), &
-
-                  ! Inputs, time-varying
-                  env       = this%agsys_inst%agsys_environmental_inputs, &
-                  soil_cond = this%agsys_inst%agsys_soil_condition, &
-
-                  ! Outputs
-                  days_after_sowing = this%agsys_inst%days_after_sowing_patch(p), &
-                  current_stage     = this%agsys_inst%current_stage_patch(p), &
-                  days_in_phase     = this%agsys_inst%days_in_phase_patch(p,:), &
-                  tt_in_phase       = this%agsys_inst%acc_thermal_time_in_phase_patch(p,:), &
-                  days_after_phase  = this%agsys_inst%days_after_phase_patch(p,:), &
-                  tt_after_phase    = this%agsys_inst%acc_thermal_time_after_phase_patch(p,:), &
-                  phase_target_tt   = this%agsys_inst%phase_target_thermal_time_path(p,:), &
-                  cumvd             = this%agsys_inst%acc_vernalization_days_patch(p))
+
+
+             if (this%agsys_inst%crop_alive_patch(p)) then
+                call this%agsys_inst%agsys_environmental_inputs%SetSpatiallyVaryingValues( &
+                     photoperiod    = grc%dayl(g), &
+                     tair_max       = temperature_inst%t_ref2m_max_patch(p), &
+                     tair_min       = temperature_inst%t_ref2m_min_patch(p), &
+                     tc_24hr        = temperature_inst%t_veg24_patch(p), &
+                     h2osoi_liq_24hr = this%agsys_inst%h2osoi_liq_24hr_col(c, 1:nlevsoi))
+
+                call get_soil_condition( &
+                     crop      = this%crops(crop_type)%cultivars(cultivar_type), &
+                     env       = this%agsys_inst%agsys_environmental_inputs, &
+                     soil_prop = this%agsys_inst%agsys_soil_properties, &
+                     root      = this%agsys_inst%agsys_root_properties, &
+                     soil_cond  = this%agsys_inst%agsys_soil_condition)
+
+                call AgSysRunPhenology ( &
+                     ! Inputs, time-constant
+                     crop      = this%crops(crop_type)%cultivars(cultivar_type), &
+
+                     ! Inputs, time-varying
+                     env       = this%agsys_inst%agsys_environmental_inputs, &
+                     soil_cond = this%agsys_inst%agsys_soil_condition, &
+
+                     ! Outputs
+                     crop_alive        = this%agsys_inst%crop_alive_patch(p), &
+                     days_after_sowing = this%agsys_inst%days_after_sowing_patch(p), &
+                     current_stage     = this%agsys_inst%current_stage_patch(p), &
+                     days_in_phase     = this%agsys_inst%days_in_phase_patch(p,:), &
+                     tt_in_phase       = this%agsys_inst%acc_thermal_time_in_phase_patch(p,:), &
+                     days_after_phase  = this%agsys_inst%days_after_phase_patch(p,:), &
+                     tt_after_phase    = this%agsys_inst%acc_thermal_time_after_phase_patch(p,:), &
+                     phase_target_tt   = this%agsys_inst%phase_target_thermal_time_path(p,:), &
+                     cumvd             = this%agsys_inst%acc_vernalization_days_patch(p))
+             end if
 
           end if
        end do
@@ -171,7 +177,7 @@ subroutine Init(this, bounds, patch)
     !-----------------------------------------------------------------------
 
     call ReadParams(this%crops)
-    !call InitRuntimeConstants(this%crops)
+    call InitRuntimeConstants(this%crops)
     call this%agsys_inst%Init(bounds, patch)
 
   end subroutine Init

src/agsys/ctsm_interface/AgSysParamReader.F90

@@ -32,16 +32,17 @@ subroutine ReadParams(crops)
     !
     ! !ARGUMENTS:
     type(agsys_cultivars_of_crop_type), intent(inout) :: crops(:)
-    class(agsys_crop_type_generic), allocatable :: cultivar
     !
     ! !LOCAL VARIABLES:
+    class(agsys_crop_type_generic), allocatable :: cultivar
+    integer :: crop_type
 
     character(len=*), parameter :: subname = 'ReadParams'
     !-----------------------------------------------------------------------
 
     SHR_ASSERT_FL((size(crops) == crop_type_maxval), sourcefile, __LINE__)
 
-    allocate(agsys_crop_type_photosensitive :: crops(crop_type_maize)%cultivars(1))
+    allocate(agsys_crop_type_maize :: crops(crop_type_maize)%cultivars(1))
 
     !!currently we hard-coded parameters for one maize cultivar (Pioneer_P04612XR_106)
     !!check the parameter values in the Maize.xml of APSIM
@@ -85,7 +86,16 @@ subroutine ReadParams(crops)
           cultivar%leaf_no_dead_const = -0.025_r8
           cultivar%leaf_no_dead_slope = 0.00035_r8      
       end select
-    end associate
+      end associate
+
+      ! TODO(wjs, 2019-11-15) Set other crops / cultivars. For now just initialize to
+      ! default values.
+      do crop_type = 1, crop_type_maxval
+         if (crop_type /= crop_type_maize) then
+            allocate(agsys_crop_type_generic :: crops(crop_type)%cultivars(1))
+            call crops(crop_type)%cultivars(1)%init()
+         end if
+      end do
 
   end subroutine ReadParams
 

src/agsys/ctsm_interface/AgSysRuntimeConstants.F90

@@ -10,7 +10,7 @@ module AgSysRuntimeConstants
   ! dimension sizes on restart and/or history files).
   !
   ! !USES:
-  use AgSysPhases, only : agsys_phases_type
+  use AgSysCropTypeGeneric,         only : agsys_cultivars_of_crop_type
   !
   implicit none
   private
@@ -28,24 +28,30 @@ module AgSysRuntimeConstants
 contains
 
   !-----------------------------------------------------------------------
-  subroutine InitRuntimeConstants(crop_phases)
+  subroutine InitRuntimeConstants(crops)
     !
     ! !DESCRIPTION:
     ! Initialize runtime constants in this module
     !
     ! !ARGUMENTS:
-    type(agsys_phases_type), intent(in) :: crop_phases(:) ! phases for each crop
+    type(agsys_cultivars_of_crop_type), intent(in) :: crops(:)
     !
     ! !LOCAL VARIABLES:
-    integer :: crop
+    integer :: crop_type
+    integer :: cultivar
+    integer :: num_phases_this_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
+    do crop_type = 1, ubound(crops,1)
+       ! All cultivars of a given crop type have the same number of phases, so just take
+       ! information from the first cultivar of each crop type.
+       cultivar = 1
+       num_phases_this_crop = crops(crop_type)%cultivars(1)%phases%num_phases
+       if (num_phases_this_crop > agsys_max_phases) then
+          agsys_max_phases = num_phases_this_crop
        end if
     end do
 

src/agsys/ctsm_interface/AgSysType.F90

@@ -50,6 +50,9 @@ module AgSysType
      ! how far we are from one stage to the next. Before sowing, this is 0.
      real(r8), pointer, public :: current_stage_patch(:)
 
+     ! Whether the crop is currently alive; true in the entire period from sowing to harvest
+     logical, pointer, public :: crop_alive_patch(:)
+
      ! Whether the crop has emerged yet this season
      logical, pointer, public :: emerged_patch(:)
 
@@ -151,6 +154,7 @@ subroutine InitAllocate(this, bounds)
     allocate(this%crop_type_patch(begp:endp)); this%crop_type_patch(:) = 0
     allocate(this%cultivar_patch(begp:endp)); this%cultivar_patch(:) = 0
     allocate(this%current_stage_patch(begp:endp)); this%current_stage_patch(:) = nan
+    allocate(this%crop_alive_patch(begp:endp)); this%crop_alive_patch(:) = .false.
     allocate(this%emerged_patch(begp:endp)); this%emerged_patch(:) = .false.
 
     allocate(this%days_in_phase_patch(begp:endp, 1:agsys_max_phases))
@@ -255,6 +259,7 @@ subroutine InitCold(this, bounds, patch)
     end do
 
     this%current_stage_patch(begp:endp) = 0._r8
+    this%crop_alive_patch(begp:endp) = .false.
     this%emerged_patch(begp:endp) = .false.
 
     this%days_in_phase_patch(begp:endp, :) = 0._r8

src/agsys/science/AgSysCropTypeGeneric.F90

@@ -7,18 +7,6 @@ module AgSysCropTypeGeneric
 
   implicit none
 
-  !!a container type to store the crop type instances
-  type, public :: agsys_crop_container_type
-     private
-     class(agsys_crop_type_generic), allocatable, public :: crop
-  end type agsys_crop_container_type
-
-  !!a container type to store all the cultivars of a given crop type
-  type, public :: agsys_cultivars_of_crop_type
-     private
-     class(agsys_crop_type_generic), allocatable, public :: cultivars(:)
-  end type agsys_cultivars_of_crop_type
-
   !-----------------------------------------------------------------
   !!a generic crop type holding parameters that shared by all crops
   type, public :: agsys_crop_type_generic
@@ -52,6 +40,12 @@ module AgSysCropTypeGeneric
      procedure :: get_target_tt_inductive_phase
   end type agsys_crop_type_generic
 
+  !!a container type to store all the cultivars of a given crop type
+  type, public :: agsys_cultivars_of_crop_type
+     private
+     class(agsys_crop_type_generic), allocatable, public :: cultivars(:)
+  end type agsys_cultivars_of_crop_type
+
 contains
   subroutine init(this)
     class(agsys_crop_type_generic), intent(inout) :: this
@@ -62,6 +56,9 @@ subroutine init(this)
     this%rc_tair_tt%num_pts = 0
     this%rc_sw_avail_phenol%num_pts = 0
     this%rc_sw_emerg_rate%num_pts = 0
+
+    ! This will be set later, by the specific child class
+    this%phases%num_phases = 0
   end subroutine init
 
   subroutine vernalization(this, env, cumvd)

src/agsys/science/AgSysCropTypeMaize.F90

@@ -30,6 +30,9 @@ module AgSysCropTypeMaize
   subroutine init(this)
     class(agsys_crop_type_maize), intent(inout) :: this
 
+    ! Initialize the parent class
+    call this%agsys_crop_type_photosensitive%init()
+
     this%croptype                   = crop_type_maize
 
     !!!initialize the parameters

src/agsys/science/AgSysCropTypePhotoSensitive.F90

@@ -16,7 +16,11 @@ module AgSysCropTypePhotoSensitive
 contains
   subroutine init(this)
     class(agsys_crop_type_photosensitive), intent(inout) :: this
-    this%rc_photoperiod_target_tt%num_pts = 0 
+
+    ! Initialize the parent class
+    call this%agsys_crop_type_generic%init()
+
+    this%rc_photoperiod_target_tt%num_pts = 0
   end subroutine init
 
   function get_target_tt_photosensitive_phase(this, env) result(target_tt)

src/agsys/science/AgSysEnvironmentalInputs.F90

@@ -18,15 +18,17 @@ module AgSysEnvironmentalInputs
      private
 
      ! Public data members
-     real(r8), public :: photoperiod ! same as day length [h]
-     real(r8), public :: tair_max  ! daily max air temperature [K]
-     real(r8), public :: tair_min  ! daily minimum air temperature [K]
-     real(r8), public :: tc_24hr   ! daily mean canopy temperature [K]
-     real(r8), allocatable, public :: h2osoi_liq_24hr(:)  ! daily mean soil liquid content for each soil layer [kg m-2]
+     integer , public :: calday                          ! calendar day; 1 = Jan 1
+     real(r8), public :: photoperiod                     ! same as day length [h]
+     real(r8), public :: tair_max                        ! daily max air temperature [K]
+     real(r8), public :: tair_min                        ! daily minimum air temperature [K]
+     real(r8), public :: tc_24hr                         ! daily mean canopy temperature [K]
+     real(r8), allocatable, public :: h2osoi_liq_24hr(:) ! daily mean soil liquid content for each soil layer [kg m-2]
 
    contains
      procedure, public :: Init       ! Allocate space for this instance (but don't set any values)
-     procedure, public :: SetValues  ! Set values for the current point
+     procedure, public :: SetSpatiallyConstantValues ! Set spatially-constant values for this time
+     procedure, public :: SetSpatiallyVaryingValues  ! Set values for the current point for this time
   end type agsys_environmental_inputs_type
 
 contains
@@ -38,8 +40,9 @@ subroutine Init(this, nlevsoi)
     ! Allocate space for this instance (but don't set any values)
     !
     ! This should be called once, in initialization. The purpose of separating this from
-    ! SetValues is so that we can just do the memory allocation once, rather than doing
-    ! this memory allocation repeatedly for every time step and every point.
+    ! SetSpatiallyVaryingValues is so that we can just do the memory allocation once,
+    ! rather than doing this memory allocation repeatedly for every time step and every
+    ! point.
     !
     ! !ARGUMENTS:
     class(agsys_environmental_inputs_type), intent(inout) :: this
@@ -55,7 +58,26 @@ subroutine Init(this, nlevsoi)
   end subroutine Init
 
   !-----------------------------------------------------------------------
-  subroutine SetValues(this, photoperiod, tair_max, tair_min, tc_24hr, h2osoi_liq_24hr)
+  subroutine SetSpatiallyConstantValues(this, calday)
+    !
+    ! !DESCRIPTION:
+    ! Set spatially-constant values for this time
+    !
+    ! !ARGUMENTS:
+    class(agsys_environmental_inputs_type), intent(inout) :: this
+    integer, intent(in) :: calday  ! calendar day; 1 = Jan 1
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'SetSpatiallyConstantValues'
+    !-----------------------------------------------------------------------
+
+    this%calday = calday
+
+  end subroutine SetSpatiallyConstantValues
+
+  !-----------------------------------------------------------------------
+  subroutine SetSpatiallyVaryingValues(this, photoperiod, tair_max, tair_min, tc_24hr, h2osoi_liq_24hr)
     !
     ! !DESCRIPTION:
     ! Set values for the current point
@@ -70,7 +92,7 @@ subroutine SetValues(this, photoperiod, tair_max, tair_min, tc_24hr, h2osoi_liq_
     !
     ! !LOCAL VARIABLES:
 
-    character(len=*), parameter :: subname = 'SetValues'
+    character(len=*), parameter :: subname = 'SetSpatiallyVaryingValues'
     !-----------------------------------------------------------------------
 
     this%photoperiod = photoperiod
@@ -79,6 +101,6 @@ subroutine SetValues(this, photoperiod, tair_max, tair_min, tc_24hr, h2osoi_liq_
     this%tc_24hr = tc_24hr
     this%h2osoi_liq_24hr(:) = h2osoi_liq_24hr(:)
 
-  end subroutine SetValues
+  end subroutine SetSpatiallyVaryingValues
 
 end module AgSysEnvironmentalInputs

src/agsys/science/AgSysPhenology.F90

@@ -16,12 +16,17 @@ module AgSysPhenology
   use AgSysExcepUtils,          only : iulog, endrun 
 
   implicit none
+  private
+
+  public :: AgSysRunPhenology
+
 contains
   !---------------------------------------------------------------
   !Some subroutines
   !
+
   subroutine AgSysRunPhenology(crop, env, soil_cond, &
-                               days_after_sowing, current_stage, days_in_phase, tt_in_phase, &
+                               crop_alive, days_after_sowing, current_stage, days_in_phase, tt_in_phase, &
                                days_after_phase, tt_after_phase, phase_target_tt, cumvd)
 
     !!---------------------------------------------------------
@@ -35,6 +40,7 @@ subroutine AgSysRunPhenology(crop, env, soil_cond, &
     type (agsys_soil_condition_type),         intent(in) :: soil_cond
 
     !!OUTPUTS: state variables
+    logical,  intent(inout) :: crop_alive         !!whether the crop is alive (true between planting and sowing)
     integer,  intent(inout) :: days_after_sowing  !!days after sowing, this is an accumulated number since sowing
     real(r8), intent(inout) :: current_stage      !!current stage number
     real(r8), intent(inout) :: days_in_phase(:)   !!days since start of the current phase