Merge pull request #724 from rgknox/fates_runmean_vars

Add running mean functions

biogeochem/EDCanopyStructureMod.F90

@@ -42,7 +42,7 @@ module EDCanopyStructureMod
   use PRTGenericMod,          only : repro_organ
   use PRTGenericMod,          only : struct_organ
   use PRTGenericMod,          only : SetState
-
+  use FatesRunningMeanMod,    only : ema_lpa
 
   ! CIME Globals
   use shr_log_mod           , only : errMsg => shr_log_errMsg
@@ -669,6 +669,12 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr,bc_in)
                    call InitHydrCohort(currentSite,copyc)
                 endif
 
+                ! (keep as an example)
+                ! Initialize running means
+                !allocate(copyc%tveg_lpa)
+                !call copyc%tveg_lpa%InitRMean(ema_lpa, &
+                !     init_value=currentPatch%tveg_lpa%GetMean())
+
                 call copy_cohort(currentCohort, copyc)
 
                 newarea = currentCohort%c_area - cc_loss
@@ -1121,6 +1127,13 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr)
                    if( hlm_use_planthydro.eq.itrue ) then
                       call InitHydrCohort(CurrentSite,copyc)
                    endif
+
+                   ! (keep as an example)
+                   ! Initialize running means
+                   !allocate(copyc%tveg_lpa)
+                   !call copyc%tveg_lpa%InitRMean(ema_lpa,&
+                   !     init_value=currentPatch%tveg_lpa%GetMean())
+
                    call copy_cohort(currentCohort, copyc) !makes an identical copy...
 
                    newarea = currentCohort%c_area - cc_gain !new area of existing cohort

biogeochem/EDCohortDynamicsMod.F90

@@ -17,6 +17,7 @@ module EDCohortDynamicsMod
   use FatesConstantsMod     , only : fates_unset_r8
   use FatesConstantsMod     , only : nearzero
   use FatesConstantsMod     , only : calloc_abs_error
+  use FatesRunningMeanMod       , only : ema_lpa
   use FatesInterfaceTypesMod     , only : hlm_days_per_year
   use FatesInterfaceTypesMod     , only : nleafage
   use SFParamsMod           , only : SF_val_CWD_frac
@@ -308,7 +309,13 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh,   &
     call InitPRTBoundaryConditions(new_cohort)
 
 
+    ! Allocate running mean functions
 
+    !  (Keeping as an example)
+    !! allocate(new_cohort%tveg_lpa)
+    !! call new_cohort%tveg_lpa%InitRMean(ema_lpa,init_value=patchptr%tveg_lpa%GetMean())
+
+
     ! Recuits do not have mortality rates, nor have they moved any
     ! carbon when they are created.  They will bias our statistics
     ! until they have experienced a full day.  We need a newly recruited flag.
@@ -999,6 +1006,10 @@ subroutine DeallocateCohort(currentCohort)
 
      type(ed_cohort_type),intent(inout) :: currentCohort
 
+     !  (Keeping as an example)
+     ! Remove the running mean structure
+     ! deallocate(currentCohort%tveg_lpa)
+
      ! At this point, nothing should be pointing to current Cohort
      if (hlm_use_planthydro.eq.itrue) call DeallocateHydrCohort(currentCohort)
 
@@ -1162,6 +1173,11 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in)
                                       end do
                                    end if
 
+                                   !  (Keeping as an example)
+                                   ! Running mean fuses based on number density fraction just
+                                   ! like other variables
+                                   !!call currentCohort%tveg_lpa%FuseRMean(nextc%tveg_lpa,currentCohort%n/newn)
+
                                    ! new cohort age is weighted mean of two cohorts
                                    currentCohort%coage = &
                                         (currentCohort%coage * (currentCohort%n/(currentCohort%n + nextc%n))) + &
@@ -1794,6 +1810,7 @@ subroutine copy_cohort( currentCohort,copyc )
     n%size_by_pft_class = o%size_by_pft_class
     n%coage_class     = o%coage_class
     n%coage_by_pft_class = o%coage_by_pft_class
