Merge pull request ESCOMP#10 from billsacks/agsys

Pack environmental inputs into a derived type

src/agsys/ctsm_interface/AgSysInterface.F90

@@ -21,7 +21,7 @@ module AgSysInterface
   use AgSysPhases, only : agsys_phases_type
   use AgSysParamReader, only : ReadParams, ReadPhases
   use AgSysRuntimeConstants, only : InitRuntimeConstants
-  use AgSysPlaceholder, only : DoTimeStep_Phenology_Placeholder, AgsysAbioticStress_Placeholder
+  use AgSysPlaceholder, only : DoTimeStep_Phenology_Placeholder
   !
   implicit none
   private
@@ -109,34 +109,29 @@ subroutine AgSysDriver(this, num_pcropp, filter_pcropp, &
              c = patch%column(p)
              cultivar_type = this%agsys_inst%cultivar_patch(p)
 
-             call AgsysAbioticStress_Placeholder( &
-                  ! Inputs, time-varying
-                  h2osoi_liq_24hr = this%agsys_inst%h2osoi_liq_24hr_col(c, 1:nlevsoi), &
-
-                  ! Outputs
-                  sw_avail_ratio = sw_avail_ratio, &
-                  pesw_seedlayer = pesw_seedlayer)
+             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 DoTimeStep_Phenology_Placeholder( &
                   ! Inputs, time-constant
                   crop            = this%crops(crop_type)%cultivars(cultivar_type), &
 
                   ! Inputs, time-varying
-                  photoperiod    = grc%dayl(g), &
-                  tair_max       = temperature_inst%t_ref2m_max_patch(p), &
-                  tair_min       = temperature_inst%t_ref2m_min_patch(p), &
-                  tc             = temperature_inst%t_veg24_patch(p), &
-                  sw_avail_ratio = sw_avail_ratio, &
-                  pesw_seedlayer = pesw_seedlayer, &
+                  agsys_environmental_inputs = this%agsys_inst%agsys_environmental_inputs, &
 
                   ! 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), &
+                  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,:), &
                   cumvd             = this%agsys_inst%acc_vernalization_days_patch(p))
+
           end if
        end do
     end if

src/agsys/ctsm_interface/AgSysType.F90

@@ -12,10 +12,11 @@ module AgSysType
   use clm_varpar       , only : nlevsoi
   use clm_varcon       , only : spval
   use pftconMod        , only : ntmp_corn, nirrig_tmp_corn, ntrp_corn, nirrig_trp_corn
-  use histFileMod      , only : hist_addfld1d
+  use histFileMod      , only : hist_addfld1d, hist_addfld2d
   use PatchType        , only : patch_type
   use AgSysRuntimeConstants, only : agsys_max_phases
   use AgSysConstants, only : crop_type_not_handled, crop_type_maize
+  use AgSysEnvironmentalInputs, only : agsys_environmental_inputs_type
   !
   implicit none
   private
@@ -51,8 +52,23 @@ module AgSysType
      ! Whether the crop has emerged yet this season
      logical, pointer, public :: emerged_patch(:)
 
