Merge pull request ESCOMP#13 from billsacks/agsys

Latest changes

src/agsys/ctsm_interface/AgSysInterface.F90

@@ -53,6 +53,9 @@ module AgSysInterface
    contains
      procedure, public :: AgSysDriver
      procedure, public :: Init
+     procedure, public :: InitAccBuffer
+     procedure, public :: InitAccVars
+     procedure, public :: UpdateAccVars
   end type agsys_interface_type
 
   character(len=*), parameter, private :: sourcefile = &
@@ -92,14 +95,19 @@ subroutine AgSysDriver(this, num_pcropp, filter_pcropp, &
     character(len=*), parameter :: subname = 'AgSysDriver'
     !-----------------------------------------------------------------------
 
-    ! TODO(wjs, 2019-11-08) Do climate stuff. This will involve updating AgSys-specific accumulator fields.
-
     ! Note that 24-hour temperature accumulators (which are needed here) are updated late
     ! in the driver loop, based on end_curr_day. We use beg_curr_day here so that the
     ! phenology routines are called the time step after these temperature accumulators
     ! are updated. (If we used end_curr_day here, we wouldn't have today's updated
     ! temperature accumulators, because the current routine is called before the various
     ! UpdateAccVars calls in the driver loop.)
+    !
+    ! TODO(wjs, 2019-11-15) We might want extra logic here preventing this from running
+    ! if the model has run less than a day total, because various daily averages wouldn't
+    ! be available yet in that case. That might be important if a run starts mid-day, and
+    ! may also be important to avoid weirdness on the first time step of the simulation.
+    ! (Is there a clm time manager routine that lets you check the elapsed time into the
+    ! simulation that we could use for this purpose?)
     if (is_beg_curr_day()) then
        call this%agsys_inst%agsys_environmental_inputs%SetSpatiallyConstantValues( &
             calday = floor(get_curr_calday()))
@@ -117,7 +125,7 @@ subroutine AgSysDriver(this, num_pcropp, filter_pcropp, &
                   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), &
+                  tc_24hr        = this%agsys_inst%t_veg24hr_patch(p), &
                   h2osoi_liq_24hr = this%agsys_inst%h2osoi_liq_24hr_col(c, 1:nlevsoi))
 
              if (.not. this%agsys_inst%crop_alive_patch(p)) then
@@ -192,4 +200,61 @@ subroutine Init(this, bounds, patch)
 
   end subroutine Init
 
+  !-----------------------------------------------------------------------
+  subroutine InitAccBuffer(this, bounds)
+    !
+    ! !DESCRIPTION:
+    ! Initialize accumulation buffer for all AgSys variables
+    !
+    ! !ARGUMENTS:
+    class(agsys_interface_type), intent(in) :: this
+    type(bounds_type), intent(in) :: bounds
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'InitAccBuffer'
+    !-----------------------------------------------------------------------
+
+    call this%agsys_inst%InitAccBuffer(bounds)
+
+  end subroutine InitAccBuffer
+
+  !-----------------------------------------------------------------------
+  subroutine InitAccVars(this, bounds)
+    !
+    ! !DESCRIPTION:
+    ! Initialize variables that are associated with accumulated fields
+    !
+    ! !ARGUMENTS:
+    class(agsys_interface_type), intent(inout) :: this
+    type(bounds_type), intent(in) :: bounds
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'InitAccVars'
+    !-----------------------------------------------------------------------
+
+    call this%agsys_inst%InitAccVars(bounds)
+
+  end subroutine InitAccVars
+
+  !-----------------------------------------------------------------------
+  subroutine UpdateAccVars(this, bounds)
+    !
+    ! !DESCRIPTION:
+    ! Update accumulated variables
+    !
+    ! !ARGUMENTS:
+    class(agsys_interface_type), intent(inout) :: this
+    type(bounds_type), intent(in) :: bounds
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'UpdateAccVars'
+    !-----------------------------------------------------------------------
+
+    call this%agsys_inst%UpdateAccVars(bounds)
+
+  end subroutine UpdateAccVars
+
 end module AgSysInterface