+
     ! This transfers the PRT objects over.
     call n%prt%CopyPRTVartypes(o%prt)
 
@@ -1803,6 +1820,10 @@ subroutine copy_cohort( currentCohort,copyc )
     n%tpu25top   = o%tpu25top
     n%kp25top    = o%kp25top
 
+    !  (Keeping as an example)
+    ! Copy over running means
+    ! call n%tveg_lpa%CopyFromDonor(o%tveg_lpa)
+
     ! CARBON FLUXES
     n%gpp_acc_hold    = o%gpp_acc_hold
     n%gpp_acc         = o%gpp_acc

biogeochem/EDMortalityFunctionsMod.F90

@@ -61,7 +61,6 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor
     real(r8),intent(out) :: smort  ! size dependent senescence term
     real(r8),intent(out) :: asmort ! age dependent senescence term 
 
-    integer  :: ifp
     real(r8) :: frac  ! relativised stored carbohydrate
     real(r8) :: leaf_c_target      ! target leaf biomass kgC
     real(r8) :: store_c
@@ -172,8 +171,8 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor
     !           Eastern US carbon sink.  Glob. Change Biol., 12, 2370-2390,              
     !           doi: 10.1111/j.1365-2486.2006.01254.x                                    
 
-    ifp = cohort_in%patchptr%patchno
-    temp_in_C = bc_in%t_veg24_pa(ifp) - tfrz
+    temp_in_C = cohort_in%patchptr%tveg24%GetMean() - tfrz
+
     temp_dep_fraction  = max(0.0_r8, min(1.0_r8, 1.0_r8 - (temp_in_C - &
                          EDPftvarcon_inst%freezetol(cohort_in%pft))/frost_mort_buffer) )
     frmort    = EDPftvarcon_inst%mort_scalar_coldstress(cohort_in%pft) * temp_dep_fraction

biogeochem/EDPatchDynamicsMod.F90

@@ -45,12 +45,14 @@ module EDPatchDynamicsMod
   use FatesInterfaceTypesMod    , only : bc_in_type
   use FatesInterfaceTypesMod    , only : hlm_days_per_year
   use FatesInterfaceTypesMod    , only : numpft
+  use FatesInterfaceTypesMod    , only : hlm_stepsize
   use FatesInterfaceTypesMod    , only : hlm_use_sp
   use FatesInterfaceTypesMod    , only : hlm_use_nocomp
   use FatesInterfaceTypesMod    , only : hlm_use_fixed_biogeog
   use FatesGlobals         , only : endrun => fates_endrun
   use FatesConstantsMod    , only : r8 => fates_r8
   use FatesConstantsMod    , only : itrue, ifalse
+  use FatesConstantsMod    , only : t_water_freeze_k_1atm
   use FatesPlantHydraulicsMod, only : InitHydrCohort
   use FatesPlantHydraulicsMod, only : AccumulateMortalityWaterStorage
   use FatesPlantHydraulicsMod, only : DeallocateHydrCohort
@@ -87,7 +89,8 @@ module EDPatchDynamicsMod
   use SFParamsMod,            only : SF_VAL_CWD_FRAC
   use EDParamsMod,            only : logging_event_code
   use EDParamsMod,            only : logging_export_frac
-
+  use FatesRunningMeanMod,    only : ema_24hr, fixed_24hr, ema_lpa
+
   ! CIME globals
   use shr_infnan_mod       , only : nan => shr_infnan_nan, assignment(=)
   use shr_log_mod          , only : errMsg => shr_log_errMsg
@@ -565,7 +568,7 @@ subroutine spawn_patches( currentSite, bc_in)
           allocate(new_patch_primary)
 
           call create_patch(currentSite, new_patch_primary, age, &
-                site_areadis_primary, primaryforest,fates_unset_int)
+                site_areadis_primary, primaryforest, fates_unset_int)
 
           ! Initialize the litter pools to zero, these
           ! pools will be populated by looping over the existing patches