-     real(r8), pointer, public :: days_in_phase_patch(:,:)  ! number of days in each phase [phase, patch] (note different dimension order than typical for CTSM)
-     real(r8), pointer, public :: days_after_phase_patch(:,:)  ! number of days after each phase [phase, patch] (note different dimension order than typical for CTSM)
+     ! TODO(wjs, 2019-11-13) 2-d variables are stored with patch (or column) as the first
+     ! dimension. This follows CTSM conventions, but is less efficient for AgSys, where we
+     ! operate on a single point at a time. Ideally, for the sake of performance, we would
+     ! switch this dimension order for 2-d AgSys variables. However, currently that
+     ! prevents us from doing history output for the given 2-d variable. (Also, although
+     ! the restart routines seem to be general enough to handle variables with either
+     ! dimension ordering - via switchdim - it doesn't look like there are any variables
+     ! that currently have switchdim = .false., so some changes may be needed to support
+     ! that robustly, in the restart utilities and/or in init_interp.) At some point, we
+     ! should consider putting in place the necessary generalization for the history
+     ! infrastructure to handle 2-d variables with this switched dimension ordering (as
+     ! long as this can be done in a way that doesn't hurt the efficiency of the history
+     ! infrastructure). At that point, we should change the variables here so that patch
+     ! (or column) is the second dimension.
+
+     real(r8), pointer, public :: days_in_phase_patch(:,:)  ! number of days in each phase [patch, phase]
+     real(r8), pointer, public :: days_after_phase_patch(:,:)  ! number of days after each phase [patch, phase]
 
      ! TODO(wjs, 2019-11-01) We may not need all of these - i.e., it maybe unnecessary to
      ! have both an emerged_thermal_time and thermal time for each phase (since the former
@@ -61,15 +77,18 @@ module AgSysType
      ! example, to store the thermal time just for the previous phase. (If we can avoid
      ! supporting this full generality, that could be good to avoid restart file bloat.)
      real(r8), pointer, public :: acc_emerged_thermal_time_patch(:)  ! accumulated thermal time since emergence (deg-days)
-     real(r8), pointer, public :: acc_thermal_time_in_phase_patch(:,:) ! accumulated thermal time in each phase (deg-days) [phase, patch] (note different dimension order than typical for CTSM)
-     real(r8), pointer, public :: acc_thermal_time_after_phase_patch(:,:) ! accumulated thermal time after each phase (deg-days) [phase, patch] (note different dimension order than typical for CTSM)
+     real(r8), pointer, public :: acc_thermal_time_in_phase_patch(:,:) ! accumulated thermal time in each phase (deg-days) [patch, phase]
+     real(r8), pointer, public :: acc_thermal_time_after_phase_patch(:,:) ! accumulated thermal time after each phase (deg-days) [patch, phase]
 
-     real(r8), pointer, public :: acc_vernalization_days_patch(:) ! accumulated vernalization days (for crops with vernalization) (unit: days) [phase, patch] (note different dimension order than typical for CTSM)
+     real(r8), pointer, public :: acc_vernalization_days_patch(:) ! accumulated vernalization days (for crops with vernalization) (unit: days)
 
      real(r8), pointer, public :: h2osoi_liq_24hr_col(:,:)  ! 24-hour average h2osoi_liq (kg/m2), just over 1:nlevsoi
 
      integer, pointer, public :: days_after_sowing_patch(:)
-
+
+     ! We store an instance of this so that we only need to allocate memory for it once,
+     ! in initialization
+     type(agsys_environmental_inputs_type), public :: agsys_environmental_inputs
 
    contains
      procedure, public :: Init
@@ -131,15 +150,15 @@ subroutine InitAllocate(this, bounds)
     allocate(this%current_stage_patch(begp:endp)); this%current_stage_patch(:) = nan
     allocate(this%emerged_patch(begp:endp)); this%emerged_patch(:) = .false.
 
-    allocate(this%days_in_phase_patch(1:agsys_max_phases, begp:endp))
+    allocate(this%days_in_phase_patch(begp:endp, 1:agsys_max_phases))
     this%days_in_phase_patch(:,:) = nan
-    allocate(this%days_after_phase_patch(1:agsys_max_phases, begp:endp))
+    allocate(this%days_after_phase_patch(begp:endp, 1:agsys_max_phases))
     this%days_after_phase_patch(:,:) = nan
 
     allocate(this%acc_emerged_thermal_time_patch(begp:endp)) ; this%acc_emerged_thermal_time_patch(:) = nan
-    allocate(this%acc_thermal_time_in_phase_patch(1:agsys_max_phases, begp:endp))
+    allocate(this%acc_thermal_time_in_phase_patch(begp:endp, 1:agsys_max_phases))
     this%acc_thermal_time_in_phase_patch(:,:) = nan
-    allocate(this%acc_thermal_time_after_phase_patch(1:agsys_max_phases, begp:endp))
+    allocate(this%acc_thermal_time_after_phase_patch(begp:endp, 1:agsys_max_phases))
     this%acc_thermal_time_after_phase_patch(:,:) = nan
 
     allocate(this%acc_vernalization_days_patch(begp:endp)); this%acc_vernalization_days_patch(:) = nan
@@ -148,6 +167,9 @@ subroutine InitAllocate(this, bounds)
 
     allocate(this%days_after_sowing_patch(begp:endp))  ; this%days_after_sowing_patch(:) = 0
 
+    call this%agsys_environmental_inputs%Init( &
+         nlevsoi = nlevsoi)
+
     end associate
   end subroutine InitAllocate
 
@@ -176,6 +198,12 @@ subroutine InitHistory(this, bounds)
          avgflag='I', long_name='Current phenological stage number (at end of history period)', &
          ptr_patch=this%current_stage_patch)
 
