diff --git a/src/agsys/ctsm_interface/AgSysInterface.F90 b/src/agsys/ctsm_interface/AgSysInterface.F90 index 43b021161c..4b4fd046b8 100644 --- a/src/agsys/ctsm_interface/AgSysInterface.F90 +++ b/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 diff --git a/src/agsys/ctsm_interface/AgSysType.F90 b/src/agsys/ctsm_interface/AgSysType.F90 index eb37ba67b2..ae5e080ee8 100644 --- a/src/agsys/ctsm_interface/AgSysType.F90 +++ b/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 diff --git a/src/agsys/science/AgSysEnvironmentalInputs.F90 b/src/agsys/science/AgSysEnvironmentalInputs.F90 new file mode 100644 index 0000000000..42d76749a5 --- /dev/null +++ b/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 diff --git a/src/agsys/science/AgSysPhases.F90 b/src/agsys/science/AgSysPhases.F90 index 805664fd6c..4e97639e99 100644 --- a/src/agsys/science/AgSysPhases.F90 +++ b/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 diff --git a/src/agsys/science/AgSysPlaceholder.F90 b/src/agsys/science/AgSysPlaceholder.F90 index f06c1d6e69..ffd89ff3de 100644 --- a/src/agsys/science/AgSysPlaceholder.F90 +++ b/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 diff --git a/src/main/clm_instMod.F90 b/src/main/clm_instMod.F90 index 89e51651ce..f4c11a5c11 100644 --- a/src/main/clm_instMod.F90 +++ b/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) diff --git a/src/main/histFileMod.F90 b/src/main/histFileMod.F90 index 23690ad51f..377f79f55d 100644 --- a/src/main/histFileMod.F90 +++ b/src/main/histFileMod.F90 @@ -12,7 +12,7 @@ module histFileMod use shr_sys_mod , only : shr_sys_flush use spmdMod , only : masterproc use abortutils , only : endrun - use clm_varctl , only : iulog, use_vertsoilc, use_fates + use clm_varctl , only : iulog, use_vertsoilc, use_fates, use_crop_agsys use clm_varcon , only : spval, ispval, dzsoi_decomp use clm_varcon , only : grlnd, nameg, namel, namec, namep, nameCohort use decompMod , only : get_proc_bounds, get_proc_global, bounds_type @@ -188,7 +188,7 @@ module histFileMod character(len=max_chars) :: units ! units character(len=hist_dim_name_length) :: type1d ! pointer to first dimension type from data type (nameg, etc) character(len=hist_dim_name_length) :: type1d_out ! hbuf first dimension type from data type (nameg, etc) - character(len=hist_dim_name_length) :: type2d ! hbuf second dimension type ["levgrnd","levlak","numrad","ltype","natpft","cft","glc_nec","elevclas","subname(n)"] + character(len=hist_dim_name_length) :: type2d ! hbuf second dimension type ["levgrnd","levlak","numrad","ltype","natpft","cft","glc_nec","elevclas","agsys_phases","subname(n)"] integer :: beg1d ! on-node 1d clm pointer start index integer :: end1d ! on-node 1d clm pointer end index integer :: num1d ! size of clm pointer first dimension (all nodes) @@ -1057,7 +1057,7 @@ subroutine hist_update_hbuf(bounds) integer :: f ! field index integer :: num2d ! size of second dimension (e.g. number of vertical levels) character(len=*),parameter :: subname = 'hist_update_hbuf' - character(len=hist_dim_name_length) :: type2d ! hbuf second dimension type ["levgrnd","levlak","numrad","ltype","natpft","cft","glc_nec","elevclas","subname(n)"] + character(len=hist_dim_name_length) :: type2d ! hbuf second dimension type ["levgrnd","levlak","numrad","ltype","natpft","cft","glc_nec","elevclas","agsys_phases","subname(n)"] !----------------------------------------------------------------------- do t = 1,ntapes @@ -1887,6 +1887,7 @@ subroutine htape_create (t, histrest) use clm_varpar , only : nlevgrnd, nlevsno, nlevlak, nlevurb, numrad, nlevcan, nvegwcs,nlevsoi use clm_varpar , only : natpft_size, cft_size, maxpatch_glcmec, nlevdecomp_full use landunit_varcon , only : max_lunit + use AgSysRuntimeConstants, only : agsys_max_phases use clm_varctl , only : caseid, ctitle, fsurdat, finidat, paramfile use clm_varctl , only : version, hostname, username, conventions, source use domainMod , only : ldomain @@ -2069,6 +2070,10 @@ subroutine htape_create (t, histrest) call ncd_defdim(lnfid, 'fates_levcnlfpf', nlevleaf * nclmax * maxveg_fates, dimid) end if + if (use_crop_agsys) then + call ncd_defdim(lnfid, 'agsys_phases', agsys_max_phases, dimid) + end if + if ( .not. lhistrest )then call ncd_defdim(lnfid, 'hist_interval', 2, hist_interval_dimid) call ncd_defdim(lnfid, 'time', ncd_unlimited, time_dimid) @@ -4764,6 +4769,7 @@ subroutine hist_addfld2d (fname, type2d, units, avgflag, long_name, type1d_out, use clm_varpar , only : nlevgrnd, nlevsno, nlevlak, numrad, nlevdecomp_full, nlevcan, nvegwcs,nlevsoi use clm_varpar , only : natpft_size, cft_size, maxpatch_glcmec use landunit_varcon , only : max_lunit + use AgSysRuntimeConstants, only : agsys_max_phases ! ! !ARGUMENTS: character(len=*), intent(in) :: fname ! field name @@ -4883,13 +4889,15 @@ subroutine hist_addfld2d (fname, type2d, units, avgflag, long_name, type1d_out, case ('levsno') num2d = nlevsno case ('nlevcan') - num2d = nlevcan + num2d = nlevcan case ('nvegwcs') - num2d = nvegwcs + num2d = nvegwcs + case ('agsys_phases') + num2d = agsys_max_phases case default write(iulog,*) trim(subname),' ERROR: unsupported 2d type ',type2d, & ' currently supported types for multi level fields are: ', & - '[levgrnd,levsoi,levlak,numrad,levdcmp,levtrc,ltype,natpft,cft,glc_nec,elevclas,levsno,nvegwcs]' + '[levgrnd,levsoi,levlak,numrad,levdcmp,levtrc,ltype,natpft,cft,glc_nec,elevclas,levsno,nvegwcs,agsys_phases]' call endrun(msg=errMsg(sourcefile, __LINE__)) end select