@@ -583,7 +586,6 @@ subroutine spawn_patches( currentSite, bc_in)
 
        endif
 
-
        ! next create patch to receive secondary forest area
        if ( site_areadis_secondary .gt. nearzero) then
           allocate(new_patch_secondary)
@@ -674,6 +676,14 @@ subroutine spawn_patches( currentSite, bc_in)
                      new_patch, patch_site_areadis,bc_in)
              endif
 
+
+             ! Copy any means or timers from the original patch to the new patch
+             ! These values will inherit all info from the original patch
+             ! --------------------------------------------------------------------------
+             call new_patch%tveg24%CopyFromDonor(currentPatch%tveg24)
+             call new_patch%tveg_lpa%CopyFromDonor(currentPatch%tveg_lpa)
+
+
              ! --------------------------------------------------------------------------
              ! The newly formed patch from disturbance (new_patch), has now been given 
              ! some litter from dead plants and pre-existing litter from the donor patches.
@@ -693,7 +703,12 @@ subroutine spawn_patches( currentSite, bc_in)
                  nc%prt => null()
                  call InitPRTObject(nc%prt)
                  call InitPRTBoundaryConditions(nc)
-
+
+                 !  (Keeping as an example)
+                 ! Allocate running mean functions
+                 !allocate(nc%tveg_lpa)
+                 !call nc%tveg_lpa%InitRMean(ema_lpa,init_value=new_patch%tveg_lpa%GetMean())
+
                  call zero_cohort(nc)
 
                  ! nc is the new cohort that goes in the disturbed patch (new_patch)... currentCohort
@@ -1988,6 +2003,8 @@ end subroutine mortality_litter_fluxes
 
   subroutine create_patch(currentSite, new_patch, age, areap, label,nocomp_pft)
 
+    use FatesInterfaceTypesMod, only : hlm_current_tod,hlm_current_date,hlm_reference_date
+
     !
     ! !DESCRIPTION:
     !  Set default values for creating a new patch
@@ -2001,6 +2018,12 @@ subroutine create_patch(currentSite, new_patch, age, areap, label,nocomp_pft)
     real(r8), intent(in) :: areap                ! initial area of this patch in m2. 
     integer, intent(in)  :: label                ! anthropogenic disturbance label
     integer, intent(in)  :: nocomp_pft
+
+
+    ! Until bc's are pointed to by sites give veg a default temp [K]
+    real(r8), parameter :: temp_init_veg = 15._r8+t_water_freeze_k_1atm 
+
+
     ! !LOCAL VARIABLES:
     !---------------------------------------------------------------------
     integer :: el                                ! element loop index
@@ -2016,7 +2039,11 @@ subroutine create_patch(currentSite, new_patch, age, areap, label,nocomp_pft)
     allocate(new_patch%sabs_dif(hlm_numSWb))
     allocate(new_patch%fragmentation_scaler(currentSite%nlevsoil))
 
-
+    allocate(new_patch%tveg24)
+    call new_patch%tveg24%InitRMean(fixed_24hr,init_value=temp_init_veg,init_offset=real(hlm_current_tod,r8) )
+    allocate(new_patch%tveg_lpa)
+    call new_patch%tveg_lpa%InitRmean(ema_lpa,init_value=temp_init_veg)
+
     ! Litter
     ! Allocate, Zero Fluxes, and Initialize to "unset" values
 
@@ -2532,6 +2559,10 @@ subroutine fuse_2_patches(csite, dp, rp)
        write(fates_log(),*) 'trying to fuse patches with different anthro_disturbance_label values'
        call endrun(msg=errMsg(sourcefile, __LINE__))
     endif
+
+    ! Weighted mean of the running means
+    call rp%tveg24%FuseRMean(dp%tveg24,rp%area*inv_sum_area)
+    call rp%tveg_lpa%FuseRMean(dp%tveg_lpa,rp%area*inv_sum_area)
 
     rp%fuel_eff_moist       = (dp%fuel_eff_moist*dp%area + rp%fuel_eff_moist*rp%area) * inv_sum_area
     rp%livegrass            = (dp%livegrass*dp%area + rp%livegrass*rp%area) * inv_sum_area
