From c7d7af6297fe541b07e34bdd4c446150acd14ff5 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Fri, 15 Nov 2019 10:45:48 -0700 Subject: [PATCH] Add 24-hour accumulators needed by AgSys Ideally these would be end-of-day averages rather than 24-hour averages (which can be different depending on when the run starts). But we don't have the functionality for end-of-day averages at this point (see ESCOMP/ctsm#839), so for now we're using 24-hour averages. --- src/agsys/ctsm_interface/AgSysInterface.F90 | 71 +++++++++- src/agsys/ctsm_interface/AgSysType.F90 | 121 ++++++++++++++++-- .../science/AgSysEnvironmentalInputs.F90 | 6 +- src/biogeophys/TemperatureType.F90 | 4 +- src/main/clm_driver.F90 | 4 + src/main/clm_initializeMod.F90 | 6 +- src/main/clm_instMod.F90 | 4 + 7 files changed, 201 insertions(+), 15 deletions(-) diff --git a/src/agsys/ctsm_interface/AgSysInterface.F90 b/src/agsys/ctsm_interface/AgSysInterface.F90 index 409f0982c8..454bb39f00 100644 --- a/src/agsys/ctsm_interface/AgSysInterface.F90 +++ b/src/agsys/ctsm_interface/AgSysInterface.F90 @@ -55,6 +55,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 = & @@ -94,14 +97,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 do fp = 1, num_pcropp p = filter_pcropp(fp) @@ -116,7 +124,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)) call DoTimeStep_Phenology_Placeholder( & @@ -170,4 +178,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 diff --git a/src/agsys/ctsm_interface/AgSysType.F90 b/src/agsys/ctsm_interface/AgSysType.F90 index ae5e080ee8..dfe1d23635 100644 --- a/src/agsys/ctsm_interface/AgSysType.F90 +++ b/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 @@ -82,7 +84,12 @@ module AgSysType 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(:) @@ -95,6 +102,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 = & @@ -163,7 +174,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 @@ -262,14 +274,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 diff --git a/src/agsys/science/AgSysEnvironmentalInputs.F90 b/src/agsys/science/AgSysEnvironmentalInputs.F90 index 42d76749a5..7dbcce0091 100644 --- a/src/agsys/science/AgSysEnvironmentalInputs.F90 +++ b/src/agsys/science/AgSysEnvironmentalInputs.F90 @@ -24,6 +24,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 :: SetValues ! Set values for the current point @@ -50,6 +53,7 @@ subroutine Init(this, nlevsoi) character(len=*), parameter :: subname = 'Init' !----------------------------------------------------------------------- + this%nlevsoi = nlevsoi allocate(this%h2osoi_liq_24hr(nlevsoi)) end subroutine Init @@ -77,7 +81,7 @@ subroutine SetValues(this, photoperiod, tair_max, tair_min, tc_24hr, h2osoi_liq_ 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 SetValues diff --git a/src/biogeophys/TemperatureType.F90 b/src/biogeophys/TemperatureType.F90 index 5c69733f8a..c3af098f1a 100644 --- a/src/biogeophys/TemperatureType.F90 +++ b/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', & diff --git a/src/main/clm_driver.F90 b/src/main/clm_driver.F90 index 501dffc9ec..6f210b0deb 100644 --- a/src/main/clm_driver.F90 +++ b/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 diff --git a/src/main/clm_initializeMod.F90 b/src/main/clm_initializeMod.F90 index b7549eb282..3148926824 100644 --- a/src/main/clm_initializeMod.F90 +++ b/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 !------------------------------------------------------------ diff --git a/src/main/clm_instMod.F90 b/src/main/clm_instMod.F90 index f4c11a5c11..4760bcdce7 100644 --- a/src/main/clm_instMod.F90 +++ b/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)