+    this%acc_thermal_time_in_phase_patch(begp:endp,:) = spval
+    call hist_addfld2d(fname='AGSYS_ACC_THERMAL_TIME_IN_PHASE', units='deg-days', &
+         type2d='agsys_phases', &
+         avgflag='I', long_name='Accumulated thermal time in each phase (at end of history period)', &
+         ptr_patch=this%acc_thermal_time_in_phase_patch)
+
     end associate
 
   end subroutine InitHistory
@@ -226,12 +254,12 @@ subroutine InitCold(this, bounds, patch)
     this%current_stage_patch(begp:endp) = 0._r8
     this%emerged_patch(begp:endp) = .false.
 
-    this%days_in_phase_patch(:, begp:endp) = 0._r8
-    this%days_after_phase_patch(:, begp:endp) = 0._r8
+    this%days_in_phase_patch(begp:endp, :) = 0._r8
+    this%days_after_phase_patch(begp:endp, :) = 0._r8
 
     this%acc_emerged_thermal_time_patch(begp:endp) = 0._r8
-    this%acc_thermal_time_in_phase_patch(:, begp:endp) = 0._r8
-    this%acc_thermal_time_after_phase_patch(:, begp:endp) = 0._r8
+    this%acc_thermal_time_in_phase_patch(begp:endp, :) = 0._r8
+    this%acc_thermal_time_after_phase_patch(begp:endp, :) = 0._r8
     this%acc_vernalization_days_patch(:) = 0._r8
 
     ! TODO(wjs, 2019-11-12) We may be able to remove this initialization once we properly

src/agsys/science/AgSysEnvironmentalInputs.F90