@@ -2870,9 +2901,13 @@ subroutine dealloc_patch(cpatch)
        deallocate(cpatch%sabs_dir)
        deallocate(cpatch%sabs_dif)
        deallocate(cpatch%fragmentation_scaler)
-
     end if
 
+
+    ! Deallocate any running means
+    deallocate(cpatch%tveg24)
+    deallocate(cpatch%tveg_lpa)
+
     return
   end subroutine dealloc_patch
 

biogeochem/EDPhysiologyMod.F90

@@ -761,7 +761,7 @@ subroutine phenology( currentSite, bc_in )
     temp_in_C = 0._r8
     cpatch => CurrentSite%oldest_patch
     do while(associated(cpatch))
-       temp_in_C = temp_in_C + bc_in%t_veg24_pa(cpatch%patchno)*cpatch%area
+       temp_in_C = temp_in_C + cpatch%tveg24%GetMean()*cpatch%area
        cpatch => cpatch%younger
     end do
     temp_in_C = temp_in_C * area_inv - tfrz
@@ -2469,7 +2469,6 @@ subroutine fragmentation_scaler( currentPatch, bc_in)
     logical  :: use_century_tfunc = .false.
     logical  :: use_hlm_soil_scalar = .true. ! Use hlm input decomp fraction scalars
     integer  :: j
-    integer  :: ifp                          ! Index of a FATES Patch "ifp"
     real(r8) :: t_scalar                     ! temperature scalar
     real(r8) :: w_scalar                     ! moisture scalar
     real(r8) :: catanf                       ! hyperbolic temperature function from CENTURY
@@ -2480,7 +2479,6 @@ subroutine fragmentation_scaler( currentPatch, bc_in)
     catanf(t1) = 11.75_r8 +(29.7_r8 / pi) * atan( pi * 0.031_r8  * ( t1 - 15.4_r8 ))
     catanf_30 = catanf(30._r8)
 
-    ifp = currentPatch%patchno
     if(currentPatch%nocomp_pft_label.ne.0)then
 
        ! Use the hlm temp and moisture decomp fractions by default
@@ -2492,19 +2490,19 @@ subroutine fragmentation_scaler( currentPatch, bc_in)
        else
 
          if ( .not. use_century_tfunc ) then
-         !calculate rate constant scalar for soil temperature,assuming that the base rate constants
-         !are assigned for non-moisture limiting conditions at 25C.
-            if (bc_in%t_veg24_pa(ifp)  >=  tfrz) then
-              t_scalar = q10_mr**((bc_in%t_veg24_pa(ifp)-(tfrz+25._r8))/10._r8)
-                  !  Q10**((t_soisno(c,j)-(tfrz+25._r8))/10._r8)
+            !calculate rate constant scalar for soil temperature,assuming that the base rate constants
+            !are assigned for non-moisture limiting conditions at 25C.
+            if (currentPatch%tveg24%GetMean()  >=  tfrz) then
+               t_scalar = q10_mr**((currentPatch%tveg24%GetMean()-(tfrz+25._r8))/10._r8)
+               !  Q10**((t_soisno(c,j)-(tfrz+25._r8))/10._r8)
             else