src/agsys/ctsm_interface/AgSysParamReader.F90

@@ -52,49 +52,51 @@ subroutine ReadParams(crops)
       ! can access params like this:
       ! cultivar_params%shoot_lag
       call cultivar%init()
+      allocate(cultivar%max_days_from_sowing_to_end_of_phase(11))
+      cultivar%max_days_from_sowing_to_end_of_phase(:) = 0
+      cultivar%p_sowing_depth = 50._r8  ! [mm], this is information about management and can be put
+      ! into a separate structure later [TODO (2019-11-12, pb)]
+      cultivar%p_shoot_lag    = 15._r8  ! [degree-days]
+      cultivar%p_shoot_rate   = 0.6_r8  ! [degree-days per mm depth]
+      cultivar%p_pesw_germ    = 0._r8   ! soil moisture threshold for seed germination
+
+      cultivar%rc_tair_tt%x=[0._r8, 18._r8, 26._r8, 34._r8, 44._r8]
+      cultivar%rc_tair_tt%y=[0._r8, 10._r8, 18._r8, 26._r8, 0._r8]
+      cultivar%rc_tair_tt%num_pts=5
+
+      cultivar%rc_sw_avail_phenol%x=[0._r8, 0.16_r8]
+      cultivar%rc_sw_avail_phenol%y=[0._r8, 1._r8]
+      cultivar%rc_sw_avail_phenol%num_pts=2
+
+      cultivar%rc_sw_emerg_rate%x=[0._r8, 1._r8]
+      cultivar%rc_sw_emerg_rate%y=[1._r8, 1._r8] !for maize, this is deactivated
+      cultivar%rc_sw_emerg_rate%num_pts=2
       select type(cultivar)
