Merge branch 'agsys_accumulators' into agsys

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/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 = &
@@ -172,7 +183,8 @@ subroutine InitAllocate(this, bounds)
 
     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
 
@@ -285,14 +297,107 @@ subroutine InitCold(this, bounds, patch)
     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

src/agsys/science/AgSysEnvironmentalInputs.F90

@@ -25,6 +25,9 @@ module AgSysEnvironmentalInputs
      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]
 
+     ! Private data members
+     integer :: nlevsoi
+
    contains
      procedure, public :: Init       ! Allocate space for this instance (but don't set any values)
      procedure, public :: SetSpatiallyConstantValues ! Set spatially-constant values for this time
@@ -53,6 +56,7 @@ subroutine Init(this, nlevsoi)
     character(len=*), parameter :: subname = 'Init'
     !-----------------------------------------------------------------------
 
+    this%nlevsoi = nlevsoi
     allocate(this%h2osoi_liq_24hr(nlevsoi))
 
   end subroutine Init
@@ -99,7 +103,7 @@ subroutine SetSpatiallyVaryingValues(this, photoperiod, tair_max, tair_min, tc_2
     this%tair_max = tair_max
     this%tair_min = tair_min
     this%tc_24hr = tc_24hr
-    this%h2osoi_liq_24hr(:) = h2osoi_liq_24hr(:)
+    this%h2osoi_liq_24hr(1:this%nlevsoi) = h2osoi_liq_24hr(1:this%nlevsoi)
 
   end subroutine SetSpatiallyVaryingValues
 

src/biogeophys/TemperatureType.F90

@@ -1139,12 +1139,12 @@ subroutine InitAccBuffer (this, bounds)
 
     this%t_veg24_patch(bounds%begp:bounds%endp) = spval
     call init_accum_field (name='T_VEG24', units='K',                                              &
-         desc='24hr average of vegetation temperature',  accum_type='runmean', accum_period=-1,    &
+         desc='24hr running mean of vegetation temperature',  accum_type='runmean', accum_period=-1,    &
          subgrid_type='pft', numlev=1, init_value=0._r8)
 
     this%t_veg240_patch(bounds%begp:bounds%endp) = spval
     call init_accum_field (name='T_VEG240', units='K',                                             &
-         desc='240hr average of vegetation temperature',  accum_type='runmean', accum_period=-10,  &
+         desc='240hr running mean of vegetation temperature',  accum_type='runmean', accum_period=-10,  &
          subgrid_type='pft', numlev=1, init_value=0._r8)
 
     call init_accum_field(name='TREFAV', units='K', &

src/main/clm_driver.F90

@@ -1244,6 +1244,10 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro
                temperature_inst%t_ref2m_patch, temperature_inst%t_soisno_col)
        end if
 
+       if (use_crop_agsys) then
+          call agsys_interface_inst%UpdateAccVars(bounds_proc)
+       end if
+
        call t_stopf('accum')
     end if
 

src/main/clm_initializeMod.F90

@@ -277,7 +277,7 @@ subroutine initialize2( )
     use clm_varcon            , only : spval
     use clm_varctl            , only : finidat, finidat_interp_source, finidat_interp_dest, fsurdat
     use clm_varctl            , only : use_century_decomp, single_column, scmlat, scmlon, use_cn, use_fates
-    use clm_varctl            , only : use_crop, ndep_from_cpl
+    use clm_varctl            , only : use_crop, use_crop_agsys, ndep_from_cpl
     use clm_varorb            , only : eccen, mvelpp, lambm0, obliqr
     use clm_time_manager      , only : get_step_size_real, get_curr_calday
     use clm_time_manager      , only : get_curr_date, get_nstep, advance_timestep 
@@ -668,6 +668,10 @@ subroutine initialize2( )
        call crop_inst%initAccVars(bounds_proc)
     end if
 
+    if (use_crop_agsys) then
+       call agsys_interface_inst%initAccVars(bounds_proc)
+    end if
+
     !------------------------------------------------------------
     ! Read monthly vegetation
     !------------------------------------------------------------

src/main/clm_instMod.F90

@@ -470,6 +470,10 @@ subroutine clm_instInit(bounds)
        call crop_inst%InitAccBuffer(bounds)
     end if
 
+    if (use_crop_agsys) then
+       call agsys_interface_inst%InitAccBuffer(bounds)
+    end if
+
     call print_accum_fields()
 
     call ncd_pio_closefile(params_ncid)