-              t_scalar = (q10_mr**(-25._r8/10._r8))*(q10_froz**((bc_in%t_veg24_pa(ifp)-tfrz)/10._r8))
-                     !Q10**(-25._r8/10._r8))*(froz_q10**((t_soisno(c,j)-tfrz)/10._r8)
+               t_scalar = (q10_mr**(-25._r8/10._r8))*(q10_froz**((currentPatch%tveg24%GetMean()-tfrz)/10._r8))
+               !  Q10**(-25._r8/10._r8))*(froz_q10**((t_soisno(c,j)-tfrz)/10._r8)
             endif
          else
             ! original century uses an arctangent function to calculate the
             ! temperature dependence of decomposition
-            t_scalar = max(catanf(bc_in%t_veg24_pa(ifp)-tfrz)/catanf_30,0.01_r8)
+            t_scalar = max(catanf(currentPatch%tveg24%GetMean()-tfrz)/catanf_30,0.01_r8)
          endif
 
          !Moisture Limitations

fire/SFMainMod.F90

@@ -143,7 +143,7 @@ subroutine  fire_danger_index ( currentSite, bc_in)
 
     iofp = currentSite%oldest_patch%patchno
 
-    temp_in_C  = bc_in%t_veg24_pa(iofp) - tfrz
+    temp_in_C  = currentSite%oldest_patch%tveg24%GetMean() - tfrz
     rainfall   = bc_in%precip24_pa(iofp)*sec_per_day
     rh         = bc_in%relhumid24_pa(iofp)
 

main/EDInitMod.F90

@@ -74,7 +74,8 @@ module EDInitMod
   use PRTGenericMod,          only : nitrogen_element
   use PRTGenericMod,          only : phosphorus_element
   use PRTGenericMod,          only : SetState
-
+  use FatesSizeAgeTypeIndicesMod,only : get_age_class_index
+
   ! CIME GLOBALS
   use shr_log_mod               , only : errMsg => shr_log_errMsg
 
@@ -140,7 +141,9 @@ subroutine init_site_vars( site_in, bc_in, bc_out )
     endif
 
     allocate(site_in%use_this_pft(1:numpft))
+    allocate(site_in%area_by_age(1:nlevage))
 
+
     ! SP mode
     allocate(site_in%sp_tlai(1:numpft))
     allocate(site_in%sp_tsai(1:numpft))
@@ -251,6 +254,8 @@ subroutine zero_site( site_in )
 
     site_in%area_pft(:) = 0._r8
     site_in%use_this_pft(:) = fates_unset_int
+    site_in%area_by_age(:) = 0._r8
+
   end subroutine zero_site
 
   ! ============================================================================
@@ -431,7 +436,8 @@ subroutine init_patches( nsites, sites, bc_in)
     integer  :: s
     integer  :: el
     real(r8) :: age !notional age of this patch
-
+    integer  :: ageclass
+
     ! dummy locals
     real(r8) :: biomass_stock
     real(r8) :: litter_stock
@@ -897,12 +903,13 @@ subroutine init_cohorts( site_in, patch_in, bc_in)
        endif !use_this_pft
     enddo !numpft
 
-    ! Zero the mass flux pools of the new cohorts
-    !    temp_cohort => patch_in%tallest
-    !    do while(associated(temp_cohort))
-    !       call temp_cohort%prt%ZeroRates()
-    !       temp_cohort => temp_cohort%shorter
-    !    end do
+    ! (Keeping as an example)
+    ! Pass patch level temperature to the new cohorts (this is a nominal 15C right now)
+    !temp_cohort => patch_in%tallest
+    !do while(associated(temp_cohort))
+        !call temp_cohort%tveg_lpa%UpdateRmean(patch_in%tveg_lpa%GetMean())
+        !temp_cohort => temp_cohort%shorter
+    !end do
 
     call fuse_cohorts(site_in, patch_in,bc_in)
     call sort_cohorts(patch_in)

main/EDMainMod.F90

@@ -683,6 +683,8 @@ subroutine ed_update_site( currentSite, bc_in, bc_out )
 
     call TotalBalanceCheck(currentSite,final_check_id)
 
+    currentSite%area_by_age(:) = 0._r8
+
     currentPatch => currentSite%oldest_patch
     do while(associated(currentPatch))
 
@@ -694,8 +696,12 @@ subroutine ed_update_site( currentSite, bc_in, bc_out )
         ! This cohort count is used in the photosynthesis loop
         call count_cohorts(currentPatch)
 
-
+        ! Update the total area of by patch age class array 
+        currentSite%area_by_age(currentPatch%age_class) = &
+             currentSite%area_by_age(currentPatch%age_class) + currentPatch%area
+
         currentPatch => currentPatch%younger
+
     enddo
 
     ! The HLMs need to know about nutrient demand, and/or