-        class is (agsys_crop_type_generic)
-          cultivar%p_sowing_depth = 50._r8  ! [mm], this is information about management and can be put 
-                                            ! into a separate structure later [TODO (2019-11-12, pb)]
-          cultivar%p_shoot_lag    = 15._r8  ! [degree-days]
-          cultivar%p_shoot_rate   = 0.6_r8  ! [degree-days per mm depth]
-          cultivar%p_pesw_germ    = 0._r8   ! soil moisture threshold for seed germination
-
-          cultivar%rc_tair_tt%x=[0._r8, 18._r8, 26._r8, 34._r8, 44._r8]
-          cultivar%rc_tair_tt%y=[0._r8, 10._r8, 18._r8, 26._r8, 0._r8]
-          cultivar%rc_tair_tt%num_pts=5
-
-          cultivar%rc_sw_avail_phenol%x=[0._r8, 0.16_r8]
-          cultivar%rc_sw_avail_phenol%y=[0._r8, 1._r8]
-          cultivar%rc_sw_avail_phenol%num_pts=2
-
-          cultivar%rc_sw_emerg_rate%x=[0._r8, 1._r8]
-          cultivar%rc_sw_emerg_rate%y=[1._r8, 1._r8] !for maize, this is deactivated
-          cultivar%rc_sw_emerg_rate%num_pts=2
-
-        class is (agsys_crop_type_photosensitive)
-          cultivar%rc_photoperiod_target_tt%x=[0._r8, 12.5_r8, 20._r8]
-          cultivar%rc_photoperiod_target_tt%y=[0._r8, 0._r8, 0._r8]
-          cultivar%rc_photoperiod_target_tt%num_pts=3
-        class is (agsys_crop_type_maize)
-          cultivar%target_tt_emerg_to_endjuv = 275._r8
-          cultivar%target_tt_flag_to_flower  = 1._r8
-          cultivar%target_tt_flower_to_maturity = 812._r8
-          cultivar%target_tt_flower_to_start_grain = 170._r8
-          cultivar%target_tt_maturity_to_ripe = 1._r8
-
-          cultivar%leaf_init_rate = 23.2_r8
-          cultivar%leaf_no_seed = 6._r8
-          cultivar%leaf_no_at_emerg = 0.5_r8
-          cultivar%leaf_no_min = 5.0_r8
-          cultivar%leaf_no_max = 40._r8
-          cultivar%leaf_no_critical = 11._r8
-          cultivar%leaf_appearance_rate_early = 65._r8
-          cultivar%leaf_appearance_rate_late = 36._r8
-
-          cultivar%potential_kernel_weight = 300._r8
-          cultivar%leaf_no_dead_const = -0.025_r8
-          cultivar%leaf_no_dead_slope = 0.00035_r8      
+      class is (agsys_crop_type_photosensitive)
+         cultivar%rc_photoperiod_target_tt%x=[0._r8, 12.5_r8, 20._r8]
+         cultivar%rc_photoperiod_target_tt%y=[0._r8, 0._r8, 0._r8]
+         cultivar%rc_photoperiod_target_tt%num_pts=3
+      end select
+      select type(cultivar)
+      class is (agsys_crop_type_maize)
+         cultivar%target_tt_emerg_to_endjuv = 275._r8
+         cultivar%target_tt_flag_to_flower  = 1._r8
+         cultivar%target_tt_flower_to_maturity = 812._r8
+         cultivar%target_tt_flower_to_start_grain = 170._r8
+         cultivar%target_tt_maturity_to_ripe = 1._r8
+
+         cultivar%leaf_init_rate = 23.2_r8
+         cultivar%leaf_no_seed = 6._r8
+         cultivar%leaf_no_at_emerg = 0.5_r8
+         cultivar%leaf_no_min = 5.0_r8
+         cultivar%leaf_no_max = 40._r8
+         cultivar%leaf_no_critical = 11._r8
+         cultivar%leaf_appearance_rate_early = 65._r8
+         cultivar%leaf_appearance_rate_late = 36._r8
+
+         cultivar%potential_kernel_weight = 300._r8
+         cultivar%leaf_no_dead_const = -0.025_r8
+         cultivar%leaf_no_dead_slope = 0.00035_r8      
       end select
       end associate
 

src/agsys/ctsm_interface/AgSysType.F90

@@ -9,10 +9,12 @@ module AgSysType
   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 clm_varpar       , only : nlevsoi
-  use clm_varcon       , only : spval
+  use clm_varpar       , only : nlevgrnd, nlevsoi
+  use clm_varcon       , only : spval, tfrz
   use pftconMod        , only : ntmp_corn, nirrig_tmp_corn, ntrp_corn, nirrig_trp_corn
   use histFileMod      , only : hist_addfld1d, hist_addfld2d
+  use clm_time_manager , only : get_nstep
+  use accumulMod       , only : init_accum_field, extract_accum_field, update_accum_field
   use PatchType        , only : patch_type
   use AgSysRuntimeConstants,    only : agsys_max_phases
   use AgSysConstants,           only : crop_type_not_handled, crop_type_maize
@@ -86,7 +88,12 @@ module AgSysType
      real(r8), pointer, public :: phase_target_thermal_time_patch(:,:)  !target thermal time to finish a phase (deg-days) [path, phase]
      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
+     real(r8), pointer, public :: t_veg24hr_patch(:)  ! 24-hour average vegetation (canopy) temperature (K)
+
+     ! TODO(wjs, 2019-11-15) Consider changing this to only be over 1:nlevsoi, rather than
+     ! going all the way up to nlevgrnd. This will require adding handling of variables
+     ! with nlevsoi dimension to the restart file and init_interp multilevel.
+     real(r8), pointer, public :: h2osoi_liq_24hr_col(:,:)  ! 24-hour average h2osoi_liq (kg/m2), over 1:nlevgrnd
 
      integer, pointer, public :: days_after_sowing_patch(:)
 