@@ -0,0 +1,84 @@
+module AgSysEnvironmentalInputs
+
+  !-----------------------------------------------------------------------
+  ! !DESCRIPTION:
+  ! Derived type holding environmental inputs sent from the host model into AgSys routines
+  !
+  ! The variables here are purely inputs into AgSys (not set by AgSys)
+  !
+  ! !USES:
+  use AgSysKinds, only : r8
+
+  implicit none
+  private
+
+  ! !PUBLIC TYPES:
+
+  type, public :: agsys_environmental_inputs_type
+     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]
+
+   contains
+     procedure, public :: Init       ! Allocate space for this instance (but don't set any values)
+     procedure, public :: SetValues  ! Set values for the current point
+  end type agsys_environmental_inputs_type
+
+contains
+
+  !-----------------------------------------------------------------------
+  subroutine Init(this, nlevsoi)
+    !
+    ! !DESCRIPTION:
+    ! 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.
+    !
+    ! !ARGUMENTS:
+    class(agsys_environmental_inputs_type), intent(inout) :: this
+    integer, intent(in) :: nlevsoi ! number of soil layers
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'Init'
+    !-----------------------------------------------------------------------
+
+    allocate(this%h2osoi_liq_24hr(nlevsoi))
+
+  end subroutine Init
+
+  !-----------------------------------------------------------------------
+  subroutine SetValues(this, photoperiod, tair_max, tair_min, tc_24hr, h2osoi_liq_24hr)
+    !
+    ! !DESCRIPTION:
+    ! Set values for the current point
+    !
+    ! !ARGUMENTS:
+    class(agsys_environmental_inputs_type), intent(inout) :: this
+    real(r8), intent(in) :: photoperiod ! same as day length [h]
+    real(r8), intent(in) :: tair_max  ! daily max air temperature [K]
+    real(r8), intent(in) :: tair_min  ! daily minimum air temperature [K]
+    real(r8), intent(in) :: tc_24hr   ! daily mean canopy temperature [K]
+    real(r8), intent(in) :: h2osoi_liq_24hr(:)  ! daily mean soil liquid content for each soil layer [kg m-2]
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'SetValues'
+    !-----------------------------------------------------------------------
+
+    this%photoperiod = photoperiod
+    this%tair_max = tair_max
+    this%tair_min = tair_min
+    this%tc_24hr = tc_24hr
+    this%h2osoi_liq_24hr(:) = h2osoi_liq_24hr(:)
+
+  end subroutine SetValues
+
+end module AgSysEnvironmentalInputs

src/agsys/science/AgSysPhases.F90

@@ -67,7 +67,7 @@ module AgSysPhases
      private
 
      ! Public data members
-     integer, allocatable, public :: num_phases ! number of phases for this crop (arrays are dimensioned to be this large)
+     integer, 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

src/agsys/science/AgSysPlaceholder.F90

@@ -4,29 +4,24 @@ module AgSysPlaceholder
 
   use AgSysKinds , only : r8
   use AgSysPhases, only : agsys_phases_type
+  use AgSysEnvironmentalInputs, only : agsys_environmental_inputs_type
 
   implicit none
   private
 
   public :: DoTimeStep_Phenology_Placeholder
-  public :: AgsysAbioticStress_Placeholder
 
 contains
 
   subroutine DoTimeStep_Phenology_Placeholder(crop, &
-       photoperiod, tair_max, tair_min, tc, sw_avail_ratio, pesw_seedlayer, &
+       agsys_environmental_inputs, &
        days_after_sowing, current_stage, days_in_phase, tt_in_phase, &
        days_after_phase, tt_after_phase, cumvd)
     ! Inputs, time-constant
     type(agsys_crop_type_generic), intent(in) :: crop
 
     ! Inputs, time-varying
-    real(r8), intent(in) :: photoperiod ! same as day length [h] (exists: grc%dayl)
-    real(r8), intent(in) :: tair_max  ! daily max [K] (exists, at least by end of driver loop)
-    real(r8), intent(in) :: tair_min  ! daily minimum [K] (exists)
-    real(r8), intent(in) :: tc        ! daily mean canopy temperature [K] (exists: t_veg24_patch)
-    real(r8), intent(in) :: sw_avail_ratio ! soil water available ratio; this will be calculated earlier by a different AgSys routine (which needs daily average h2osoi_liq in each column, and time-constant soil properties)
-    real(r8), intent(in) :: pesw_seedlayer ! see comment for sw_avail_ratio
+    type(agsys_environmental_inputs_type), intent(in) :: agsys_environmental_inputs
 
     ! Outputs
     integer, intent(inout) :: days_after_sowing
@@ -38,13 +33,4 @@ subroutine DoTimeStep_Phenology_Placeholder(crop, &
     real(r8), intent(inout) :: cumvd ! cumulative vernalization days (ignored for crops without vernalization)
   end subroutine DoTimeStep_Phenology_Placeholder
 
-  subroutine AgsysAbioticStress_Placeholder(h2osoi_liq_24hr, &
-       ! And some other inputs (soil parameters, etc.)
-       sw_avail_ratio, pesw_seedlayer)
-    real(r8), intent(in) :: h2osoi_liq_24hr(:)  ! 1 .. num_soil_layers
-    ! And some other inputs (soil parameters, etc.)
-    real(r8), intent(out) :: sw_avail_ratio
-    real(r8), intent(out) :: pesw_seedlayer
-  end subroutine AgsysAbioticStress_Placeholder
-
 end module AgSysPlaceholder

src/main/clm_instMod.F90

@@ -9,7 +9,7 @@ module clm_instMod
   use decompMod       , only : bounds_type
   use clm_varpar      , only : ndecomp_pools, nlevdecomp_full
   use clm_varctl      , only : use_cn, use_c13, use_c14, use_lch4, use_cndv, use_fates
-  use clm_varctl      , only : use_century_decomp, use_crop, snow_cover_fraction_method, paramfile
+  use clm_varctl      , only : use_century_decomp, use_crop, use_crop_agsys, snow_cover_fraction_method, paramfile
   use clm_varcon      , only : bdsno, c13ratio, c14ratio
   use landunit_varcon , only : istice_mec, istsoil
   use perf_mod        , only : t_startf, t_stopf
@@ -438,6 +438,10 @@ subroutine clm_instInit(bounds)
        call clm_fates%Init(bounds)
     end if
 
+    if (use_crop_agsys) then
+       call agsys_interface_inst%Init(bounds, patch)
+    end if
+
     deallocate (h2osno_col)
     deallocate (snow_depth_col)