Merge pull request ESCOMP#14 from billsacks/agsys

Fix accumulators and add design document

doc/design/agsys.rst

@@ -0,0 +1,56 @@
+.. sectnum::
+
+.. contents::
+
+============================
+ Overview and general notes
+============================
+
+This document describes some key design aspects of the code under src/agsys.
+
+An important driver of the AgSys design was the desire to have the core science code be
+model-independent, so that the same code could be used in CTSM, Noah-MP, etc. Code in the
+``science`` directory does not depend at all on CTSM; it only depends on other code in
+``science`` and code in the ``ctsm_wrappers`` directory. The wrappers would be replaced by
+different wrappers if you plugged this into a different model.
+
+Partly to avoid dependence on CTSM-specific data structures and concepts, code in
+``science`` operates on single points. Normally we try to avoid single-point routines for
+performance reasons, but in this case, much of what's done just happens once per day, not
+every time step, so this seemed less important. And given the heavily-conditional nature
+of the code, it wouldn't be vectorizable even if it were written as loops.
+
+===================================================================
+ Fundamental data structures: single-point or structures of arrays
+===================================================================
+
+We considered whether to have the fundamental data structures be single-point data
+structures or structures of arrays as is done elsewhere throughout CTSM. We went with the
+latter because we need these structures of arrays to use CTSM infrastructure for
+diagnostics (history), restart, accumulators, etc. We could have used arrays of structures
+as our fundamental data, which would have simplified the argument lists in the science
+code, but this would have required copy-outs of all data in the interface layer for the
+sake of diagnostics and restarts, which felt messier and more error-prone.
+
+Note, though, that we *do* copy various input arguments into a single-point structure
+before calling the AgSys routines, in order to simplify argument lists, and to allow a
+consistent interface between different crops' implementations of a polymorphic subroutine
+(without resorting to passing the superset of all inputs needed by any crop's
+implementation).
+
+There are arguments for and against storing AgSys's state variables in a similar
+single-ponit structure, and doing a copy-out (into CTSM's patch-level arrays) at the end
+of the AgSys calls for the sake of diagnostics and restart. An advantage is that it would
+simplify the interface to (and possibly within) the science code, and it would make the
+science code more usable outside the interface layer (since it would have its own state
+variable structures). It also feels more right to have the internal agsys code define and
+use its own state variables, rather than having these state variables defined and
+individually passed. However, we wouldn't really be able to hide the variables needed
+internally in AgSys, since we'd still need duplicates of all of those variables in the
+interface layer. And in some ways I like passing output variables one-by-one, so you can
+clearly see what's being set by each routine. Also, this would mean having each variable
+stored persistently in two places, which is a bit confusing, not to mention wasteful of
+memory. So, given the pros and cons in both directions, and the fact that we started with
+the implementation of having each inout variable passed individually (so it would take
+more work to change it), for now we're sticking with having each inout variable passed
+individually and stored fundamentally in a structure of arrays.

src/agsys/ctsm_interface/AgSysInterface.F90

@@ -23,6 +23,8 @@ module AgSysInterface
   use AgSysCropTypeGeneric, only : agsys_cultivars_of_crop_type
   use AgSysRoot, only : get_soil_condition
   use AgSysPhenology, only : AgSysDeterminePlanting, AgSysRunPhenology
+  use TemperatureType, only : temperature_type
+  use WaterStateBulkType, only : waterstatebulk_type
   !
   implicit none
   private
@@ -239,21 +241,23 @@ subroutine InitAccVars(this, bounds)
   end subroutine InitAccVars
 
   !-----------------------------------------------------------------------
-  subroutine UpdateAccVars(this, bounds)
+  subroutine UpdateAccVars(this, bounds, temperature_inst, waterstatebulk_inst)
     !
     ! !DESCRIPTION:
     ! Update accumulated variables
     !
     ! !ARGUMENTS:
     class(agsys_interface_type), intent(inout) :: this
     type(bounds_type), intent(in) :: bounds
+    type(temperature_type), intent(in) :: temperature_inst
+    type(waterstatebulk_type), intent(in) :: waterstatebulk_inst
     !
     ! !LOCAL VARIABLES:
 
     character(len=*), parameter :: subname = 'UpdateAccVars'
     !-----------------------------------------------------------------------
 
-    call this%agsys_inst%UpdateAccVars(bounds)
+    call this%agsys_inst%UpdateAccVars(bounds, temperature_inst, waterstatebulk_inst)
 
   end subroutine UpdateAccVars
 

src/agsys/ctsm_interface/AgSysType.F90

@@ -20,6 +20,8 @@ module AgSysType
   use AgSysConstants,           only : crop_type_not_handled, crop_type_maize
   use AgSysEnvironmentalInputs, only : agsys_environmental_inputs_type
   use AgSysRoot,                only : agsys_soil_property_type, agsys_soil_condition_type, agsys_root_type     
+  use TemperatureType, only : temperature_type
+  use WaterStateBulkType, only : waterstatebulk_type
   !
   implicit none
   private
@@ -368,7 +370,7 @@ subroutine InitAccVars(this, bounds)
   end subroutine InitAccVars
 
   !-----------------------------------------------------------------------
-  subroutine UpdateAccVars(this, bounds)
+  subroutine UpdateAccVars(this, bounds, temperature_inst, waterstatebulk_inst)
     !
     ! !DESCRIPTION:
     ! Update accumulation variables for this instance
@@ -380,22 +382,26 @@ subroutine UpdateAccVars(this, bounds)
     ! !ARGUMENTS:
     class(agsys_type), intent(inout) :: this
     type(bounds_type), intent(in) :: bounds
+    type(temperature_type), intent(in) :: temperature_inst
+    type(waterstatebulk_type), intent(in) :: waterstatebulk_inst
     !
     ! !LOCAL VARIABLES:
     integer  :: nstep
     integer  :: ier                   ! error status
-    real(r8), pointer :: rbufslp(:)
-    real(r8), pointer :: rbufmlc(:,:)
+    real(r8), pointer :: temp_ptr(:,:)
 
     character(len=*), parameter :: subname = 'UpdateAccVars'
     !-----------------------------------------------------------------------
 
     nstep = get_nstep()
 
-    call update_accum_field('AGTVEG24', this%t_veg24hr_patch, nstep)
+    call update_accum_field('AGTVEG24', temperature_inst%t_veg_patch, nstep)
     call extract_accum_field('AGTVEG24', this%t_veg24hr_patch, nstep)
 
-    call update_accum_field('AGH2OS24', this%h2osoi_liq_24hr_col, nstep)
+    ! Need to use a temporary pointer in order to only include levels from 1:nlevgrnd
+    ! (so, excluding snow levels)
+    temp_ptr(bounds%begc:, 1:) => waterstatebulk_inst%h2osoi_liq_col(:, 1:)
+    call update_accum_field('AGH2OS24', temp_ptr, nstep)
     call extract_accum_field('AGH2OS24', this%h2osoi_liq_24hr_col, nstep)
 
   end subroutine UpdateAccVars

src/main/clm_driver.F90

@@ -1245,7 +1245,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro
        end if
 
        if (use_crop_agsys) then
-          call agsys_interface_inst%UpdateAccVars(bounds_proc)
+          call agsys_interface_inst%UpdateAccVars(bounds_proc, temperature_inst, water_inst%waterstatebulk_inst)
        end if
 
        call t_stopf('accum')