@@ -101,6 +108,10 @@ module AgSysType
      procedure, private :: InitAllocate
      procedure, private :: InitHistory
      procedure, private :: InitCold
+
+     procedure, public :: InitAccBuffer
+     procedure, public :: InitAccVars
+     procedure, public :: UpdateAccVars
   end type agsys_type
 
   character(len=*), parameter, private :: sourcefile = &
@@ -167,10 +178,13 @@ subroutine InitAllocate(this, bounds)
     this%acc_thermal_time_in_phase_patch(:,:) = nan
     allocate(this%acc_thermal_time_after_phase_patch(begp:endp, 1:agsys_max_phases))
     this%acc_thermal_time_after_phase_patch(:,:) = nan
+    allocate(this%phase_target_thermal_time_patch(begp:endp, 1:agsys_max_phases))
+    this%phase_target_thermal_time_patch(:,:) = nan
 
     allocate(this%acc_vernalization_days_patch(begp:endp)); this%acc_vernalization_days_patch(:) = nan
 
-    allocate(this%h2osoi_liq_24hr_col(begc:endc, 1:nlevsoi)); this%h2osoi_liq_24hr_col(:,:) = nan
+    allocate(this%t_veg24hr_patch(begp:endp)); this%t_veg24hr_patch(:) = nan
+    allocate(this%h2osoi_liq_24hr_col(begc:endc, 1:nlevgrnd)); this%h2osoi_liq_24hr_col(:,:) = nan
 
     allocate(this%days_after_sowing_patch(begp:endp))  ; this%days_after_sowing_patch(:) = 0
 
@@ -247,6 +261,13 @@ subroutine InitCold(this, bounds, patch)
        if ( ptype == ntmp_corn .or. ptype == nirrig_tmp_corn .or. &
             ptype == ntrp_corn .or. ptype == nirrig_trp_corn) then
           this%crop_type_patch(p) = crop_type_maize
+
+          ! TODO(wjs, 2019-11-15) For now, for variables output to history file, just
+          ! initialize over the crops that we're handling. Eventually we should
+          ! initialize these for all crops.
+          this%current_stage_patch(p) = 0._r8
+          this%acc_thermal_time_in_phase_patch(p, :) = 0._r8
+
           ! TODO(wjs, 2019-11-12) Handle more crop types here
        else
           this%crop_type_patch(p) = crop_type_not_handled
@@ -258,26 +279,125 @@ subroutine InitCold(this, bounds, patch)
        end associate
     end do
 
-    this%current_stage_patch(begp:endp) = 0._r8
+    ! ------------------------------------------------------------------------
+    ! For variables not output to history file, initialize over all points
+    !
+    ! TODO(wjs, 2019-11-15) We may just want to initialize these over crop points. It
+    ! would be best if all variable initialization happened the same way here, regardless
+    ! of whether the variable is output to history files.
+    ! ------------------------------------------------------------------------
+
     this%crop_alive_patch(begp:endp) = .false.
     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%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_vernalization_days_patch(:) = 0._r8
 
-    ! TODO(wjs, 2019-11-12) We may be able to remove this initialization once we properly
-    ! initialize the accumulator field related to this variable
-    this%h2osoi_liq_24hr_col(begc:endc, :) = 0._r8
-
     this%days_after_sowing_patch(begp:endp) = 0
 
     end associate
 
   end subroutine InitCold
 
+  !-----------------------------------------------------------------------
+  subroutine InitAccBuffer(this, bounds)
+    !
+    ! !DESCRIPTION:
+    ! Initialize accumulation buffer for all required accumulated fields
+    !
+    ! !ARGUMENTS:
+    class(agsys_type), intent(in) :: this
+    type(bounds_type), intent(in) :: bounds
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'InitAccBuffer'
+    !-----------------------------------------------------------------------
+
+    ! TODO(wjs, 2019-11-14) Change these 24-hour accumulators to instead be end-of-day
+    ! daily average accumulation once that is implemented.
+    !
+    ! BUG(wjs, 2019-11-14, ESCOMP/ctsm#839) Waiting on implementation of end-of-day daily
+    ! average accumulation
+
+    ! NOTE(wjs, 2019-11-15) This differs from t_veg24_patch in temperature_type in that
+    ! it uses timeavg rather than runmean interpolation. This name is constrained by the
+    ! current 8-character limit. It stands for AgSys ("AG") vegetation temperature
+    ! ("TVEG"), 24-hour average ("24").
+    call init_accum_field(name='AGTVEG24', units='K', &
+         desc='24hr average of vegetation temperature', accum_type='timeavg', accum_period=-1, &
+         subgrid_type='pft', numlev=1, init_value=0._r8)
+
+    ! NOTE(wjs, 2019-11-15) This name is constrained by the current 8-character limit. It
+    ! stands for AgSys ("AG") h2osoi_liq ("H2OS"), 24-hour average ("24").
+    call init_accum_field(name='AGH2OS24', units='mm H2O', &
+         desc='24hr average of h2osoi_liq', accum_type='timeavg', accum_period=-1, &
+         subgrid_type='column', numlev=nlevgrnd, type2d='levgrnd', init_value=0._r8)
+
+  end subroutine InitAccBuffer
+
+  !-----------------------------------------------------------------------
+  subroutine InitAccVars(this, bounds)
+    !
+    ! !DESCRIPTION:
+    ! Initialize accumulation variables for this instance
+    !
+    ! BUG(wjs, 2019-11-15, ESCOMP/ctsm#30) This needs to be called from outside a clump
+    ! loop (currently it operates on all elements of the arrays, rather than just
+    ! begp:endp or begc:endc).
+    !
+    ! !ARGUMENTS:
+    class(agsys_type), intent(inout) :: this
+    type(bounds_type), intent(in) :: bounds
+    !
+    ! !LOCAL VARIABLES:
+    integer  :: nstep
+
+    character(len=*), parameter :: subname = 'InitAccVars'
+    !-----------------------------------------------------------------------
+
+    nstep = get_nstep()
+
+    call extract_accum_field('AGTVEG24', this%t_veg24hr_patch, nstep)
+    call extract_accum_field('AGH2OS24', this%h2osoi_liq_24hr_col, nstep)
+
+  end subroutine InitAccVars
+
+  !-----------------------------------------------------------------------
+  subroutine UpdateAccVars(this, bounds)
+    !
+    ! !DESCRIPTION:
+    ! Update accumulation variables for this instance
+    !
+    ! BUG(wjs, 2019-11-15, ESCOMP/ctsm#30) This needs to be called from outside a clump
+    ! loop (currently it operates on all elements of the arrays, rather than just
+    ! begp:endp or begc:endc).
+    !
+    ! !ARGUMENTS:
+    class(agsys_type), intent(inout) :: this
+    type(bounds_type), intent(in) :: bounds
+    !
+    ! !LOCAL VARIABLES:
+    integer  :: nstep
+    integer  :: ier                   ! error status
+    real(r8), pointer :: rbufslp(:)
+    real(r8), pointer :: rbufmlc(:,:)
+
+    character(len=*), parameter :: subname = 'UpdateAccVars'
+    !-----------------------------------------------------------------------
+
+    nstep = get_nstep()
+
+    call update_accum_field('AGTVEG24', this%t_veg24hr_patch, nstep)
+    call extract_accum_field('AGTVEG24', this%t_veg24hr_patch, nstep)
+
+    call update_accum_field('AGH2OS24', this%h2osoi_liq_24hr_col, nstep)
+    call extract_accum_field('AGH2OS24', this%h2osoi_liq_24hr_col, nstep)
+
+  end subroutine UpdateAccVars
+
 end module AgSysType