diff --git a/core/biogeochem/POP.F90 b/core/biogeochem/POP.F90 index 5c62059cb..b0487ed34 100755 --- a/core/biogeochem/POP.F90 +++ b/core/biogeochem/POP.F90 @@ -82,60 +82,59 @@ MODULE POP_Constants ! REAL(dp),PARAMETER:: CrowdingFactor = 0.0128 ! REAL(dp),PARAMETER:: ALPHA_CPC = 3.0 - REAL(dp), PARAMETER :: FULTON_ALPHA = 3.5_dp ! recruitment scalar alpha in Fulton (1991) - REAL(dp), PARAMETER :: DENSINDIV_MAX = 0.2_dp ! 0.5 ! Maximum density of individuals within a cohort indiv/m2 - REAL(dp), PARAMETER :: DENSINDIV_MIN = 1.0e-9_dp ! - REAL(dp), PARAMETER :: Kbiometric = 50.0_dp ! Constant in height-diameter relationship - REAL(dp), PARAMETER :: WD = 300.0_dp ! Wood density kgC/m3 + REAL(dp) :: FULTON_ALPHA = 3.5_dp ! recruitment scalar alpha in Fulton (1991) + REAL(dp) :: DENSINDIV_MAX = 0.2_dp ! 0.5 ! Maximum density of individuals within a cohort indiv/m2 + REAL(dp) :: DENSINDIV_MIN = 1.0e-9_dp ! + REAL(dp) :: Kbiometric = 50.0_dp ! Constant in height-diameter relationship + REAL(dp) :: WD = 300.0_dp ! Wood density kgC/m3 ! threshold growth efficiency for enhanced mortality (higher value gives higher biomass turnover) - REAL(dp), PARAMETER :: GROWTH_EFFICIENCY_MIN = 0.009_dp ! 0.0095 ! 0.0089 ! 0.0084 - REAL(dp), PARAMETER :: Pmort = 5.0_dp ! exponent in mortality formula - REAL(dp), PARAMETER :: MORT_MAX = 0.3_dp ! upper asymptote for enhanced mortality - REAL(dp), PARAMETER :: THETA_recruit = 0.95_dp ! shape parameter in recruitment equation - REAL(dp), PARAMETER :: CMASS_STEM_INIT = 1.0e-4_dp ! initial biomass kgC/m2 - REAL(dp), PARAMETER :: POWERbiomass = 0.67_dp ! exponent for biomass in proportion to which cohorts preempt resources - REAL(dp), PARAMETER :: POWERGrowthEfficiency = 0.67_dp - REAL(dp), PARAMETER :: CrowdingFactor = 0.043_dp ! 0.043 ! 0.039 !0.029 ! 0.033 - REAL(dp), PARAMETER :: ALPHA_CPC = 3.5_dp - REAL(dp), PARAMETER :: k_allom1 = 200.0_dp ! crown area = k_allom1 * diam ** k_rp - REAL(dp), PARAMETER :: k_rp = 1.67_dp ! constant in crown area relation to tree diameter - REAL(dp), PARAMETER :: ksapwood = 0.05_dp ! rate constant for conversion of sapwood to heartwood (y-1) - REAL(dp), PARAMETER :: Q=7.0_dp ! governs rate of increase of mortality with age (2=exponential) - REAL, PARAMETER :: rshootfrac = 0.63 - REAL(dp), PARAMETER :: shootfrac = real(rshootfrac,dp) - REAL(dp), PARAMETER :: CtoNw = 400.0_dp - REAL(dp), PARAMETER :: CtoNl = 60.0_dp - REAL(dp), PARAMETER :: CtoNr = 70.0_dp - REAL(dp), PARAMETER :: N_EXTENT = 2.0_dp ! multiple of crown diameters within which tree competes with other cohorts + REAL(dp) :: GROWTH_EFFICIENCY_MIN = 0.009_dp ! 0.0095 ! 0.0089 ! 0.0084 + REAL(dp) :: Pmort = 5.0_dp ! exponent in mortality formula + REAL(dp) :: MORT_MAX = 0.3_dp ! upper asymptote for enhanced mortality + REAL(dp) :: THETA_recruit = 0.95_dp ! shape parameter in recruitment equation + REAL(dp) :: CMASS_STEM_INIT = 1.0e-4_dp ! initial biomass kgC/m2 + REAL(dp) :: POWERbiomass = 0.67_dp ! exponent for biomass in proportion to which cohorts preempt resources + REAL(dp) :: POWERGrowthEfficiency = 0.67_dp + REAL(dp) :: CrowdingFactor = 0.043_dp ! 0.043 ! 0.039 !0.029 ! 0.033 + REAL(dp) :: ALPHA_CPC = 3.5_dp + REAL(dp) :: k_allom1 = 200.0_dp ! crown area = k_allom1 * diam ** k_rp + REAL(dp) :: k_rp = 1.67_dp ! constant in crown area relation to tree diameter + REAL(dp) :: ksapwood = 0.05_dp ! rate constant for conversion of sapwood to heartwood (y-1) + REAL(dp) :: Q=7.0_dp ! governs rate of increase of mortality with age (2=exponential) + REAL(dp) :: shootfrac = 0.63_dp + REAL :: rshootfrac = 0.63 + REAL(dp) :: CtoNw = 400.0_dp + REAL(dp) :: CtoNl = 60.0_dp + REAL(dp) :: CtoNr = 70.0_dp REAL(dp), PARAMETER :: EPS = 1.0e-12_dp - INTEGER(i4b), PARAMETER :: NLAYER = 1 ! number of vertical veg layers (1 is currently the only option) - INTEGER(i4b), PARAMETER :: NCOHORT_MAX = 20 ! maximum number of cohorts - INTEGER(i4b), PARAMETER :: NDISTURB = 1 ! number of disturbance regimes (1 (total only) or 2 (partial and total)) - INTEGER(i4b), PARAMETER :: PATCH_REPS = 10 ! higher number reduces 'noise' - INTEGER(i4b), PARAMETER :: NAGE_MAX = 1 ! number of maxium ages - INTEGER(i4b), PARAMETER :: PATCH_REPS1 = 60 ! number of first dist years - INTEGER(i4b), PARAMETER :: PATCH_REPS2 = 1 ! number of second dist years - INTEGER(i4b), PARAMETER :: NPATCH = PATCH_REPS1*PATCH_REPS2 - INTEGER(i4b), PARAMETER :: NPATCH1D = NPATCH - INTEGER(i4b), PARAMETER :: NPATCH2D = NPATCH - INTEGER(i4b), PARAMETER :: HEIGHT_BINS = 12 ! number of height categories to keep track of for diagnostics - REAL(dp), PARAMETER :: BIN_POWER = 1.4_dp ! bins have muscles - ! Time base factor (to be multiplied by mean dist interval to give TIMEBASE) - ! for sampling disturbance probabilities from Poisson distribution - INTEGER(i4b), PARAMETER :: TIMEBASE_FACTOR=50 + INTEGER(i4b) :: NLayer = 1 + INTEGER(i4b) :: NCOHORT_MAX = 20 ! maximum number of cohorts + INTEGER(i4b) :: NDISTURB = 1 ! number of disturbance regimes (1 (total only) or 2 (partial and total)) + INTEGER(i4b) :: PATCH_REPS = 10 ! higher number reduces 'noise' + INTEGER(i4b) :: PATCH_REPS1 = 60 ! number of first dist years + INTEGER(i4b) :: PATCH_REPS2 = 1 ! number of second dist years + INTEGER(i4b) :: NPATCH + INTEGER(i4b) :: NPATCH2D + INTEGER(i4b) :: HEIGHT_BINS = 12 ! number of height categories to keep track of for diagnostics + REAL(dp) :: BIN_POWER = 1.4_dp ! bins have muscles + + ! This seems like it should be meaningful for sampling the distribution + !! Time base factor (to be multiplied by mean dist interval to give TIMEBASE) + !! for sampling disturbance probabilities from Poisson distribution + !INTEGER(i4b) :: TIMEBASE_FACTOR=50 REAL(dp), PARAMETER :: PI=3.14159265358979323846264_dp ! 0 == default; 1 = top-end allometry (requires precip as input to POPSTEP); 2 = Allometry following Williams 2005, Model 5b - INTEGER(i4b), PARAMETER :: ALLOM_SWITCH = 2 + INTEGER(i4b) :: ALLOM_SWITCH = 2 ! 0 == binnned max height variable; 1 = continuous (needs lots of memory); 2 = binned by integer heights - INTEGER(i4b), PARAMETER :: MAX_HEIGHT_SWITCH = 2 - INTEGER(i4b), PARAMETER :: RESOURCE_SWITCH = 1 ! 0 = default; 1 fraction net resource uptake - INTEGER(i4b), PARAMETER :: RECRUIT_SWITCH = 1 ! 0 = default, 1 = Pgap-dependence - INTEGER(i4b), PARAMETER :: INTERP_SWITCH = 1 ! 0 = sum over weighted patches, 1 = sum over interpolated patches - INTEGER(i4b), PARAMETER :: SMOOTH_SWITCH = 0 ! smooth disturbance flux - INTEGER(i4b), PARAMETER :: NYEAR_WINDOW = 5 ! one-side of smoothing window (y) - INTEGER(i4b), PARAMETER :: NYEAR_SMOOTH = 2*NYEAR_WINDOW + 1 ! smoothing window (y) - INTEGER(i4b), PARAMETER :: NYEAR_HISTORY = NYEAR_SMOOTH-NYEAR_WINDOW - INTEGER(i4b), PARAMETER :: AGEMAX = 1000 + INTEGER(i4b) :: MAX_HEIGHT_SWITCH = 2 + INTEGER(i4b) :: RESOURCE_SWITCH = 1 ! 0 = default; 1 fraction net resource uptake + INTEGER(i4b) :: RECRUIT_SWITCH = 1 ! 0 = default, 1 = Pgap-dependence + INTEGER(i4b) :: INTERP_SWITCH = 1 ! 0 = sum over weighted patches, 1 = sum over interpolated patches + INTEGER(i4b) :: SMOOTH_SWITCH = 0 ! smooth disturbance flux + INTEGER(i4b) :: NYEAR_WINDOW = 5 ! one-side of smoothing window (y) + INTEGER(i4b) :: NYEAR_SMOOTH ! smoothing window (y) + INTEGER(i4b) :: NYEAR_HISTORY + INTEGER(i4b) :: AGEMAX = 1000 END MODULE POP_Constants @@ -146,8 +145,8 @@ END MODULE POP_Constants MODULE POP_Types USE TYPEdef, ONLY: dp, i4b - USE POP_Constants, ONLY: NCOHORT_MAX, NLAYER, HEIGHT_BINS, NDISTURB, NPATCH, NPATCH2D, & - NYEAR_HISTORY, AGEMAX + USE POP_Constants, ONLY: NCOHORT_MAX, HEIGHT_BINS, NDISTURB, NPATCH, NPATCH2D, & + NYEAR_HISTORY, AGEMAX, NLayer IMPLICIT NONE @@ -175,17 +174,17 @@ MODULE POP_Types REAL(dp) :: Croot END TYPE Cohort - TYPE Layer - TYPE(Cohort), DIMENSION(NCOHORT_MAX) :: Cohort + TYPE layer + TYPE(Cohort), DIMENSION(:), ALLOCATABLE :: Cohort INTEGER(i4b) :: ncohort ! number of cohorts with density >0 REAL(dp) :: biomass ! layer biomass REAL(dp) :: density ! layer tree density REAL(dp) :: hmean ! layer mean tree height (weighted mean over patches) REAL(dp) :: hmax ! layer max tree height - END TYPE Layer + END TYPE layer TYPE Patch - TYPE(Layer), DIMENSION(NLAYER) :: Layer + TYPE(layer), DIMENSION(:), ALLOCATABLE :: Layer REAL(dp) :: factor_recruit REAL(dp) :: pgap REAL(dp) :: lai @@ -206,9 +205,9 @@ MODULE POP_Types REAL(dp) :: sapwood_area_loss REAL(dp) :: growth ! biomass growth in each patch due to stem increment REAL(dp) :: area_growth ! basal area growth in each patch due to stem increment - INTEGER(i4b) :: disturbance_interval(NDISTURB) ! prescribed disturbance(s) interval for this patch - INTEGER(i4b) :: first_disturbance_year(NDISTURB) - INTEGER(i4b) :: age(NDISTURB) ! number of years since last disturbance(s) + INTEGER(i4b), DIMENSION(:), ALLOCATABLE :: disturbance_interval ! prescribed disturbance(s) interval for this patch + INTEGER(i4b), DIMENSION(:), ALLOCATABLE :: first_disturbance_year + INTEGER(i4b), DIMENSION(:), ALLOCATABLE :: age ! number of years since last disturbance(s) INTEGER(i4b) :: id REAL(dp) :: frac_NPP REAL(dp) :: frac_respiration @@ -217,25 +216,26 @@ MODULE POP_Types END TYPE Patch TYPE Landscape - TYPE(Patch), DIMENSION(NPATCH2D) :: patch - REAL(dp), DIMENSION(NPATCH2D) :: freq ! patch weighting - REAL(dp), DIMENSION(NPATCH2D) :: freq_old ! patch weighting (previous time-step) - REAL(dp), DIMENSION(NPATCH2D) :: fire_freq ! - REAL(dp), DIMENSION(NPATCH2D) :: fire_freq_old ! - REAL(dp), DIMENSION(NPATCH2D) :: cat_freq ! - REAL(dp), DIMENSION(NPATCH2D) :: cat_freq_old ! - REAL(dp), DIMENSION(NPATCH2D,NDISTURB) :: freq_ranked_age_unique ! unique age weighting - INTEGER(i4b), DIMENSION(NPATCH2D, NDISTURB) :: ranked_age_unique ! unique age - INTEGER(i4b), DIMENSION(NDISTURB) :: n_age ! number of unique ages - REAL(dp), DIMENSION(NLAYER) :: biomass ! landscape stem biomass (weighted mean over patches) - REAL(dp), DIMENSION(NLAYER) :: density ! landscape tree density (weighted mean over patches) - REAL(dp), DIMENSION(NLAYER) :: hmean ! landscape mean treen height (weighted mean over patches) - REAL(dp), DIMENSION(NLAYER) :: hmax ! landscape max tree height - REAL(dp), DIMENSION(HEIGHT_BINS) :: cmass_stem_bin ! biomass by height bin - REAL(dp), DIMENSION(HEIGHT_BINS) :: densindiv_bin ! density by height bin - REAL(dp), DIMENSION(HEIGHT_BINS) :: height_bin ! mean height in each bin - REAL(dp), DIMENSION(HEIGHT_BINS) :: diameter_bin ! mean diameter in each bin - CHARACTER(100), DIMENSION(HEIGHT_BINS) :: bin_labels ! text strings for bin bounds + ! Dangerous having type and variable with same name + TYPE(Patch), DIMENSION(:), ALLOCATABLE :: patch + REAL(dp), DIMENSION(:), ALLOCATABLE :: freq ! patch weighting + REAL(dp), DIMENSION(:), ALLOCATABLE :: freq_old ! patch weighting (previous time-step) + REAL(dp), DIMENSION(:), ALLOCATABLE :: fire_freq ! + REAL(dp), DIMENSION(:), ALLOCATABLE :: fire_freq_old ! + REAL(dp), DIMENSION(:), ALLOCATABLE :: cat_freq ! + REAL(dp), DIMENSION(:), ALLOCATABLE :: cat_freq_old ! + REAL(dp), DIMENSION(:,:), ALLOCATABLE :: freq_ranked_age_unique ! unique age weighting + INTEGER(i4b), DIMENSION(:,:), ALLOCATABLE :: ranked_age_unique ! unique age + INTEGER(i4b), DIMENSION(:), ALLOCATABLE :: n_age ! number of unique ages + REAL(dp) :: biomass ! landscape stem biomass (weighted mean over patches) + REAL(dp) :: density ! landscape tree density (weighted mean over patches) + REAL(dp) :: hmean ! landscape mean treen height (weighted mean over patches) + REAL(dp) :: hmax ! landscape max tree height + REAL(dp), DIMENSION(:), ALLOCATABLE :: cmass_stem_bin ! biomass by height bin + REAL(dp), DIMENSION(:), ALLOCATABLE :: densindiv_bin ! density by height bin + REAL(dp), DIMENSION(:), ALLOCATABLE :: height_bin ! mean height in each bin + REAL(dp), DIMENSION(:), ALLOCATABLE :: diameter_bin ! mean diameter in each bin + CHARACTER(100), DIMENSION(:), ALLOCATABLE :: bin_labels ! text strings for bin bounds REAL(dp) :: cmass_sum ! landscape biomass REAL(dp) :: cmass_sum_old ! landscape biomass REAL(dp) :: cheartwood_sum ! landscape biomass (heart wood) @@ -266,10 +266,10 @@ MODULE POP_Types REAL(dp) :: smoothing_buffer_cat REAL(dp) :: fire_mortality_smoothed REAL(dp) :: cat_mortality_smoothed - REAL(dp), DIMENSION(NYEAR_HISTORY) :: fire_mortality_history - REAL(dp), DIMENSION(NYEAR_HISTORY) :: cat_mortality_history - REAL(dp), DIMENSION(AGEMAX) :: freq_age ! age weighting (by age in y: 0:AGE_MAX-1) - REAL(dp), DIMENSION(AGEMAX) :: biomass_age + REAL(dp), DIMENSION(:), ALLOCATABLE :: fire_mortality_history + REAL(dp), DIMENSION(:), ALLOCATABLE :: cat_mortality_history + REAL(dp), DIMENSION(:), ALLOCATABLE :: freq_age ! age weighting (by age in y: 0:AGE_MAX-1) + REAL(dp), DIMENSION(:), ALLOCATABLE :: biomass_age END TYPE Landscape TYPE POP_TYPE @@ -292,6 +292,7 @@ MODULE POPModule USE TYPEdef, ONLY: sp, i4b USE POP_Types USE POP_Constants + !USE common_module, ONLY: handle_iostat IMPLICIT NONE @@ -340,126 +341,127 @@ SUBROUTINE ZeroPOP(POP,n) endif DO g=a, b - POP%pop_grid(g)%freq = 0.0_dp ! patch weighting - POP%pop_grid(g)%freq_old = 0.0_dp ! patch weighting - POP%pop_grid(g)%fire_freq = 0.0_dp - POP%pop_grid(g)%fire_freq_old = 0.0_dp - POP%pop_grid(g)%cat_freq = 0.0_dp - POP%pop_grid(g)%cat_freq_old = 0.0_dp - POP%pop_grid(g)%freq_ranked_age_unique = 0.0_dp - POP%pop_grid(g)%ranked_age_unique = 0 - POP%pop_grid(g)%n_age = 0 - POP%pop_grid(g)%biomass = 0.0_dp ! landscape stem biomass (weighted mean over patches) - POP%pop_grid(g)%density = 0.0_dp ! landscape tree density (weighted mean over patches) - POP%pop_grid(g)%hmean = 0.0_dp ! landscape mean treen height (weighted mean over patches) - POP%pop_grid(g)%hmax = 0.0_dp ! landscape max tree height - POP%pop_grid(g)%cmass_stem_bin = 0.0_dp ! biomass by height bin - POP%pop_grid(g)%densindiv_bin = 0.0_dp ! density by height bin - POP%pop_grid(g)%height_bin = 0.0_dp ! mean height in each bin - POP%pop_grid(g)%diameter_bin = 0.0_dp ! mean diameter in each bin - POP%pop_grid(g)%bin_labels = ' ' ! text strings for bin bounds - POP%pop_grid(g)%cmass_sum = 0.0_dp ! landscape biomass - POP%pop_grid(g)%cmass_sum_old = 0.0_dp ! landscape biomass - POP%pop_grid(g)%cheartwood_sum = 0.0_dp ! landscape biomass - POP%pop_grid(g)%csapwood_sum = 0.0_dp ! landscape biomass - POP%pop_grid(g)%csapwood_sum_old = 0.0_dp ! landscape biomass - POP%pop_grid(g)%densindiv = 0.0_dp ! landscape density of individuals - POP%pop_grid(g)%height_mean = 0.0_dp - POP%pop_grid(g)%height_max = 0.0_dp - POP%pop_grid(g)%basal_area = 0.0_dp - POP%pop_grid(g)%sapwood_loss = 0.0_dp - POP%pop_grid(g)%sapwood_area_loss = 0.0_dp - POP%pop_grid(g)%stress_mortality = 0.0_dp ! (kg C m-2 y-1) - POP%pop_grid(g)%crowding_mortality = 0.0_dp ! (kg C m-2 y-1) - POP%pop_grid(g)%fire_mortality = 0.0_dp ! (kg C m-2 y-1) - POP%pop_grid(g)%cat_mortality = 0.0_dp ! (kg C m-2 y-1) - POP%pop_grid(g)%res_mortality = 0.0_dp ! (kg C m-2 y-1) - POP%pop_grid(g)%growth = 0.0_dp - POP%pop_grid(g)%area_growth = 0.0_dp - POP%pop_grid(g)%crown_cover = 0.0_dp - POP%pop_grid(g)%crown_area = 0.0_dp - POP%pop_grid(g)%crown_volume = 0.0_dp - POP%pop_grid(g)%sapwood_area = 0.0_dp - POP%pop_grid(g)%sapwood_area_old = 0.0_dp - POP%pop_grid(g)%Kclump = 1.0_dp - POP%pop_grid(g)%npatch_active = 0 - POP%pop_grid(g)%smoothing_buffer = 0.0_dp - POP%pop_grid(g)%smoothing_buffer_cat = 0.0_dp - POP%pop_grid(g)%fire_mortality_smoothed = 0.0_dp - POP%pop_grid(g)%cat_mortality_smoothed = 0.0_dp - POP%pop_grid(g)%fire_mortality_history = 0.0_dp - POP%pop_grid(g)%cat_mortality_history = 0.0_dp - POP%pop_grid(g)%freq_age = 0.0_dp - IF (PRESENT(n)) THEN - POP%pop_grid(g)%freq_age(1) = 1.0_dp - ENDIF - POP%pop_grid(g)%biomass_age = 0.0_dp - - DO k=1, NPATCH2D - POP%pop_grid(g)%patch(k)%factor_recruit = 0.0_dp - POP%pop_grid(g)%patch(k)%pgap = 0.0_dp - POP%pop_grid(g)%patch(k)%lai = 0.0_dp - POP%pop_grid(g)%patch(k)%biomass = 0.0_dp ! total biomass in patch - POP%pop_grid(g)%patch(k)%biomass_old = 0.0_dp - POP%pop_grid(g)%patch(k)%sapwood = 0.0_dp ! total biomass in patch (sapwood) - POP%pop_grid(g)%patch(k)%heartwood = 0.0_dp ! total biomass in patch (heartwood) - POP%pop_grid(g)%patch(k)%sapwood_old = 0.0_dp - POP%pop_grid(g)%patch(k)%sapwood_area = 0.0_dp - POP%pop_grid(g)%patch(k)%sapwood_area_old = 0.0_dp - POP%pop_grid(g)%patch(k)%stress_mortality = 0.0_dp ! biomass lost in each patch due to stress - POP%pop_grid(g)%patch(k)%fire_mortality = 0.0_dp ! biomass lost in each patch due to fire partial dist - POP%pop_grid(g)%patch(k)%cat_mortality = 0.0_dp ! biomass lost in each patch due to fire partial dist - POP%pop_grid(g)%patch(k)%crowding_mortality = 0.0_dp - POP%pop_grid(g)%patch(k)%cpc = 0.0_dp - POP%pop_grid(g)%patch(k)%mortality = 0.0_dp - POP%pop_grid(g)%patch(k)%sapwood_loss = 0.0_dp - POP%pop_grid(g)%patch(k)%sapwood_area_loss = 0.0_dp - POP%pop_grid(g)%patch(k)%growth = 0.0_dp ! biomass growth in each patch due stem increment - POP%pop_grid(g)%patch(k)%area_growth = 0.0_dp - POP%pop_grid(g)%patch(k)%disturbance_interval = 0 ! prescribed disturbance(s) interval for this patch - POP%pop_grid(g)%patch(k)%first_disturbance_year = 0 - POP%pop_grid(g)%patch(k)%age = 0 ! number of years since last disturbance(s) - POP%pop_grid(g)%patch(k)%id = 0 - POP%pop_grid(g)%patch(k)%frac_NPP = 0.0_dp - POP%pop_grid(g)%patch(k)%frac_respiration = 0.0_dp - POP%pop_grid(g)%patch(k)%frac_light_uptake = 0.0_dp - POP%pop_grid(g)%patch(k)%fire_top_kill_density = 0.0_dp - - DO l=1, NLAYER - POP%pop_grid(g)%patch(k)%Layer(L)%ncohort = 0 ! number of cohorts with density >0.0_dp - POP%pop_grid(g)%patch(k)%Layer(L)%biomass = 0.0_dp ! layer biomass - POP%pop_grid(g)%patch(k)%Layer(L)%density = 0.0_dp ! layer tree density - POP%pop_grid(g)%patch(k)%Layer(L)%hmean = 0.0_dp ! layer mean tree height (weighted mean over patches) - POP%pop_grid(g)%patch(k)%Layer(L)%hmax = 0.0_dp ! layer max tree height - - DO c=1, NCOHORT_MAX - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%id = 0 - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%age = 0 ! cohort age - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%biomass = 0.0_dp ! cohort biomass - ! landscape tree density (weighted mean over patches) - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%density = 0.0_dp - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%frac_resource_uptake = 0.0_dp - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%frac_light_uptake = 0.0_dp - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%frac_interception = 0.0_dp - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%frac_respiration = 0.0_dp - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%frac_NPP = 0.0_dp - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%respiration_scalar = 0.0_dp - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%crown_area = 0.0_dp - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%Pgap = 0.0_dp - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%height = 0.0_dp - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%diameter = 0.0_dp - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%sapwood = 0.0_dp ! cohort sapwood - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%heartwood = 0.0_dp ! cohort heartwood - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%sapwood_area = 0.0_dp - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%basal_area = 0.0_dp - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%LAI = 0.0_dp - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%Cleaf = 0.0_dp - POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%Croot = 0.0_dp - ENDDO ! NCOHORT_MAX - - ENDDO ! NLAYER - - ENDDO ! NPATCH2D + POP%pop_grid(g)%freq = 0.0_dp ! patch weighting + POP%pop_grid(g)%freq_old = 0.0_dp ! patch weighting + POP%pop_grid(g)%fire_freq = 0.0_dp + POP%pop_grid(g)%fire_freq_old = 0.0_dp + POP%pop_grid(g)%cat_freq = 0.0_dp + POP%pop_grid(g)%cat_freq_old = 0.0_dp + POP%pop_grid(g)%freq_ranked_age_unique = 0.0_dp + POP%pop_grid(g)%ranked_age_unique = 0 + POP%pop_grid(g)%n_age = 0 + POP%pop_grid(g)%biomass = 0.0_dp ! landscape stem biomass (weighted mean over patches) + POP%pop_grid(g)%density = 0.0_dp ! landscape tree density (weighted mean over patches) + POP%pop_grid(g)%hmean = 0.0_dp ! landscape mean treen height (weighted mean over patches) + POP%pop_grid(g)%hmax = 0.0_dp ! landscape max tree height + POP%pop_grid(g)%cmass_stem_bin = 0.0_dp ! biomass by height bin + POP%pop_grid(g)%densindiv_bin = 0.0_dp ! density by height bin + POP%pop_grid(g)%height_bin = 0.0_dp ! mean height in each bin + POP%pop_grid(g)%diameter_bin = 0.0_dp ! mean diameter in each bin + POP%pop_grid(g)%bin_labels = ' ' ! text strings for bin bounds + POP%pop_grid(g)%cmass_sum = 0.0_dp ! landscape biomass + POP%pop_grid(g)%cmass_sum_old = 0.0_dp ! landscape biomass + POP%pop_grid(g)%cheartwood_sum = 0.0_dp ! landscape biomass + POP%pop_grid(g)%csapwood_sum = 0.0_dp ! landscape biomass + POP%pop_grid(g)%csapwood_sum_old = 0.0_dp ! landscape biomass + POP%pop_grid(g)%densindiv = 0.0_dp ! landscape density of individuals + POP%pop_grid(g)%height_mean = 0.0_dp + POP%pop_grid(g)%height_max = 0.0_dp + POP%pop_grid(g)%basal_area = 0.0_dp + POP%pop_grid(g)%sapwood_loss = 0.0_dp + POP%pop_grid(g)%sapwood_area_loss = 0.0_dp + POP%pop_grid(g)%stress_mortality = 0.0_dp ! (kg C m-2 y-1) + POP%pop_grid(g)%crowding_mortality = 0.0_dp ! (kg C m-2 y-1) + POP%pop_grid(g)%fire_mortality = 0.0_dp ! (kg C m-2 y-1) + POP%pop_grid(g)%cat_mortality = 0.0_dp ! (kg C m-2 y-1) + POP%pop_grid(g)%res_mortality = 0.0_dp ! (kg C m-2 y-1) + POP%pop_grid(g)%growth = 0.0_dp + POP%pop_grid(g)%area_growth = 0.0_dp + POP%pop_grid(g)%crown_cover = 0.0_dp + POP%pop_grid(g)%crown_area = 0.0_dp + POP%pop_grid(g)%crown_volume = 0.0_dp + POP%pop_grid(g)%sapwood_area = 0.0_dp + POP%pop_grid(g)%sapwood_area_old = 0.0_dp + POP%pop_grid(g)%Kclump = 1.0_dp + POP%pop_grid(g)%npatch_active = 0 + POP%pop_grid(g)%smoothing_buffer = 0.0_dp + POP%pop_grid(g)%smoothing_buffer_cat = 0.0_dp + POP%pop_grid(g)%fire_mortality_smoothed = 0.0_dp + POP%pop_grid(g)%cat_mortality_smoothed = 0.0_dp + POP%pop_grid(g)%fire_mortality_history = 0.0_dp + POP%pop_grid(g)%cat_mortality_history = 0.0_dp + POP%pop_grid(g)%freq_age = 0.0_dp + IF (PRESENT(n)) THEN + POP%pop_grid(g)%freq_age(1) = 1.0_dp + ENDIF + POP%pop_grid(g)%biomass_age = 0.0_dp + + DO k=1, NPATCH2D + POP%pop_grid(g)%patch(k)%factor_recruit = 0.0_dp + POP%pop_grid(g)%patch(k)%pgap = 0.0_dp + POP%pop_grid(g)%patch(k)%lai = 0.0_dp + POP%pop_grid(g)%patch(k)%biomass = 0.0_dp ! total biomass in patch + POP%pop_grid(g)%patch(k)%biomass_old = 0.0_dp + POP%pop_grid(g)%patch(k)%sapwood = 0.0_dp ! total biomass in patch (sapwood) + POP%pop_grid(g)%patch(k)%heartwood = 0.0_dp ! total biomass in patch (heartwood) + POP%pop_grid(g)%patch(k)%sapwood_old = 0.0_dp + POP%pop_grid(g)%patch(k)%sapwood_area = 0.0_dp + POP%pop_grid(g)%patch(k)%sapwood_area_old = 0.0_dp + POP%pop_grid(g)%patch(k)%stress_mortality = 0.0_dp ! biomass lost in each patch due to stress + POP%pop_grid(g)%patch(k)%fire_mortality = 0.0_dp ! biomass lost in each patch due to fire partial dist + POP%pop_grid(g)%patch(k)%cat_mortality = 0.0_dp ! biomass lost in each patch due to fire partial dist + POP%pop_grid(g)%patch(k)%crowding_mortality = 0.0_dp + POP%pop_grid(g)%patch(k)%cpc = 0.0_dp + POP%pop_grid(g)%patch(k)%mortality = 0.0_dp + POP%pop_grid(g)%patch(k)%sapwood_loss = 0.0_dp + POP%pop_grid(g)%patch(k)%sapwood_area_loss = 0.0_dp + POP%pop_grid(g)%patch(k)%growth = 0.0_dp ! biomass growth in each patch due stem increment + POP%pop_grid(g)%patch(k)%area_growth = 0.0_dp + POP%pop_grid(g)%patch(k)%disturbance_interval = 0 ! prescribed disturbance(s) interval for this patch + POP%pop_grid(g)%patch(k)%first_disturbance_year = 0 + POP%pop_grid(g)%patch(k)%age = 0 ! number of years since last disturbance(s) + POP%pop_grid(g)%patch(k)%id = 0 + POP%pop_grid(g)%patch(k)%frac_NPP = 0.0_dp + POP%pop_grid(g)%patch(k)%frac_respiration = 0.0_dp + POP%pop_grid(g)%patch(k)%frac_light_uptake = 0.0_dp + POP%pop_grid(g)%patch(k)%fire_top_kill_density = 0.0_dp + + ! Iterate over layers + DO L = 1, NLayer + POP%pop_grid(g)%patch(k)%layer(L)%ncohort = 0 ! number of cohorts with density >0.0_dp + POP%pop_grid(g)%patch(k)%layer(L)%biomass = 0.0_dp ! layer biomass + POP%pop_grid(g)%patch(k)%layer(L)%density = 0.0_dp ! layer tree density + POP%pop_grid(g)%patch(k)%layer(L)%hmean = 0.0_dp ! layer mean tree height (weighted mean over patches) + POP%pop_grid(g)%patch(k)%layer(L)%hmax = 0.0_dp ! layer max tree height + + ! Iterate over cohorts + DO c=1, NCOHORT_MAX + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%id = 0 + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%age = 0 ! cohort age + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%biomass = 0.0_dp ! cohort biomass + ! landscape tree density (weighted mean over patches) + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%density = 0.0_dp + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%frac_resource_uptake = 0.0_dp + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%frac_light_uptake = 0.0_dp + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%frac_interception = 0.0_dp + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%frac_respiration = 0.0_dp + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%frac_NPP = 0.0_dp + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%respiration_scalar = 0.0_dp + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%crown_area = 0.0_dp + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%Pgap = 0.0_dp + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%height = 0.0_dp + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%diameter = 0.0_dp + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%sapwood = 0.0_dp ! cohort sapwood + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%heartwood = 0.0_dp ! cohort heartwood + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%sapwood_area = 0.0_dp + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%basal_area = 0.0_dp + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%LAI = 0.0_dp + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%Cleaf = 0.0_dp + POP%pop_grid(g)%patch(k)%layer(L)%cohort(c)%Croot = 0.0_dp + END DO ! Finish iterating over cohorts + END DO ! Finish iterating over layers + + ENDDO ! NPATCH2D ENDDO ! pop_grid%np @@ -471,7 +473,7 @@ END SUBROUTINE ZeroPOP SUBROUTINE InitPOP2D_Poisson(POP, mean_disturbance_interval, m) ! Initialises vector of patches with maximum age correpondding to 95% of pdf - ! Starting year: uniform distribution up to maximum age + ! Starting year uniform distribution up to maximum age IMPLICIT NONE @@ -625,7 +627,7 @@ SUBROUTINE POPStep(POP, StemNPP, disturbance_interval, disturbance_intensity,LAI REAL(dp), INTENT(IN), OPTIONAL :: frac_intensity1(:), precip(:) REAL(dp), INTENT(IN), OPTIONAL :: StemNPP_pot(:) - INTEGER(i4b) :: idisturb,np,g + INTEGER(i4b) :: idisturb,np,g, P INTEGER(i4b), allocatable :: it(:) !INTEGER, INTENT(IN) :: wlogn @@ -634,9 +636,26 @@ SUBROUTINE POPStep(POP, StemNPP, disturbance_interval, disturbance_intensity,LAI np = SIZE(POP%POP_grid) allocate(it(np)) - do g=1, np - it(g) = maxval(pop%pop_grid(g)%patch(:)%age(1)) + 1 - enddo + ! I think this gets the maximum age in cell? + ! So we can initialise all the values to 0 and do repeated checking + it(:) = 0 + + IterateCells: DO g = 1, np + IteratePatches: DO P = 1, NPATCH2D + IF (POP%POP_grid(g)%patch(P)%age(1) > it(g)) THEN + it(g) = POP%POP_grid(p)%patch(P)%age(1) + END IF + END DO IteratePatches + END DO IterateCells + + it(:) = it(:) + 1 + + ! Original code for reference + !do g=1, np + !it(g) = maxval(pop%pop_grid(g)%patch(:)%age(1)) + 1 + !enddo + + ! DO idisturb = 1,NDISTURB ! CALL GetUniqueAgeFrequencies(POP, disturbance_interval, idisturb) ! ENDDO @@ -744,25 +763,25 @@ SUBROUTINE PatchAnnualDynamics(pop, StemNPP, NPPtoGPP, it, StemNPP_pot, precip) if (NPATCH2D >1 .and. it(j) > 1 .and. RESOURCE_SWITCH>0) then DO k=1,NPATCH2D freq = pop%pop_grid(j)%freq(pop%pop_grid(j)%patch(k)%id) - nc = pop%pop_grid(j)%patch(k)%Layer(1)%ncohort + nc = pop%pop_grid(j)%patch(k)%layer(1)%ncohort DO c=1,nc - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_light_uptake = & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_light_uptake = & pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_interception ! defined in terms of Pgap ! total autotrophic resp, summed over all cohorts and patches tmp_respiration = tmp_respiration + & - freq*pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%respiration_scalar + freq*pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%respiration_scalar ENDDO tmp_light = tmp_light + freq*(1.0_dp - pop%pop_grid(j)%patch(k)%Pgap) ENDDO IF (tmp_respiration .gt. 1.0e-8_dp .and. tmp_light .gt. 1.0e-8_dp) then DO k=1,NPATCH2D ! fraction respiration and un-normalised NPP for each patch - nc = pop%pop_grid(j)%patch(k)%Layer(1)%ncohort + nc = pop%pop_grid(j)%patch(k)%layer(1)%ncohort freq = pop%pop_grid(j)%freq(pop%pop_grid(j)%patch(k)%id) ! frac autotrophic resp pop%pop_grid(j)%patch(k)%frac_respiration = & - sum(pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:nc)%respiration_scalar)/tmp_respiration + sum(pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:nc)%respiration_scalar)/tmp_respiration ! frac gpp pop%pop_grid(j)%patch(k)%frac_light_uptake = & @@ -801,11 +820,11 @@ SUBROUTINE PatchAnnualDynamics(pop, StemNPP, NPPtoGPP, it, StemNPP_pot, precip) tmp_light = 0.0_dp tmp_respiration = 0.0_dp tmp_fracNPP = 0.0_dp - if (pop%pop_grid(j)%patch(k)%Layer(1)%ncohort>1) then - nc = pop%pop_grid(j)%patch(k)%Layer(1)%ncohort - DO c=1,pop%pop_grid(j)%patch(k)%Layer(1)%ncohort + if (pop%pop_grid(j)%patch(k)%layer(1)%ncohort>1) then + nc = pop%pop_grid(j)%patch(k)%layer(1)%ncohort + DO c=1,pop%pop_grid(j)%patch(k)%layer(1)%ncohort - cmass_stem = pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%biomass + cmass_stem = pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass densindiv = pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density ! get initial basal area @@ -823,85 +842,85 @@ SUBROUTINE PatchAnnualDynamics(pop, StemNPP, NPPtoGPP, it, StemNPP_pot, precip) crown_area = densindiv*PI*(((k_allom1 * diam ** k_rp )/PI)**0.5_dp)**2 endif - tmp = tmp + (pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%biomass/ & ! sum over all cohorts - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%density)**POWERbiomass * & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%density + tmp = tmp + (pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass/ & ! sum over all cohorts + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density)**POWERbiomass * & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density - tmp_light = tmp_light + pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_interception + tmp_light = tmp_light + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_interception - tmp_respiration = tmp_respiration + pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%respiration_scalar + tmp_respiration = tmp_respiration + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%respiration_scalar - tmp2(c) = sum((pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c:nc)%biomass/ & ! sum over all cohorts c:nc - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c:nc)%density)**POWERbiomass * & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c:nc)%density) + tmp2(c) = sum((pop%pop_grid(j)%patch(k)%layer(1)%cohort(c:nc)%biomass/ & ! sum over all cohorts c:nc + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c:nc)%density)**POWERbiomass * & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c:nc)%density) ENDDO ! un-normalised fractional gross resource uptake: weighted combination of components ! where cohort competes with older cohorts and where it does not - DO c=1,pop%pop_grid(j)%patch(k)%Layer(1)%ncohort + DO c=1,pop%pop_grid(j)%patch(k)%layer(1)%ncohort if (RESOURCE_SWITCH ==1) then - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_interception = & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_interception/tmp_light + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_interception = & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_interception/tmp_light else - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_interception = 1.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_interception = 1.0_dp endif ENDDO !normalised fractional gross resource uptake - nc = pop%pop_grid(j)%patch(k)%Layer(1)%ncohort - DO c=1,pop%pop_grid(j)%patch(k)%Layer(1)%ncohort + nc = pop%pop_grid(j)%patch(k)%layer(1)%ncohort + DO c=1,pop%pop_grid(j)%patch(k)%layer(1)%ncohort !normalised fractional gross resource uptake - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_light_uptake = & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_interception/ & - sum(pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:nc)%frac_interception) + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_light_uptake = & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_interception/ & + sum(pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:nc)%frac_interception) ENDDO ! fraction respiration and un-normalised NPP - DO c=1,pop%pop_grid(j)%patch(k)%Layer(1)%ncohort - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_respiration = & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%respiration_scalar/tmp_respiration + DO c=1,pop%pop_grid(j)%patch(k)%layer(1)%ncohort + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_respiration = & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%respiration_scalar/tmp_respiration - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_NPP = & - max(pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_light_uptake*(1.0_dp/NPPtoGPP(j)) - & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_respiration*(1.0_dp/NPPtoGPP(j)-1.0_dp),0.0_dp) + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_NPP = & + max(pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_light_uptake*(1.0_dp/NPPtoGPP(j)) - & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_respiration*(1.0_dp/NPPtoGPP(j)-1.0_dp),0.0_dp) - tmp_fracNPP = tmp_fracNPP + pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_NPP + tmp_fracNPP = tmp_fracNPP + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_NPP ENDDO ! normalised fraction NPP - DO c=1,pop%pop_grid(j)%patch(k)%Layer(1)%ncohort - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_NPP = & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_NPP/tmp_fracNPP + DO c=1,pop%pop_grid(j)%patch(k)%layer(1)%ncohort + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_NPP = & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_NPP/tmp_fracNPP ENDDO ! fraction net resource uptake - DO c=1,pop%pop_grid(j)%patch(k)%Layer(1)%ncohort + DO c=1,pop%pop_grid(j)%patch(k)%layer(1)%ncohort if (RESOURCE_SWITCH==0) then ! default net fraction resource uptake - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_resource_uptake = & - (pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%biomass/ & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%density)**POWERbiomass * & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%density/tmp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_resource_uptake = & + (pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass/ & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density)**POWERbiomass * & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density/tmp elseif (RESOURCE_SWITCH==1) then ! fraction net resource uptake = fraction NPP - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_resource_uptake = & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_NPP * & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_resource_uptake = & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_NPP * & pop%pop_grid(j)%patch(k)%frac_NPP endif @@ -910,14 +929,14 @@ SUBROUTINE PatchAnnualDynamics(pop, StemNPP, NPPtoGPP, it, StemNPP_pot, precip) else c = 1 - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_NPP = 1 + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_NPP = 1 tmp_fracNPP = 1 - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_respiration = 1 - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_light_uptake =1 - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_resource_uptake = 1 + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_respiration = 1 + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_light_uptake =1 + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_resource_uptake = 1 if (RESOURCE_SWITCH==1) then - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_resource_uptake = pop%pop_grid(j)%patch(k)%frac_NPP + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_resource_uptake = pop%pop_grid(j)%patch(k)%frac_NPP endif endif @@ -934,11 +953,11 @@ SUBROUTINE PatchAnnualDynamics(pop, StemNPP, NPPtoGPP, it, StemNPP_pot, precip) pop%pop_grid(j)%patch(k)%biomass_old = pop%pop_grid(j)%patch(k)%biomass pop%pop_grid(j)%patch(k)%growth = 0.0_dp pop%pop_grid(j)%patch(k)%area_growth = 0.0_dp - nc = pop%pop_grid(j)%patch(k)%Layer(1)%ncohort + nc = pop%pop_grid(j)%patch(k)%layer(1)%ncohort freq = pop%pop_grid(j)%freq(pop%pop_grid(j)%patch(k)%id) - DO c=1,pop%pop_grid(j)%patch(k)%Layer(1)%ncohort + DO c=1,pop%pop_grid(j)%patch(k)%layer(1)%ncohort - cmass_stem = pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%biomass + cmass_stem = pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass densindiv = pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density ! get initial basal area @@ -950,12 +969,12 @@ SUBROUTINE PatchAnnualDynamics(pop, StemNPP, NPPtoGPP, it, StemNPP_pot, precip) ENDIF ! increment biomass in cohort - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%biomass = & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%biomass + & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_resource_uptake*StemNPP(j,1) + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass = & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass + & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_resource_uptake*StemNPP(j,1) - cmass_stem = pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%biomass - tmp = tmp + freq*pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_resource_uptake + cmass_stem = pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass + tmp = tmp + freq*pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_resource_uptake ! get incremented basal area IF ( PRESENT(precip) ) THEN @@ -967,29 +986,29 @@ SUBROUTINE PatchAnnualDynamics(pop, StemNPP, NPPtoGPP, it, StemNPP_pot, precip) ! increment sapwood in cohort - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood = & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood + & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_resource_uptake*StemNPP(j,1) + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood = & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood + & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_resource_uptake*StemNPP(j,1) ! increment heartwood in cohort - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%heartwood = & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%heartwood + & - ksapwood*pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%heartwood = & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%heartwood + & + ksapwood*pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood ! keep track of patch-level sapwood loss pop%pop_grid(j)%patch(k)%sapwood_loss = pop%pop_grid(j)%patch(k)%sapwood_loss + & - ksapwood*pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood + ksapwood*pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood ! decrease sapwood in cohort (accounting for loss to heartwood) - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood = & - (1.0_dp - ksapwood)*pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood = & + (1.0_dp - ksapwood)*pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood - !if ( pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%density.gt.1e-9) then + !if ( pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density.gt.1e-9) then ! patch biomass increment pop%pop_grid(j)%patch(k)%growth = pop%pop_grid(j)%patch(k)%growth + & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_resource_uptake*StemNPP(j,1) + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_resource_uptake*StemNPP(j,1) ! patch sapwood area increment pop%pop_grid(j)%patch(k)%area_growth = pop%pop_grid(j)%patch(k)%area_growth + & @@ -998,13 +1017,13 @@ SUBROUTINE PatchAnnualDynamics(pop, StemNPP, NPPtoGPP, it, StemNPP_pot, precip) ! endif ! increment cohort age - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%age = & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%age + 1 + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%age = & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%age + 1 ENDDO - ! Layer biomass (summed over cohorts) - nc = pop%pop_grid(j)%patch(k)%Layer(1)%ncohort - pop%pop_grid(j)%patch(k)%Layer(1)%biomass = SUM(pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:nc)%biomass) + ! layer biomass (summed over cohorts) + nc = pop%pop_grid(j)%patch(k)%layer(1)%ncohort + pop%pop_grid(j)%patch(k)%layer(1)%biomass = SUM(pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:nc)%biomass) ENDDO @@ -1023,17 +1042,17 @@ SUBROUTINE PatchAnnualDynamics(pop, StemNPP, NPPtoGPP, it, StemNPP_pot, precip) pop%pop_grid(j)%patch(k)%crowding_mortality = 0.0_dp pop%pop_grid(j)%patch(k)%mortality = 0.0_dp crown_area = 0.0_dp - DO c=1,pop%pop_grid(j)%patch(k)%Layer(1)%ncohort - cmass_stem = pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%biomass - cmass_stem_inc=StemNPP(j,1)*pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_resource_uptake + DO c=1,pop%pop_grid(j)%patch(k)%layer(1)%ncohort + cmass_stem = pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass + cmass_stem_inc=StemNPP(j,1)*pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_resource_uptake if (present(StemNPP_pot)) then - growth_efficiency=pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_resource_uptake* & + growth_efficiency=pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%frac_resource_uptake* & StemNPP_pot(j) /(cmass_stem**(POWERGrowthEfficiency)) else growth_efficiency=cmass_stem_inc/(cmass_stem**(POWERGrowthEfficiency)) endif - ! growth_efficiency=cmass_stem_inc/(pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood**(POWERGrowthEfficiency)) + ! growth_efficiency=cmass_stem_inc/(pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood**(POWERGrowthEfficiency)) mort=MORT_MAX/(1.0_dp+(growth_efficiency/GROWTH_EFFICIENCY_MIN)**Pmort) ! mort = 0 ! test @@ -1043,10 +1062,10 @@ SUBROUTINE PatchAnnualDynamics(pop, StemNPP, NPPtoGPP, it, StemNPP_pot, precip) IF (pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%diameter*100_dp .GT. 1.0_dp) THEN if (ALLOM_SWITCH.eq.1) then ! assumes crown radius (m) = 0.14 * dbh (cm) - crown_area = crown_area + pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%density* & + crown_area = crown_area + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density* & PI*(pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%diameter*100.0_dp*0.14_dp)**2 else - crown_area = crown_area + pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%density* & + crown_area = crown_area + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density* & k_allom1 * pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%diameter ** k_rp endif ELSE @@ -1070,57 +1089,57 @@ SUBROUTINE PatchAnnualDynamics(pop, StemNPP, NPPtoGPP, it, StemNPP_pot, precip) mort*cmass_stem pop%pop_grid(j)%patch(k)%sapwood_loss = pop%pop_grid(j)%patch(k)%sapwood_loss + & - mort*pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood + mort*pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood pop%pop_grid(j)%patch(k)%sapwood_area_loss = pop%pop_grid(j)%patch(k)%sapwood_area_loss + & - mort*pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood_area + mort*pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood_area - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%biomass = cmass_stem*(1.0_dp-mort) - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%heartwood = & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%heartwood *(1.0_dp-mort) - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood = & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood *(1.0_dp-mort) + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass = cmass_stem*(1.0_dp-mort) + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%heartwood = & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%heartwood *(1.0_dp-mort) + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood = & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood *(1.0_dp-mort) - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%density = & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%density*(1.0_dp-mort) - IF (pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%density.LT.DENSINDIV_MIN) THEN + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density = & + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density*(1.0_dp-mort) + IF (pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density.LT.DENSINDIV_MIN) THEN ! remove cohort - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%density = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%id = 0 + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%id = 0 pop%pop_grid(j)%patch(k)%stress_mortality = pop%pop_grid(j)%patch(k)%stress_mortality + & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%biomass + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass pop%pop_grid(j)%patch(k)%sapwood_loss = pop%pop_grid(j)%patch(k)%sapwood_loss + & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood pop%pop_grid(j)%patch(k)%sapwood_area_loss = pop%pop_grid(j)%patch(k)%sapwood_area_loss + & - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood_area - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%biomass = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%heartwood = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood_area + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%heartwood = 0.0_dp ELSE - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%id = 1 + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%id = 1 !COMMLN Why is id here 1 instead of c or sth useful? Call it differently nc = nc+1 ivec(nc)=c ENDIF ENDDO ! SHUFFLE if necessary to remove zero-density cohorts - IF (nc.LT.pop%pop_grid(j)%patch(k)%Layer(1)%ncohort) THEN - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:nc)=pop%pop_grid(j)%patch(k)%Layer(1)%cohort(ivec(1:nc)) - pop%pop_grid(j)%patch(k)%Layer(1)%ncohort = nc - - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%density = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%id = 0 - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%biomass = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%sapwood = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%sapwood_area = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%basal_area = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%heartwood = 0.0_dp + IF (nc.LT.pop%pop_grid(j)%patch(k)%layer(1)%ncohort) THEN + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:nc)=pop%pop_grid(j)%patch(k)%layer(1)%cohort(ivec(1:nc)) + pop%pop_grid(j)%patch(k)%layer(1)%ncohort = nc + + pop%pop_grid(j)%patch(k)%layer(1)%cohort(nc+1:NCOHORT_MAX)%density = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(nc+1:NCOHORT_MAX)%id = 0 + pop%pop_grid(j)%patch(k)%layer(1)%cohort(nc+1:NCOHORT_MAX)%biomass = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(nc+1:NCOHORT_MAX)%sapwood = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(nc+1:NCOHORT_MAX)%sapwood_area = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(nc+1:NCOHORT_MAX)%basal_area = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(nc+1:NCOHORT_MAX)%heartwood = 0.0_dp ENDIF - ! Layer biomass (summed over cohorts) - nc = pop%pop_grid(j)%patch(k)%Layer(1)%ncohort - pop%pop_grid(j)%patch(k)%Layer(1)%biomass = SUM(pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:nc)%biomass) + ! layer biomass (summed over cohorts) + nc = pop%pop_grid(j)%patch(k)%layer(1)%ncohort + pop%pop_grid(j)%patch(k)%layer(1)%biomass = SUM(pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:nc)%biomass) ENDDO ENDDO @@ -1156,7 +1175,7 @@ SUBROUTINE GetUniqueAgeFrequencies(pop, disturbance_interval, idisturb) TYPE(POP_TYPE), INTENT(INOUT) :: POP INTEGER(i4b), INTENT(IN) :: disturbance_interval(:,:), idisturb - INTEGER(i4b) :: g, i,j,k,agecopy,idcopy + INTEGER(i4b) :: g, i,j,k,agecopy,idcopy, P REAL(dp), ALLOCATABLE :: midpoint(:) INTEGER(i4b), ALLOCATABLE :: ranked_age(:), ranked_age_init(:) INTEGER(i4b) :: age_tmp @@ -1164,7 +1183,7 @@ SUBROUTINE GetUniqueAgeFrequencies(pop, disturbance_interval, idisturb) REAL(dp), ALLOCATABLE :: tmp(:), freq_tmp(:), freq_tmp1(:) REAL(dp) :: freq INTEGER(i4b) :: n_age ! number of unique ages - INTEGER(i4b) :: npatch_active ! number of active patches + !INTEGER(i4b) :: NPATCH2D ! number of active patches REAL(dp):: disturbance_freq INTEGER(i4b) :: i_max, Poisson_age(1000), np REAL(dp):: CumPoisson_weight(1000) @@ -1176,23 +1195,27 @@ SUBROUTINE GetUniqueAgeFrequencies(pop, disturbance_interval, idisturb) np = SIZE(POP%POP_grid) DO g=1,np - npatch_active = NPATCH2D - IF (.NOT.ALLOCATED(midpoint)) ALLOCATE(midpoint(npatch_active)) - IF (.NOT.ALLOCATED(counter)) ALLOCATE(counter(npatch_active)) - IF (.NOT.ALLOCATED(ranked_age)) ALLOCATE(ranked_age(npatch_active)) - IF (.NOT.ALLOCATED(ranked_age_init)) ALLOCATE(ranked_age_init(npatch_active)) - IF (.NOT.ALLOCATED(ranked_age_id)) ALLOCATE(ranked_age_id(npatch_active)) - IF (.NOT.ALLOCATED(ranked_age_unique_id)) ALLOCATE(ranked_age_unique_id(npatch_active)) - IF (.NOT.ALLOCATED(tmp)) ALLOCATE(tmp(npatch_active)) - IF (.NOT.ALLOCATED(freq_tmp)) ALLOCATE(freq_tmp(npatch_active)) - IF (.NOT.ALLOCATED(freq_tmp1)) ALLOCATE(freq_tmp1(npatch_active)) + IF (.NOT.ALLOCATED(midpoint)) ALLOCATE(midpoint(NPATCH2D)) + IF (.NOT.ALLOCATED(counter)) ALLOCATE(counter(NPATCH2D)) + IF (.NOT.ALLOCATED(ranked_age)) ALLOCATE(ranked_age(NPATCH2D)) + IF (.NOT.ALLOCATED(ranked_age_init)) ALLOCATE(ranked_age_init(NPATCH2D)) + IF (.NOT.ALLOCATED(ranked_age_id)) ALLOCATE(ranked_age_id(NPATCH2D)) + IF (.NOT.ALLOCATED(ranked_age_unique_id)) ALLOCATE(ranked_age_unique_id(NPATCH2D)) + IF (.NOT.ALLOCATED(tmp)) ALLOCATE(tmp(NPATCH2D)) + IF (.NOT.ALLOCATED(freq_tmp)) ALLOCATE(freq_tmp(NPATCH2D)) + IF (.NOT.ALLOCATED(freq_tmp1)) ALLOCATE(freq_tmp1(NPATCH2D)) ! rank patches in order of age pop%pop_grid(g)%ranked_age_unique(:, idisturb) = 0 - ranked_age_init = pop%pop_grid(g)%patch%age(idisturb) - ranked_age = pop%pop_grid(g)%patch%age(idisturb) - ranked_age_id = pop%pop_grid(g)%patch%id + ! Previous code required compile-time knowledge of array sizes + + IteratePatches: DO P = 1, NPATCH2D + ranked_age_init(P) = pop%pop_grid(g)%patch(P)%age(idisturb) + ranked_age(P) = pop%pop_grid(g)%patch(P)%age(idisturb) + ranked_age_id(P) = pop%pop_grid(g)%patch(P)%id + END DO IteratePatches + ranked_age_unique_id = 0 freq_tmp = 0.0_dp freq = 0.0_dp @@ -1200,8 +1223,8 @@ SUBROUTINE GetUniqueAgeFrequencies(pop, disturbance_interval, idisturb) midpoint = 0.0_dp - DO i = 1, npatch_active -1 - DO j = i+1, npatch_active + DO i = 1, NPATCH2D -1 + DO j = i+1, NPATCH2D IF (ranked_age(i).GT.ranked_age(j)) THEN agecopy = ranked_age(i) idcopy = ranked_age_id(i) @@ -1216,7 +1239,7 @@ SUBROUTINE GetUniqueAgeFrequencies(pop, disturbance_interval, idisturb) ! subset to unique ages k=0 age_tmp = -1 - DO i = 1, npatch_active + DO i = 1, NPATCH2D IF (ranked_age(i).NE.age_tmp) k = k+1 pop%pop_grid(g)%ranked_age_unique(k, idisturb) = ranked_age(i) ranked_age_unique_id(k) = ranked_age_id(i) @@ -1280,7 +1303,7 @@ SUBROUTINE GetUniqueAgeFrequencies(pop, disturbance_interval, idisturb) ENDDO ENDDO - pop%pop_grid(g)%freq_ranked_age_unique(1:npatch_active,idisturb) = freq_tmp + pop%pop_grid(g)%freq_ranked_age_unique(1:NPATCH2D,idisturb) = freq_tmp pop%pop_grid(g)%n_age(idisturb) = n_age DEALLOCATE (bound) @@ -1300,7 +1323,7 @@ SUBROUTINE GetPatchFrequencies(pop) TYPE(POP_TYPE), INTENT(INOUT) :: POP - INTEGER(i4b) :: n1, n2, g, REPCOUNT, np, idist + INTEGER(i4b) :: n1, n2, g, np, idist, P, repcount REAL(dp) :: sum_freq np = SIZE(Pop%pop_grid) @@ -1310,30 +1333,71 @@ SUBROUTINE GetPatchFrequencies(pop) DO idist = 1, NDISTURB IF (idist.EQ.1) THEN DO n1=1,pop%pop_grid(g)%n_age(1) - repcount = COUNT(pop%pop_grid(g)%patch(:)%age(1).EQ.pop%pop_grid(g)%ranked_age_unique(n1,1)) - WHERE (pop%pop_grid(g)%patch(:)%age(1).EQ.pop%pop_grid(g)%ranked_age_unique(n1,1)) - pop%pop_grid(g)%freq = pop%pop_grid(g)%freq_ranked_age_unique(n1,1) /REAL(repcount,dp) - ENDWHERE + repcount = 0 + ! count number of occurrences of specific age I think? + ! This loop replaces the COUNT function + IteratePatches_1: DO P = 1, NPATCH2D + IF (POP%POP_grid(g)%patch(P)%age(1) == POP%POP_grid(g)%ranked_age_unique(n1, 1)) THEN + repcount = repcount + 1 + END IF + END DO IteratePatches_1 + + ! Use number of occurrences to average? + ! This loop replaces the WHERE block + IteratePatches_2: DO P = 1, NPATCH2D + IF (POP%POP_grid(g)%patch(P)%age(1) == POP%POP_grid(g)%ranked_age_unique(n1, 1)) THEN + POP%POP_grid(g)%freq(P) = POP%POP_grid(g)%freq_ranked_age_unique(n1, 1) / REAL(repcount, dp) + END IF + END DO IteratePatches_2 + + ! Original code for reference during PR + !repcount = COUNT(pop%pop_grid(g)%patch(:)%age(1).EQ.pop%pop_grid(g)%ranked_age_unique(n1,1)) + !WHERE (pop%pop_grid(g)%patch(:)%age(1).EQ.pop%pop_grid(g)%ranked_age_unique(n1,1)) + !pop%pop_grid(g)%freq = pop%pop_grid(g)%freq_ranked_age_unique(n1,1) /REAL(repcount,dp) + !ENDWHERE ENDDO ELSEIF (idist.EQ.2) THEN ! first calculate weights for patches with age(2)>age(1) DO n1=1,pop%pop_grid(g)%n_age(1) - - DO n2=1,pop%pop_grid(g)%n_age(idist) - repcount = COUNT((pop%pop_grid(g)%patch(1:NPATCH)%age(1) .EQ. & - pop%pop_grid(g)%ranked_age_unique(n1,1)).AND. & - (pop%pop_grid(g)%patch(1:NPATCH)%age(idist) .EQ. & - pop%pop_grid(g)%ranked_age_unique(n2,idist))) - WHERE ((pop%pop_grid(g)%patch(1:NPATCH)%age(1).EQ.pop%pop_grid(g)%ranked_age_unique(n1,1)).AND. & - (pop%pop_grid(g)%patch(1:NPATCH)%age(idist).EQ.pop%pop_grid(g)%ranked_age_unique(n2,idist))) - pop%pop_grid(g)%freq(1:NPATCH) = pop%pop_grid(g)%freq_ranked_age_unique(n1,1)* & - pop%pop_grid(g)%freq_ranked_age_unique(n2,idist) & - /REAL(repcount,dp) - - ENDWHERE + repcount = 0 + ! This loop replaces the COUNT function + IteratePatches_3: DO P = 1, NPATCH2D + IF ((POP%POP_grid(g)%patch(P)%age(1) ==& + POP%POP_grid(g)%ranked_age_unique(n1, 1)) .AND.& + (POP%POP_grid(g)%patch(P)%age(idist) ==& + POP%POP_grid(g)%ranked_age_unique(n2, idist))) THEN + repcount = repcount + 1 + END IF + END DO IteratePatches_3 + + ! Replaces the WHERE block + IteratePatches_4: DO P = 1, NPATCH2D + IF ((POP%POP_grid(g)%patch(P)%age(1) ==& + POP%POP_grid(g)%ranked_age_unique(n1, 1)) .AND.& + (POP%POP_grid(g)%patch(P)%age(idist) ==& + POP%POP_grid(g)%ranked_age_unique(n2, idist))) THEN + + POP%POP_grid(g)%freq(P) = & + POP%pop_grid(g)%freq_ranked_age_unique(n1, 1) *& + POP%pop_grid(g)%freq_ranked_age_unique(n2, idist) /& + REAL(repcount, dp) + END IF + + END DO IteratePatches_4 + !repcount = COUNT((pop%pop_grid(g)%patch(1:NPATCH)%age(1) .EQ. & + !pop%pop_grid(g)%ranked_age_unique(n1,1)).AND. & + !(pop%pop_grid(g)%patch(1:NPATCH)%age(idist) .EQ. & + !pop%pop_grid(g)%ranked_age_unique(n2,idist))) + !WHERE ((pop%pop_grid(g)%patch(1:NPATCH)%age(1).EQ.pop%pop_grid(g)%ranked_age_unique(n1,1)).AND. & + !(pop%pop_grid(g)%patch(1:NPATCH)%age(idist).EQ.pop%pop_grid(g)%ranked_age_unique(n2,idist))) + !pop%pop_grid(g)%freq(1:NPATCH) = pop%pop_grid(g)%freq_ranked_age_unique(n1,1)* & + !pop%pop_grid(g)%freq_ranked_age_unique(n2,idist) & + !/REAL(repcount,dp) + + !ENDWHERE ENDDO @@ -1854,7 +1918,7 @@ SUBROUTINE Patch_partial_disturb(pop,idisturb,intensity,frac_intensity1) Psurvival*pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%heartwood pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density = & Psurvival*pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density - IF (pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%density.LT.DENSINDIV_MIN) THEN + IF (pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density.LT.DENSINDIV_MIN) THEN ! remove cohort pop%pop_grid(j)%patch(k)%fire_mortality = pop%pop_grid(j)%patch(k)%fire_mortality + & pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass @@ -1862,31 +1926,31 @@ SUBROUTINE Patch_partial_disturb(pop,idisturb,intensity,frac_intensity1) pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood pop%pop_grid(j)%patch(k)%sapwood_area_loss = pop%pop_grid(j)%patch(k)%sapwood_area_loss + & pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood_area - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%density = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%id = 0 - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%biomass = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%heartwood = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood_area = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%id = 0 + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%heartwood = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood_area = 0.0_dp ELSE - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%id = 1 + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%id = 1 nc = nc+1 ivec(nc)=c ENDIF ENDDO ! SHUFFLE if necessary to remove zero-density cohorts - IF (nc.LT.pop%pop_grid(j)%patch(k)%Layer(1)%ncohort) THEN - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:nc)=pop%pop_grid(j)%patch(k)%Layer(1)%cohort(ivec(1:nc)) - pop%pop_grid(j)%patch(k)%Layer(1)%ncohort = nc - - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%density = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%id = 0 - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%biomass = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%sapwood = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%sapwood_area = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%heartwood = 0.0_dp + IF (nc.LT.pop%pop_grid(j)%patch(k)%layer(1)%ncohort) THEN + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:nc)=pop%pop_grid(j)%patch(k)%layer(1)%cohort(ivec(1:nc)) + pop%pop_grid(j)%patch(k)%layer(1)%ncohort = nc + + pop%pop_grid(j)%patch(k)%layer(1)%cohort(nc+1:NCOHORT_MAX)%density = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(nc+1:NCOHORT_MAX)%id = 0 + pop%pop_grid(j)%patch(k)%layer(1)%cohort(nc+1:NCOHORT_MAX)%biomass = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(nc+1:NCOHORT_MAX)%sapwood = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(nc+1:NCOHORT_MAX)%sapwood_area = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(nc+1:NCOHORT_MAX)%heartwood = 0.0_dp ENDIF pop%pop_grid(j)%patch(k)%age(idisturb) = 0 @@ -1921,9 +1985,9 @@ SUBROUTINE Patch_partial_disturb2(pop,idisturb) DO j=1,np DO k=1,NPATCH2D pop%pop_grid(j)%patch(k)%cat_mortality = 0.0_dp - ! Layer biomass (summed over cohorts) - nc = pop%pop_grid(j)%patch(k)%Layer(1)%ncohort - pop%pop_grid(j)%patch(k)%Layer(1)%biomass = SUM(pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:nc)%biomass) + ! layer biomass (summed over cohorts) + nc = pop%pop_grid(j)%patch(k)%layer(1)%ncohort + pop%pop_grid(j)%patch(k)%layer(1)%biomass = SUM(pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:nc)%biomass) IF (((pop%pop_grid(j)%patch(k)%first_disturbance_year(idisturb).NE.0).AND. & (pop%pop_grid(j)%patch(k)%first_disturbance_year(idisturb).EQ.pop%pop_grid(j)%patch(k)%age(idisturb))).OR. & @@ -1937,8 +2001,8 @@ SUBROUTINE Patch_partial_disturb2(pop,idisturb) pop%pop_grid(j)%patch(k)%cat_mortality = 0.0_dp DO c = 1, pop%pop_grid(j)%patch(k)%layer(1)%ncohort ! kill fraction of each cohort, up to 80% of patch biomass - if (pop%pop_grid(j)%patch(k)%cat_mortality < 0.8_dp * pop%pop_grid(j)%patch(k)%Layer(1)%biomass ) then - Pmort = min( (0.8_dp*pop%pop_grid(j)%patch(k)%Layer(1)%biomass-pop%pop_grid(j)%patch(k)%fire_mortality) & + if (pop%pop_grid(j)%patch(k)%cat_mortality < 0.8_dp * pop%pop_grid(j)%patch(k)%layer(1)%biomass ) then + Pmort = min( (0.8_dp*pop%pop_grid(j)%patch(k)%layer(1)%biomass-pop%pop_grid(j)%patch(k)%fire_mortality) & /pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass, 1.0_dp) else Pmort = 0.0_dp @@ -1959,7 +2023,7 @@ SUBROUTINE Patch_partial_disturb2(pop,idisturb) Psurvival*pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%heartwood pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density = & Psurvival*pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density - IF (pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%density.LT.DENSINDIV_MIN) THEN + IF (pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density.LT.DENSINDIV_MIN) THEN ! remove cohort pop%pop_grid(j)%patch(k)%cat_mortality = pop%pop_grid(j)%patch(k)%cat_mortality + & pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass @@ -1967,31 +2031,31 @@ SUBROUTINE Patch_partial_disturb2(pop,idisturb) pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood pop%pop_grid(j)%patch(k)%sapwood_area_loss = pop%pop_grid(j)%patch(k)%sapwood_area_loss + & pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood_area - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%density = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%id = 0 - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%biomass = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%heartwood = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood_area = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%id = 0 + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%heartwood = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood_area = 0.0_dp ELSE - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%id = 1 + pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%id = 1 nc = nc+1 ivec(nc)=c ENDIF ENDDO ! SHUFFLE if necessary to remove zero-density cohorts - IF (nc.LT.pop%pop_grid(j)%patch(k)%Layer(1)%ncohort) THEN - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:nc)=pop%pop_grid(j)%patch(k)%Layer(1)%cohort(ivec(1:nc)) - pop%pop_grid(j)%patch(k)%Layer(1)%ncohort = nc - - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%density = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%id = 0 - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%biomass = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%sapwood = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%sapwood_area = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%heartwood = 0.0_dp + IF (nc.LT.pop%pop_grid(j)%patch(k)%layer(1)%ncohort) THEN + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:nc)=pop%pop_grid(j)%patch(k)%layer(1)%cohort(ivec(1:nc)) + pop%pop_grid(j)%patch(k)%layer(1)%ncohort = nc + + pop%pop_grid(j)%patch(k)%layer(1)%cohort(nc+1:NCOHORT_MAX)%density = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(nc+1:NCOHORT_MAX)%id = 0 + pop%pop_grid(j)%patch(k)%layer(1)%cohort(nc+1:NCOHORT_MAX)%biomass = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(nc+1:NCOHORT_MAX)%sapwood = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(nc+1:NCOHORT_MAX)%sapwood_area = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(nc+1:NCOHORT_MAX)%heartwood = 0.0_dp ENDIF pop%pop_grid(j)%patch(k)%age(idisturb) = 0 @@ -2037,29 +2101,29 @@ SUBROUTINE Patch_disturb(pop,idisturb,precip) SUM(pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:nc)%sapwood) pop%pop_grid(j)%patch(k)%sapwood_area_loss = pop%pop_grid(j)%patch(k)%sapwood_area_loss + & SUM(pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:nc)%sapwood_area) - pop%pop_grid(j)%patch(k)%layer(1:NLayer)%ncohort = 0 - pop%pop_grid(j)%patch(k)%layer(1:NLayer)%biomass = 0.0_dp - pop%pop_grid(j)%patch(k)%layer(1:NLayer)%density = 0.0_dp - pop%pop_grid(j)%patch(k)%layer(1:NLayer)%hmean = 0.0_dp - pop%pop_grid(j)%patch(k)%layer(1:NLayer)%hmax = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%ncohort = 0 + pop%pop_grid(j)%patch(k)%layer(1)%biomass = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%density = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%hmean = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%hmax = 0.0_dp pop%pop_grid(j)%patch(k)%age(idisturb) = 0 pop%pop_grid(j)%patch(k)%first_disturbance_year(idisturb) = 0 - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%density = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%id = 0 - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%biomass = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%sapwood = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%sapwood_area = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%heartwood = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%age = 0 - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%lai = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%height = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%pgap = 1.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_resource_uptake = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_light_uptake = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_interception = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_respiration = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_NPP = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%density = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%id = 0 + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%biomass = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%sapwood = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%sapwood_area = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%heartwood = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%age = 0 + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%lai = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%height = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%pgap = 1.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%frac_resource_uptake = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%frac_light_uptake = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%frac_interception = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%frac_respiration = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%frac_NPP = 0.0_dp pop%pop_grid(j)%patch(k)%area_growth = 0.0_dp pop%pop_grid(j)%patch(k)%pgap = 1.0_dp ! understorey recruitment @@ -2078,29 +2142,29 @@ SUBROUTINE Patch_disturb(pop,idisturb,precip) pop%pop_grid(j)%patch(k)%sapwood_area_loss = pop%pop_grid(j)%patch(k)%sapwood_area_loss + & SUM(pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:nc)%sapwood_area) pop%pop_grid(j)%patch(k)%cat_mortality = SUM(pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:nc)%biomass) - pop%pop_grid(j)%patch(k)%layer(1:NLayer)%ncohort = 0 - pop%pop_grid(j)%patch(k)%layer(1:NLayer)%biomass = 0.0_dp - pop%pop_grid(j)%patch(k)%layer(1:NLayer)%density = 0.0_dp - pop%pop_grid(j)%patch(k)%layer(1:NLayer)%hmean = 0.0_dp - pop%pop_grid(j)%patch(k)%layer(1:NLayer)%hmax = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%ncohort = 0 + pop%pop_grid(j)%patch(k)%layer(1)%biomass = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%density = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%hmean = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%hmax = 0.0_dp pop%pop_grid(j)%patch(k)%age(idisturb) = 0 pop%pop_grid(j)%patch(k)%first_disturbance_year(idisturb) = 0 - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%density = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%id = 0 - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%biomass = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%sapwood = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%sapwood_area = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%heartwood = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%age = 0 - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%lai = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%height = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%pgap = 1.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_resource_uptake = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_light_uptake = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_interception = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_respiration = 0.0_dp - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_NPP = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%density = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%id = 0 + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%biomass = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%sapwood = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%sapwood_area = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%heartwood = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%age = 0 + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%lai = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%height = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%pgap = 1.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%frac_resource_uptake = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%frac_light_uptake = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%frac_interception = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%frac_respiration = 0.0_dp + pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:NCOHORT_MAX)%frac_NPP = 0.0_dp pop%pop_grid(j)%patch(k)%area_growth = 0.0_dp pop%pop_grid(j)%patch(k)%pgap = 1.0_dp ! understorey recruitment @@ -2138,7 +2202,7 @@ SUBROUTINE layer_recruitment(pop,precip) DO j=1,np DO k=1,NPATCH2D IF (RECRUIT_SWITCH==0) THEN - pop%pop_grid(j)%patch(k)%factor_recruit = EXP(-0.6_dp*((pop%pop_grid(j)%patch(k)%Layer(1)%biomass)**(0.6667_dp))) + pop%pop_grid(j)%patch(k)%factor_recruit = EXP(-0.6_dp*((pop%pop_grid(j)%patch(k)%layer(1)%biomass)**(0.6667_dp))) ELSEIF (RECRUIT_SWITCH==1) THEN pop%pop_grid(j)%patch(k)%factor_recruit = max(pop%pop_grid(j)%patch(k)%pgap,1.0e-3_dp) ENDIF @@ -2149,17 +2213,17 @@ SUBROUTINE layer_recruitment(pop,precip) densindiv=DENSINDIV_MAX*mu + pop%pop_grid(j)%patch(k)%fire_top_kill_density cmass=CMASS_STEM_INIT*densindiv/DENSINDIV_MAX - !write(5599,*), pop%pop_grid(j)%patch(k)%fire_top_kill_density, densindiv, pop%pop_grid(j)%patch(k)%Layer(1)%ncohort + !write(5599,*), pop%pop_grid(j)%patch(k)%fire_top_kill_density, densindiv, pop%pop_grid(j)%patch(k)%layer(1)%ncohort !COMMLN below: should not be cohort +1 or .LE. ! IF (cmass>EPS*10.0_dp .AND. densindiv>DENSINDIV_MIN .AND. & - (pop%pop_grid(j)%patch(k)%Layer(1)%ncohort+1).LT.NCOHORT_MAX) THEN + (pop%pop_grid(j)%patch(k)%layer(1)%ncohort+1).LT.NCOHORT_MAX) THEN ! create a new cohort - pop%pop_grid(j)%patch(k)%Layer(1)%ncohort = pop%pop_grid(j)%patch(k)%Layer(1)%ncohort + 1 - ncohort = pop%pop_grid(j)%patch(k)%Layer(1)%ncohort - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(ncohort)%biomass = cmass - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(ncohort)%density = densindiv - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(ncohort)%sapwood = cmass + pop%pop_grid(j)%patch(k)%layer(1)%ncohort = pop%pop_grid(j)%patch(k)%layer(1)%ncohort + 1 + ncohort = pop%pop_grid(j)%patch(k)%layer(1)%ncohort + pop%pop_grid(j)%patch(k)%layer(1)%cohort(ncohort)%biomass = cmass + pop%pop_grid(j)%patch(k)%layer(1)%cohort(ncohort)%density = densindiv + pop%pop_grid(j)%patch(k)%layer(1)%cohort(ncohort)%sapwood = cmass IF ( PRESENT(precip) ) THEN CALL GET_ALLOMETRY( ALLOM_SWITCH, cmass, densindiv, ht, diam, basal, precip(j)) @@ -2197,7 +2261,7 @@ SUBROUTINE layer_recruitment_single_patch(pop, index, grid_index,precip) DO j=grid_index,grid_index DO k=index,index IF (RECRUIT_SWITCH==0) THEN - pop%pop_grid(j)%patch(k)%factor_recruit = EXP(-0.6_dp*((pop%pop_grid(j)%patch(k)%Layer(1)%biomass)**(0.6667_dp))) + pop%pop_grid(j)%patch(k)%factor_recruit = EXP(-0.6_dp*((pop%pop_grid(j)%patch(k)%layer(1)%biomass)**(0.6667_dp))) ELSEIF (RECRUIT_SWITCH==1) THEN !pop%pop_grid(j)%patch(k)%factor_recruit = pop%pop_grid(j)%patch(k)%pgap pop%pop_grid(j)%patch(k)%factor_recruit = 1 @@ -2208,13 +2272,13 @@ SUBROUTINE layer_recruitment_single_patch(pop, index, grid_index,precip) cmass=CMASS_STEM_INIT*densindiv/DENSINDIV_MAX IF (cmass>EPS*10.0_dp .AND. densindiv>DENSINDIV_MIN .AND. & - (pop%pop_grid(j)%patch(k)%Layer(1)%ncohort+1).LT.NCOHORT_MAX) THEN + (pop%pop_grid(j)%patch(k)%layer(1)%ncohort+1).LT.NCOHORT_MAX) THEN ! create a new cohort - pop%pop_grid(j)%patch(k)%Layer(1)%ncohort = pop%pop_grid(j)%patch(k)%Layer(1)%ncohort + 1 - ncohort = pop%pop_grid(j)%patch(k)%Layer(1)%ncohort - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(ncohort)%biomass = cmass - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(ncohort)%density = densindiv - pop%pop_grid(j)%patch(k)%Layer(1)%cohort(ncohort)%sapwood = cmass + pop%pop_grid(j)%patch(k)%layer(1)%ncohort = pop%pop_grid(j)%patch(k)%layer(1)%ncohort + 1 + ncohort = pop%pop_grid(j)%patch(k)%layer(1)%ncohort + pop%pop_grid(j)%patch(k)%layer(1)%cohort(ncohort)%biomass = cmass + pop%pop_grid(j)%patch(k)%layer(1)%cohort(ncohort)%density = densindiv + pop%pop_grid(j)%patch(k)%layer(1)%cohort(ncohort)%sapwood = cmass IF ( PRESENT(precip) ) THEN CALL GET_ALLOMETRY( ALLOM_SWITCH, cmass, densindiv, ht, diam, basal, precip(j)) @@ -2446,11 +2510,11 @@ SUBROUTINE INTERPOLATE_BIOMASS_1D(pop, disturbance_interval,it,g) INTEGER(i4b) :: nage,iage, i_min, i_max INTEGER(i4b) :: i_min_growth, i_max_growth - REAL(dp) :: disturbance_freq, tmp_min, tmp_max, tmp1_min, tmp1_max, tmp_array(NPATCH2D) - REAL(dp) :: tmp2_min, tmp2_max - REAL(dp) :: tmp3_min, tmp3_max - REAL(dp) :: tmp4_min, tmp4_max - LOGICAL :: MASK(NPATCH2D) + REAL(dp) :: disturbance_freq + REAL(dp) :: biomass_min, sapwood_min, sapwood_area_min, stress_mort_min,& + crowd_mort_min, growth_min + REAL(dp) :: biomass_max, sapwood_max, sapwood_area_max, stress_mort_max,& + crowd_mort_max, growth_max INTEGER(i4b) :: age_min, age_max INTEGER(i4b) :: age_min_growth, age_max_growth INTEGER(i4b), ALLOCATABLE :: age(:) @@ -2458,6 +2522,10 @@ SUBROUTINE INTERPOLATE_BIOMASS_1D(pop, disturbance_interval,it,g) REAL(dp), ALLOCATABLE ::csapwood_age(:), sapwood_area_age(:), growth_age(:) REAL(dp), ALLOCATABLE ::freq_age(:) + ! Patch iterator + INTEGER :: P + REAL(dp) :: NPoints + ! get interpolated biomass,sapwood, stress mortality, crowding mortality, disturbance mortality POP%pop_grid(g)%cmass_sum= 0.0_dp POP%pop_grid(g)%stress_mortality = 0.0_dp @@ -2465,7 +2533,6 @@ SUBROUTINE INTERPOLATE_BIOMASS_1D(pop, disturbance_interval,it,g) pop%pop_grid(g)%crowding_mortality = 0.0_dp pop%pop_grid(g)%csapwood_sum = 0.0_dp pop%pop_grid(g)%sapwood_area = 0.0_dp - tmp_array = 0.0_dp nage = min(POP%pop_grid(g)%patch(1)%disturbance_interval(1),it)+1 ! maximum age !nage = maxval(pop%pop_grid(g)%patch(:)%age(1)) @@ -2493,8 +2560,6 @@ SUBROUTINE INTERPOLATE_BIOMASS_1D(pop, disturbance_interval,it,g) !pop%pop_grid(g)%biomass_age(2:agemax) = pop%pop_grid(g)%biomass_age(1:agemax-1) !pop%pop_grid(g)%biomass_age(1) = 0.0 !cmass_age = pop%pop_grid(g)%biomass_age - tmp_min = 0.0_dp - tmp_max = 0.0_dp pop%pop_grid(g)%biomass_age = 0.0_dp @@ -2514,161 +2579,373 @@ SUBROUTINE INTERPOLATE_BIOMASS_1D(pop, disturbance_interval,it,g) if (sum(freq_age)>0.0_dp) freq_age = freq_age/sum(freq_age) DO iage = 1, nage + ! Can simplify logic a lot here- we just want the largest age that's less + ! than age(iage), and smallest age that's greater than age(iage) + age_min = 0 + i_min = 0 + age_max = 0 + i_max = 0 + IteratePatches: DO P = 1, NPATCH2D + IF ((POP%POP_grid(g)%patch(P)%age(1) <= age(iage)) .AND.& + (POP%POP_grid(g)%patch(P)%age(1) > age_min)) THEN + age_min = POP%POP_grid(g)%patch(P)%age(1) + i_min = P + END IF + + ! We also need to catch the first iteration here- the rest of the code + ! is configured to assume age_max=0 and i_max=0 when there are no + ! cohorts with age greater than age(iage), so we can't just start with + ! age_max=HUGE(age_max), i_max=HUGE(i_max) + IF ((POP%POP_grid(g)%patch(P)%age(1) >= age(iage)) .AND.& + ((POP%POP_grid(g)%patch(P)%age(1) < age_max) .OR.& + (age_max == 0))) THEN + age_max = POP%POP_grid(g)%patch(P)%age(1) + i_max = P + END IF + END DO IteratePatches + ! get nearest ages bracketing age(iage) - if (any(pop%pop_grid(g)%patch(:)%age(1).LE.age(iage))) then - age_min = MAXVAL(pop%pop_grid(g)%patch(:)%age(1), & - pop%pop_grid(g)%patch(:)%age(1).LE.age(iage)) - i_min = MAXLOC(pop%pop_grid(g)%patch(:)%age(1), 1, & - pop%pop_grid(g)%patch(:)%age(1).LE.age(iage)) - else - age_min = 0 - i_min = 0 - endif - if (any(pop%pop_grid(g)%patch(:)%age(1).GE.age(iage))) then - age_max = MINVAL(pop%pop_grid(g)%patch(:)%age(1), & - pop%pop_grid(g)%patch(:)%age(1).GE.age(iage)) - i_max = MINLOC(pop%pop_grid(g)%patch(:)%age(1), 1, & - pop%pop_grid(g)%patch(:)%age(1).GE.age(iage)) - else - age_max = 0 - i_max = 0 - endif + !if (any(pop%pop_grid(g)%patch(:)%age(1).LE.age(iage))) then + !age_min = MAXVAL(pop%pop_grid(g)%patch(:)%age(1), & + !pop%pop_grid(g)%patch(:)%age(1).LE.age(iage)) + !i_min = MAXLOC(pop%pop_grid(g)%patch(:)%age(1), 1, & + !pop%pop_grid(g)%patch(:)%age(1).LE.age(iage)) + !else + !age_min = 0 + !i_min = 0 + !endif + !if (any(pop%pop_grid(g)%patch(:)%age(1).GE.age(iage))) then + !age_max = MINVAL(pop%pop_grid(g)%patch(:)%age(1), & + !pop%pop_grid(g)%patch(:)%age(1).GE.age(iage)) + !i_max = MINLOC(pop%pop_grid(g)%patch(:)%age(1), 1, & + !pop%pop_grid(g)%patch(:)%age(1).GE.age(iage)) + !else + !age_max = 0 + !i_max = 0 + !endif age_min_growth = age_min age_max_growth = age_max i_min_growth = i_min i_max_growth = i_max + ! Set the values to be accumulated to 0.0 + cmass_age(iage) = 0.0 + growth_age(iage) = 0.0 + csapwood_age(iage) = 0.0 + sapwood_area_age(iage) = 0.0 + stress_mort_age(iage) = 0.0 + crowd_mort_age(iage) = 0.0 + + growth_min = 0.0 + biomass_min = 0.0 + sapwood_min = 0.0 + sapwood_area_min = 0.0 + stress_mort_min = 0.0 + crowd_mort_min = 0.0 + + growth_max = 0.0 + biomass_max = 0.0 + sapwood_max = 0.0 + sapwood_area_max = 0.0 + stress_mort_max = 0.0 + crowd_mort_max = 0.0 + if ((i_min.gt.0).and.(i_max.gt.0).and.(age_max.eq.age_min)) then ! no need to interpolate - MASK = pop%pop_grid(g)%patch(:)%age(1).eq.age_min - where (MASK) - tmp_array = 1.0_dp - elsewhere - tmp_array = 0.0_dp - endwhere - cmass_age(iage) = & - SUM(pop%pop_grid(g)%patch(:)%layer(1)%biomass,MASK)/SUM(tmp_array) - growth_age(iage) = & - SUM(pop%pop_grid(g)%patch(:)%growth,MASK)/SUM(tmp_array) - csapwood_age(iage) = SUM(pop%pop_grid(g)%patch(:)%sapwood,MASK)/SUM(tmp_array) - sapwood_area_age(iage) = & - SUM(pop%pop_grid(g)%patch(:)%sapwood_area,MASK)/SUM(tmp_array) - stress_mort_age(iage)= & - SUM(pop%pop_grid(g)%patch(:)%stress_mortality,MASK)/SUM(tmp_array) - crowd_mort_age(iage)= & - SUM(pop%pop_grid(g)%patch(:)%crowding_mortality,MASK)/SUM(tmp_array) + ! We want to find all the patches where the age equals the minimum + ! age, then do averages over those patches. So accumulate while + ! keeping track of number of points (as double for later division) + Npoints = 0.0 + IteratePatches_2: DO P = 1, NPATCH2D + IF (POP%POP_grid(g)%patch(P)%age(1) == age_min) THEN + + cmass_age(iage) = cmass_age(iage) +& + POP%POP_grid(g)%patch(P)%layer(1)%biomass + growth_age(iage) = growth_age(iage) +& + POP%POP_grid(g)%patch(P)%growth + csapwood_age(iage) = csapwood_age(iage) +& + POP%POP_grid(g)%patch(P)%sapwood + sapwood_area_age(iage) = sapwood_area_age(iage) +& + POP%POP_grid(g)%patch(P)%sapwood_area + stress_mort_age(iage) = stress_mort_age(iage) +& + POP%POP_grid(g)%patch(P)%stress_mortality + crowd_mort_age(iage) = crowd_mort_age(iage) +& + POP%POP_grid(g)%patch(P)%crowding_mortality + + Npoints = Npoints + 1 + END IF + END DO IteratePatches_2 + + ! Scale by number of patches that meet condition + cmass_age(iage) = cmass_age(iage) / Npoints + growth_age(iage) = growth_age(iage) / Npoints + csapwood_age(iage) = csapwood_age(iage) / Npoints + sapwood_area_age(iage) = sapwood_area_age(iage) / Npoints + stress_mort_age(iage) = stress_mort_age(iage) / Npoints + crowd_mort_age(iage) = crowd_mort_age(iage) / Npoints + + !MASK = pop%pop_grid(g)%patch(:)%age(1).eq.age_min + !where (MASK) + !tmp_array = 1.0_dp + !elsewhere + !tmp_array = 0.0_dp + !endwhere + !cmass_age(iage) = & + !SUM(pop%pop_grid(g)%patch(:)%layer(1)%biomass,MASK)/SUM(tmp_array) + !growth_age(iage) = & + !SUM(pop%pop_grid(g)%patch(:)%growth,MASK)/SUM(tmp_array) + !csapwood_age(iage) = SUM(pop%pop_grid(g)%patch(:)%sapwood,MASK)/SUM(tmp_array) + !sapwood_area_age(iage) = & + !SUM(pop%pop_grid(g)%patch(:)%sapwood_area,MASK)/SUM(tmp_array) + !stress_mort_age(iage)= & + !SUM(pop%pop_grid(g)%patch(:)%stress_mortality,MASK)/SUM(tmp_array) + !crowd_mort_age(iage)= & + !SUM(pop%pop_grid(g)%patch(:)%crowding_mortality,MASK)/SUM(tmp_array) else ! interpolate or extrapolate if ((i_min.eq.0).and.(i_max.gt.0)) then - ! interpolate to zero - age_min = 0 - i_min = 0 + ! Case where no patches with age less than iage + ! interpolate to zero + age_min = 0 + i_min = 0 elseif ((i_max.eq.0).and.(i_min.gt.0)) then - ! extrapolate to higher age - age_max = age_min - i_max = i_min - age_min = MAXVAL(pop%pop_grid(g)%patch(:)%age(1), & - pop%pop_grid(g)%patch(:)%age(1).LT.age_max) - i_min = MAXLOC(pop%pop_grid(g)%patch(:)%age(1),1, & - pop%pop_grid(g)%patch(:)%age(1).LT.age_max) + ! Case where no patches with age more than iage + ! extrapolate to higher age + age_max = age_min + i_max = i_min + IteratePatches_3: DO P = 1, NPATCH2D + IF (POP%POP_grid(g)%patch(P)%age(1) > age_min) THEN + age_min = POP%POP_grid(g)%patch(P)%age(1) + i_min = P + END IF + END DO IteratePatches_3 + !age_min = MAXVAL(pop%pop_grid(g)%patch(:)%age(1), & + !pop%pop_grid(g)%patch(:)%age(1).LT.age_max) + !i_min = MAXLOC(pop%pop_grid(g)%patch(:)%age(1),1, & + !pop%pop_grid(g)%patch(:)%age(1).LT.age_max) endif ! interpolate or extrapolate (growth) + ! i_min_growth and i_max_growth were set to i_min and i_max, why not + ! reuse prior logic checks? if ((i_min_growth.eq.0).and.(i_max_growth.gt.0).and.age(iage).LE.2) then - ! interpolate to zero - age_min_growth = 0 - i_min_growth = 0 + ! interpolate to zero + age_min_growth = 0 + i_min_growth = 0 elseif (((age_min_growth.LE.2).OR.(i_min_growth.eq.0)).and. & - (i_max_growth.gt.0).and.age(iage).GT.2) then - ! extrapolate to lower age - age_min_growth = age_max_growth - i_min_growth = i_max_growth - - age_max_growth = MINVAL(pop%pop_grid(g)%patch(:)%age(1), & - pop%pop_grid(g)%patch(:)%age(1).GT.age_min_growth) - i_max_growth = MINLOC(pop%pop_grid(g)%patch(:)%age(1), 1, & - pop%pop_grid(g)%patch(:)%age(1).GT.age_min_growth) + (i_max_growth.gt.0).and.age(iage).GT.2) then + ! I don't understand precisely what's going on here- what's this + ! logic checking for physically? + ! extrapolate to lower age + age_min_growth = age_max_growth + i_min_growth = i_max_growth + + IteratePatches_4: DO P = 1, NPATCH2D + IF ((POP%POP_grid(g)%patch(P)%age(1) > age_min_growth) .AND.& + (POP%POP_grid(g)%patch(P)%age(1) < age_max_growth)) THEN + + age_max_growth = POP%POP_grid(g)%patch(P)%age(1) + i_max_growth = P + END IF + END DO IteratePatches_4 + + !age_max_growth = MINVAL(pop%pop_grid(g)%patch(:)%age(1), & + !pop%pop_grid(g)%patch(:)%age(1).GT.age_min_growth) + !i_max_growth = MINLOC(pop%pop_grid(g)%patch(:)%age(1), 1, & + !pop%pop_grid(g)%patch(:)%age(1).GT.age_min_growth) elseif ((i_max_growth.eq.0).and.(i_min_growth.gt.0)) then - ! extrapolate to higher age - age_max_growth = age_min_growth - i_max_growth = i_min_growth - - age_min_growth = MAXVAL(pop%pop_grid(g)%patch(:)%age(1), & - pop%pop_grid(g)%patch(:)%age(1).LT.age_max_growth) - i_min_growth = MAXLOC(pop%pop_grid(g)%patch(:)%age(1),1, & - pop%pop_grid(g)%patch(:)%age(1).LT.age_max_growth) + ! extrapolate to higher age + age_max_growth = age_min_growth + i_max_growth = i_min_growth + + IteratePatches_5: DO P = 1, NPATCH2D + IF ((POP%POP_grid(g)%patch(P)%age(1) < age_max_growth) .AND.& + (POP%POP_grid(g)%patch(P)%age(1) > age_min_growth)) THEN + + age_min_growth = POP%POP_grid(g)%patch(P)%age(1) + i_min_growth = P + END IF + END DO IteratePatches_5 + + !age_min_growth = MAXVAL(pop%pop_grid(g)%patch(:)%age(1), & + !pop%pop_grid(g)%patch(:)%age(1).LT.age_max_growth) + !i_min_growth = MAXLOC(pop%pop_grid(g)%patch(:)%age(1),1, & + !pop%pop_grid(g)%patch(:)%age(1).LT.age_max_growth) endif if (i_min.ne.0.and.age_min.ne.0) then - MASK = pop%pop_grid(g)%patch(:)%age(1).eq.age_min - where (MASK) - tmp_array = 1.0_dp - elsewhere - tmp_array = 0.0_dp - endwhere - tmp_min = SUM(pop%pop_grid(g)%patch(:)%layer(1)%biomass,MASK)/SUM(tmp_array) - tmp1_min = SUM(pop%pop_grid(g)%patch(:)%stress_mortality,MASK)/SUM(tmp_array) - tmp2_min = SUM(pop%pop_grid(g)%patch(:)%crowding_mortality,MASK)/SUM(tmp_array) - tmp3_min = SUM(pop%pop_grid(g)%patch(:)%sapwood,MASK)/SUM(tmp_array) - tmp4_min = SUM(pop%pop_grid(g)%patch(:)%sapwood_area,MASK)/SUM(tmp_array) - else - tmp_min = 0.0_dp - tmp1_min = 0.0_dp - tmp2_min = 0.0_dp - tmp3_min = 0.0_dp - tmp4_min = 0.0_dp - endif + ! Do the minimum values + Npoints = 0.0 + IteratePatches_6: DO P = 1, NPATCH2D + IF (POP%POP_grid(g)%patch(P)%age(1) == age_min) THEN + + biomass_min = biomass_min +& + POP%POP_grid(g)%patch(P)%layer(1)%biomass + sapwood_min = sapwood_min +& + POP%POP_grid(g)%patch(P)%sapwood + sapwood_area_min = sapwood_area_min +& + POP%POP_grid(g)%patch(P)%sapwood_area + stress_mort_min = stress_mort_min +& + POP%POP_grid(g)%patch(P)%stress_mortality + crowd_mort_min = crowd_mort_min +& + POP%POP_grid(g)%patch(P)%crowding_mortality + + Npoints = Npoints + 1.0 + END IF + END DO IteratePatches_6 + + biomass_min = biomass_min / Npoints + sapwood_min = sapwood_min / Npoints + sapwood_area_min = sapwood_area_min / Npoints + stress_mort_min = stress_mort_min / Npoints + crowd_mort_min = crowd_mort_min / Npoints - MASK = pop%pop_grid(g)%patch(:)%age(1).eq.age_max - where (MASK) - tmp_array = 1.0_dp - elsewhere - tmp_array = 0.0_dp - endwhere - tmp_max = SUM(pop%pop_grid(g)%patch(:)%layer(1)%biomass,MASK)/SUM(tmp_array) - tmp1_max = SUM(pop%pop_grid(g)%patch(:)%stress_mortality,MASK)/SUM(tmp_array) - tmp2_max = SUM(pop%pop_grid(g)%patch(:)%crowding_mortality,MASK)/SUM(tmp_array) - tmp3_max = SUM(pop%pop_grid(g)%patch(:)%sapwood,MASK)/SUM(tmp_array) - tmp4_max = SUM(pop%pop_grid(g)%patch(:)%sapwood_area,MASK)/SUM(tmp_array) - - cmass_age(iage) = tmp_min + (tmp_max-tmp_min)/real(age_max-age_min,dp)* & - real(age(iage)-age_min,dp) + ELSE - stress_mort_age(iage) = tmp1_min + (tmp1_max-tmp1_min)/real(age_max-age_min,dp)* & - real(age(iage)-age_min,dp) + biomass_min = 0.0_dp + sapwood_min = 0.0_dp + sapwood_area_min = 0.0_dp + stress_mort_min = 0.0_dp + crowd_mort_min = 0.0_dp + + END IF + !MASK = pop%pop_grid(g)%patch(:)%age(1).eq.age_min + !where (MASK) + !tmp_array = 1.0_dp + !elsewhere + !tmp_array = 0.0_dp + !endwhere + !tmp_min = SUM(pop%pop_grid(g)%patch(:)%layer(1)%biomass,MASK)/SUM(tmp_array) + !tmp1_min = SUM(pop%pop_grid(g)%patch(:)%stress_mortality,MASK)/SUM(tmp_array) + !tmp2_min = SUM(pop%pop_grid(g)%patch(:)%crowding_mortality,MASK)/SUM(tmp_array) + !tmp3_min = SUM(pop%pop_grid(g)%patch(:)%sapwood,MASK)/SUM(tmp_array) + !tmp4_min = SUM(pop%pop_grid(g)%patch(:)%sapwood_area,MASK)/SUM(tmp_array) + !else + !tmp_min = 0.0_dp + !tmp1_min = 0.0_dp + !tmp2_min = 0.0_dp + !tmp3_min = 0.0_dp + !tmp4_min = 0.0_dp + !endif + + ! Do the maximum values + Npoints = 0.0 + IteratePatches_7: DO P = 1, NPATCH2D + IF (POP%POP_grid(g)%patch(P)%age(1) == age_max) THEN + + biomass_max = biomass_max +& + POP%POP_grid(g)%patch(P)%layer(1)%biomass + sapwood_max = sapwood_max +& + POP%POP_grid(g)%patch(P)%sapwood + sapwood_area_max = sapwood_area_max +& + POP%POP_grid(g)%patch(P)%sapwood_area + stress_mort_max = stress_mort_max +& + POP%POP_grid(g)%patch(P)%stress_mortality + crowd_mort_max = crowd_mort_max +& + POP%POP_grid(g)%patch(P)%crowding_mortality + + Npoints = Npoints + 1.0 + END IF + END DO IteratePatches_7 + + biomass_max = biomass_max / Npoints + sapwood_max = sapwood_max / Npoints + sapwood_area_max = sapwood_area_max / Npoints + stress_mort_max = stress_mort_max / Npoints + crowd_mort_max = crowd_mort_max / Npoints + + !MASK = pop%pop_grid(g)%patch(:)%age(1).eq.age_max + !where (MASK) + !tmp_array = 1.0_dp + !elsewhere + !tmp_array = 0.0_dp + !endwhere + !tmp_max = SUM(pop%pop_grid(g)%patch(:)%layer(1)%biomass,MASK)/SUM(tmp_array) + !tmp1_max = SUM(pop%pop_grid(g)%patch(:)%stress_mortality,MASK)/SUM(tmp_array) + !tmp2_max = SUM(pop%pop_grid(g)%patch(:)%crowding_mortality,MASK)/SUM(tmp_array) + !tmp3_max = SUM(pop%pop_grid(g)%patch(:)%sapwood,MASK)/SUM(tmp_array) + !tmp4_max = SUM(pop%pop_grid(g)%patch(:)%sapwood_area,MASK)/SUM(tmp_array) + + ! Now perform linear interpolation + cmass_age(iage) = biomass_min + (biomass_max-biomass_min) *& + REAL(age(iage) - age_min, dp) / REAL(age_max-age_min, dp) + + csapwood_age(iage) = sapwood_min + (sapwood_max-sapwood_min) *& + REAL(age(iage) - age_min, dp) / REAL(age_max-age_min, dp) + + sapwood_area_age(iage) = sapwood_area_min +& + (sapwood_area_max-sapwood_area_min) *& + REAL(age(iage) - age_min, dp) / REAL(age_max-age_min, dp) + + stress_mort_age(iage) = stress_mort_min +& + (stress_mort_max-stress_mort_min) *& + REAL(age(iage) - age_min, dp) / REAL(age_max-age_min, dp) + + crowd_mort_age(iage) = crowd_mort_min + (crowd_mort_max-crowd_mort_min) *& + REAL(age(iage) - age_min, dp) / REAL(age_max-age_min, dp) + + !stress_mort_age(iage) = tmp1_min + (tmp1_max-tmp1_min)/real(age_max-age_min,dp)* & + !real(age(iage)-age_min,dp) + + !crowd_mort_age(iage) = tmp2_min + (tmp2_max-tmp2_min)/real(age_max-age_min,dp)* & + !real(age(iage)-age_min,dp) + + !csapwood_age(iage) = tmp3_min + (tmp3_max-tmp3_min)/real(age_max-age_min,dp)* & + !real(age(iage)-age_min,dp) + + !sapwood_area_age(iage) = tmp4_min + (tmp4_max-tmp4_min)/real(age_max-age_min,dp)* & + !real(age(iage)-age_min,dp) + + ! Now do the growth + ! Minimum growth + if (i_min_growth.ne.0.and.age_min_growth.ne.0) then + Npoints = 0.0 + IteratePatches_8: DO P = 1, NPATCH2D + IF (POP%POP_grid(g)%patch(P)%age(1) == age_min_growth) THEN + + growth_min = growth_min + POP%POP_grid(g)%patch(P)%growth + + Npoints = Npoints + 1.0 + END IF + END DO IteratePatches_8 + + growth_min = growth_min / Npoints + !MASK = pop%pop_grid(g)%patch(:)%age(1).eq.age_min_growth + !where (MASK) + !tmp_array = 1.0_dp + !elsewhere + !tmp_array = 0.0_dp + !endwhere + !tmp_min = SUM(pop%pop_grid(g)%patch(:)%growth,MASK)/SUM(tmp_array) + else + growth_min = 0.0_dp + endif - crowd_mort_age(iage) = tmp2_min + (tmp2_max-tmp2_min)/real(age_max-age_min,dp)* & - real(age(iage)-age_min,dp) + ! Maximum growth + Npoints = 0.0 + IteratePatches_9: DO P = 1, NPATCH2D + IF (POP%POP_grid(g)%patch(P)%age(1) == age_max_growth) THEN - csapwood_age(iage) = tmp3_min + (tmp3_max-tmp3_min)/real(age_max-age_min,dp)* & - real(age(iage)-age_min,dp) + growth_max = growth_max + POP%POP_grid(g)%patch(P)%growth - sapwood_area_age(iage) = tmp4_min + (tmp4_max-tmp4_min)/real(age_max-age_min,dp)* & - real(age(iage)-age_min,dp) + Npoints = Npoints + 1.0 + END IF + END DO IteratePatches_9 - if (i_min_growth.ne.0.and.age_min_growth.ne.0) then - MASK = pop%pop_grid(g)%patch(:)%age(1).eq.age_min_growth - where (MASK) - tmp_array = 1.0_dp - elsewhere - tmp_array = 0.0_dp - endwhere - tmp_min = SUM(pop%pop_grid(g)%patch(:)%growth,MASK)/SUM(tmp_array) - else - tmp_min = 0.0_dp - endif + growth_max = growth_max / Npoints - MASK = pop%pop_grid(g)%patch(:)%age(1).eq.age_max_growth - where (MASK) - tmp_array = 1.0_dp - elsewhere - tmp_array = 0.0_dp - endwhere - tmp_max = SUM(pop%pop_grid(g)%patch(:)%growth,MASK)/SUM(tmp_array) + !MASK = pop%pop_grid(g)%patch(:)%age(1).eq.age_max_growth + !where (MASK) + !tmp_array = 1.0_dp + !elsewhere + !tmp_array = 0.0_dp + !endwhere + !tmp_max = SUM(pop%pop_grid(g)%patch(:)%growth,MASK)/SUM(tmp_array) - growth_age(iage) = tmp_min + (tmp_max-tmp_min)/real(age_max_growth-age_min_growth,dp)* & - real(age(iage)-age_min_growth,dp) + growth_age(iage) = growth_min + (growth_max - growth_min) *& + REAL(age(iage) - age_min_growth, dp) /& + real(age_max_growth-age_min_growth, dp) !cmass_age(iage) = cmass_age(iage) + growth_age(iage) - stress_mort_age(iage) - crowd_mort_age(iage) endif @@ -2721,6 +2998,7 @@ SUBROUTINE INTERPOLATE_BIOMASS_1D(pop, disturbance_interval,it,g) !!$ write(603,"(350e16.6)") real(age) !!$if (it==5) stop !!$endif + ! Why deallocate, it already deallocates when it goes out of scope? DEALLOCATE(age) DEALLOCATE(freq_age) DEALLOCATE(cmass_age) @@ -2741,18 +3019,21 @@ SUBROUTINE INTERPOLATE_FIREMORTALITY(pop, disturbance_interval,it,g) INTEGER(i4b) :: nage,iage, i_min, i_max INTEGER(i4b) :: i_min_growth, i_max_growth - REAL(dp) :: disturbance_freq,tmp_min, tmp_max, tmp_array(NPATCH2D) - REAL(dp) :: tmp5_min, tmp5_max - LOGICAL :: MASK(NPATCH2D) + REAL(dp) :: disturbance_freq, fire_mort_min, fire_mort_max INTEGER(i4b) :: age_min, age_max INTEGER(i4b) :: age_min_growth, age_max_growth INTEGER(i4b), ALLOCATABLE :: age(:) REAL(dp), ALLOCATABLE :: fire_mort_age(:) REAL(dp), ALLOCATABLE :: freq_age(:) + INTEGER :: P + REAL(dp) :: Npoints + + ! Why don't we interpolate just like all the other variables? + ! I don't see any different between the code here and the code in the other + ! interpolations ! get interpolated fire mortality POP%pop_grid(g)%fire_mortality = 0.0_dp - tmp_array = 0.0_dp nage = min(POP%pop_grid(g)%patch(1)%disturbance_interval(1),it)+1 ! maximum age IF (POP%pop_grid(g)%LU==2) then ! secondary forest @@ -2772,8 +3053,7 @@ SUBROUTINE INTERPOLATE_FIREMORTALITY(pop, disturbance_interval,it,g) IF(.NOT.ALLOCATED(freq_age)) ALLOCATE(freq_age(nage)) IF(.NOT.ALLOCATED(fire_mort_age)) ALLOCATE(fire_mort_age(nage)) - tmp_min = 0.0_dp - tmp_max = 0.0_dp + ! Why do we set this to 0? pop%pop_grid(g)%biomass_age = 0.0_dp IF (POP%pop_grid(g)%LU==2) then ! secondary forest @@ -2791,43 +3071,53 @@ SUBROUTINE INTERPOLATE_FIREMORTALITY(pop, disturbance_interval,it,g) ENDIF if (sum(freq_age)>0.0_dp) freq_age = freq_age/sum(freq_age) + fire_mort_min = 0.0 + fire_mort_max = 0.0 + DO iage = 1, nage - ! get nearest ages bracketing age(iage) - if (any(pop%pop_grid(g)%patch(:)%age(1).LE.age(iage))) then - age_min = MAXVAL(pop%pop_grid(g)%patch(:)%age(1), & - pop%pop_grid(g)%patch(:)%age(1).LE.age(iage)) - i_min = MAXLOC(pop%pop_grid(g)%patch(:)%age(1), 1, & - pop%pop_grid(g)%patch(:)%age(1).LE.age(iage)) - else - age_min = 0 - i_min = 0 - endif - if (any(pop%pop_grid(g)%patch(:)%age(1).GE.age(iage))) then - age_max = MINVAL(pop%pop_grid(g)%patch(:)%age(1), & - pop%pop_grid(g)%patch(:)%age(1).GE.age(iage)) - i_max = MINLOC(pop%pop_grid(g)%patch(:)%age(1), 1, & - pop%pop_grid(g)%patch(:)%age(1).GE.age(iage)) - else - age_max = 0 - i_max = 0 - endif + ! Can simplify logic a lot here- we just want the largest age that's less + ! than age(iage), and smallest age that's greater than age(iage) + age_min = 0 + i_min = 0 + age_max = 0 + i_max = 0 + IteratePatches: DO P = 1, NPATCH2D + IF ((POP%POP_grid(g)%patch(P)%age(1) <= age(iage)) .AND.& + (POP%POP_grid(g)%patch(P)%age(1) > age_min)) THEN + age_min = POP%POP_grid(g)%patch(P)%age(1) + i_min = P + END IF + + IF ((POP%POP_grid(g)%patch(P)%age(1) >= age(iage)) .AND.& + (POP%POP_grid(g)%patch(P)%age(1) < age_max)) THEN + age_max = POP%POP_grid(g)%patch(P)%age(1) + i_max = P + END IF + END DO IteratePatches age_min_growth = age_min age_max_growth = age_max i_min_growth = i_min i_max_growth = i_max + ! Zero the fire mortality for this age + fire_mort_age(iage) = 0 + if ((i_min.gt.0).and.(i_max.gt.0).and.(age_max.eq.age_min)) then ! no need to interpolate + Npoints = 0.0 + IteratePatches_2: DO P = 1, NPATCH2D + IF (POP%POP_grid(g)%patch(P)%age(1) == age_min) THEN + + fire_mort_age(iage) = fire_mort_age(iage) +& + POP%POP_grid(g)%patch(P)%fire_mortality + + Npoints = Npoints + 1.0 + END IF + END DO IteratePatches_2 + + fire_mort_age(iage) = fire_mort_age(iage) / Npoints - MASK = pop%pop_grid(g)%patch(:)%age(1).eq.age_min - where (MASK) - tmp_array = 1.0_dp - elsewhere - tmp_array = 0.0_dp - endwhere - fire_mort_age(iage)= & - SUM(pop%pop_grid(g)%patch(:)%fire_mortality,MASK)/SUM(tmp_array) else ! interpolate or extrapolate if ((i_min.eq.0).and.(i_max.gt.0)) then @@ -2838,35 +3128,64 @@ SUBROUTINE INTERPOLATE_FIREMORTALITY(pop, disturbance_interval,it,g) ! extrapolate to higher age age_max = age_min i_max = i_min - age_min = MAXVAL(pop%pop_grid(g)%patch(:)%age(1), & - pop%pop_grid(g)%patch(:)%age(1).LT.age_max) - i_min = MAXLOC(pop%pop_grid(g)%patch(:)%age(1),1, & - pop%pop_grid(g)%patch(:)%age(1).LT.age_max) + IteratePatches_3: DO P = 1, NPATCH2D + IF (POP%POP_grid(g)%patch(P)%age(1) > age_min) THEN + age_min = POP%POP_grid(g)%patch(P)%age(1) + i_min = P + END IF + END DO IteratePatches_3 endif - if (i_min.ne.0.and.age_min.ne.0) then - MASK = pop%pop_grid(g)%patch(:)%age(1).eq.age_min - where (MASK) - tmp_array = 1.0_dp - elsewhere - tmp_array = 0.0_dp - endwhere - tmp5_min = SUM(pop%pop_grid(g)%patch(:)%fire_mortality,MASK)/SUM(tmp_array) + ! Do minimum values + Npoints = 0.0 + IteratePatches_4: DO P = 1, NPATCH2D + IF (POP%POP_grid(g)%patch(P)%age(1) == age_min) THEN + + fire_mort_min = fire_mort_min +& + POP%POP_grid(g)%patch(P)%fire_mortality + + Npoints = Npoints + 1.0 + END IF + END DO IteratePatches_4 + + fire_mort_min = fire_mort_min / Npoints + + !MASK = pop%pop_grid(g)%patch(:)%age(1).eq.age_min + !where (MASK) + !tmp_array = 1.0_dp + !elsewhere + !tmp_array = 0.0_dp + !endwhere + !tmp5_min = SUM(pop%pop_grid(g)%patch(:)%fire_mortality,MASK)/SUM(tmp_array) else - tmp5_min = 0.0_dp + fire_mort_min = 0.0_dp endif - MASK = pop%pop_grid(g)%patch(:)%age(1).eq.age_max - where (MASK) - tmp_array = 1.0_dp - elsewhere - tmp_array = 0.0_dp - endwhere - tmp5_max = SUM(pop%pop_grid(g)%patch(:)%fire_mortality,MASK)/SUM(tmp_array) + ! Do the maximum values + Npoints = 0.0 + IteratePatches_5: DO P = 1, NPATCH2D + IF (POP%POP_grid(g)%patch(P)%age(1) == age_max) THEN + + fire_mort_max = fire_mort_max +& + POP%POP_grid(g)%patch(P)%fire_mortality + + Npoints = Npoints + 1.0 + END IF + END DO IteratePatches_5 + + fire_mort_max = fire_mort_max / Npoints - fire_mort_age(iage) = tmp5_min + (tmp5_max-tmp5_min)/real(age_max-age_min,dp)* & - real(age(iage)-age_min,dp) + !MASK = pop%pop_grid(g)%patch(:)%age(1).eq.age_max + !where (MASK) + !tmp_array = 1.0_dp + !elsewhere + !tmp_array = 0.0_dp + !endwhere + !tmp5_max = SUM(pop%pop_grid(g)%patch(:)%fire_mortality,MASK)/SUM(tmp_array) + + fire_mort_age(iage) = fire_mort_min + (fire_mort_max-fire_mort_min) *& + REAL(age(iage) - age_min, dp) / REAL(age_max-age_min, dp) endif POP%pop_grid(g)%fire_mortality = POP%pop_grid(g)%fire_mortality + & @@ -2874,9 +3193,10 @@ SUBROUTINE INTERPOLATE_FIREMORTALITY(pop, disturbance_interval,it,g) enddo - DEALLOCATE(age) - DEALLOCATE(freq_age) - DEALLOCATE(fire_mort_age) + ! Shouldn't need these deallocates + !DEALLOCATE(age) + !DEALLOCATE(freq_age) + !DEALLOCATE(fire_mort_age) END SUBROUTINE INTERPOLATE_FIREMORTALITY @@ -2895,59 +3215,65 @@ SUBROUTINE ADJUST_POP_FOR_FIRE(pop,disturbance_interval, burned_area, FLI) INTEGER(i4b) :: g, np, c, k, it, nc REAL(dp) :: mort, cmass_stem, dbh + INTEGER :: P np = SIZE(POP%POP_grid) mort = 0.0 DO g=1,np - POP%pop_grid(g)%fire_mortality = 0.0_dp - - - - if (burned_area(g) > 0.0_dp) then + POP%pop_grid(g)%fire_mortality = 0.0_dp - it = maxval(pop%pop_grid(g)%patch(:)%age(1)) + 1 - DO k=1,NPATCH - nc = pop%pop_grid(g)%patch(k)%Layer(1)%ncohort + if (burned_area(g) > 0.0_dp) then + + it = 0 + IteratePatches: DO P = 1, NPATCH2D + IF (POP%POP_grid(g)%patch(P)%age(1) > it) THEN + it = POP%POP_grid(g)%patch(P)%age(1) + END IF + END DO IteratePatches - pop%pop_grid(g)%patch(k)%fire_mortality = 0.0_dp - DO c=1,nc - - dbh = pop%pop_grid(g)%patch(k)%layer(1)%cohort(c)%diameter*100.0_dp - cmass_stem = pop%pop_grid(g)%patch(k)%layer(1)%cohort(c)%biomass + it = it + 1 + + !it = maxval(pop%pop_grid(g)%patch(:)%age(1)) + 1 + DO k=1,NPATCH + nc = pop%pop_grid(g)%patch(k)%layer(1)%ncohort - mort = TopKill_Collins(dbh, FLI(g)) * burned_area(g) + pop%pop_grid(g)%patch(k)%fire_mortality = 0.0_dp + DO c=1,nc - pop%pop_grid(g)%patch(k)%fire_mortality = mort* & - pop%pop_grid(g)%patch(k)%Layer(1)%cohort(c)%biomass+ & - pop%pop_grid(g)%patch(k)%fire_mortality + dbh = pop%pop_grid(g)%patch(k)%layer(1)%cohort(c)%diameter*100.0_dp + cmass_stem = pop%pop_grid(g)%patch(k)%layer(1)%cohort(c)%biomass - pop%pop_grid(g)%patch(k)%Layer(1)%cohort(c)%biomass = cmass_stem*(1.0_dp-mort) - pop%pop_grid(g)%patch(k)%Layer(1)%cohort(c)%heartwood = & - pop%pop_grid(g)%patch(k)%Layer(1)%cohort(c)%heartwood *(1.0_dp-mort) - pop%pop_grid(g)%patch(k)%Layer(1)%cohort(c)%sapwood = & - pop%pop_grid(g)%patch(k)%Layer(1)%cohort(c)%sapwood *(1.0_dp-mort) + mort = TopKill_Collins(dbh, FLI(g)) * burned_area(g) - pop%pop_grid(g)%patch(k)%fire_top_kill_density = & - pop%pop_grid(g)%patch(k)%fire_top_kill_density + & - pop%pop_grid(g)%patch(k)%Layer(1)%cohort(c)%density *mort + pop%pop_grid(g)%patch(k)%fire_mortality = mort* & + pop%pop_grid(g)%patch(k)%layer(1)%cohort(c)%biomass+ & + pop%pop_grid(g)%patch(k)%fire_mortality - pop%pop_grid(g)%patch(k)%Layer(1)%cohort(c)%density = & - pop%pop_grid(g)%patch(k)%Layer(1)%cohort(c)%density*(1.0_dp-mort) + pop%pop_grid(g)%patch(k)%layer(1)%cohort(c)%biomass = cmass_stem*(1.0_dp-mort) + pop%pop_grid(g)%patch(k)%layer(1)%cohort(c)%heartwood = & + pop%pop_grid(g)%patch(k)%layer(1)%cohort(c)%heartwood *(1.0_dp-mort) + pop%pop_grid(g)%patch(k)%layer(1)%cohort(c)%sapwood = & + pop%pop_grid(g)%patch(k)%layer(1)%cohort(c)%sapwood *(1.0_dp-mort) + pop%pop_grid(g)%patch(k)%fire_top_kill_density = & + pop%pop_grid(g)%patch(k)%fire_top_kill_density + & + pop%pop_grid(g)%patch(k)%layer(1)%cohort(c)%density *mort - ENDDO + pop%pop_grid(g)%patch(k)%layer(1)%cohort(c)%density = & + pop%pop_grid(g)%patch(k)%layer(1)%cohort(c)%density*(1.0_dp-mort) + ENDDO - nc = pop%pop_grid(g)%patch(k)%Layer(1)%ncohort - pop%pop_grid(g)%patch(k)%biomass_old = pop%pop_grid(g)%patch(k)%Layer(1)%biomass - pop%pop_grid(g)%patch(k)%Layer(1)%biomass = & - SUM(pop%pop_grid(g)%patch(k)%Layer(1)%cohort(1:nc)%biomass) + nc = pop%pop_grid(g)%patch(k)%layer(1)%ncohort + pop%pop_grid(g)%patch(k)%biomass_old = pop%pop_grid(g)%patch(k)%layer(1)%biomass + pop%pop_grid(g)%patch(k)%layer(1)%biomass = & + SUM(pop%pop_grid(g)%patch(k)%layer(1)%cohort(1:nc)%biomass) - ! need to remove cohorts with very low density? - ! This will get done at the end of the year anyway + ! need to remove cohorts with very low density? + ! This will get done at the end of the year anyway - ENDDO + ENDDO ENDIF @@ -3404,13 +3730,14 @@ SUBROUTINE SMOOTH_FLUX(POP,g,t) TYPE(POP_TYPE), INTENT(INOUT) :: POP INTEGER(i4b), INTENT(IN) :: g, t - INTEGER(i4b), PARAMETER :: SPAN = NYEAR_WINDOW + INTEGER(i4b) :: SPAN REAL(dp) :: x(NYEAR_SMOOTH), y(NYEAR_SMOOTH), a, b, r REAL(dp) :: sumflux, sumsmooth, flux(NYEAR_HISTORY), smoothed_flux REAL(dp) :: dbuf INTEGER(i4b) :: t0, tt, n, k ! update fire_mortality_history + SPAN = NYEAR_WINDOW IF (t.gt.NYEAR_HISTORY) THEN DO k = 1, NYEAR_HISTORY-1 @@ -3456,13 +3783,14 @@ SUBROUTINE SMOOTH_FLUX_cat(POP,g,t) TYPE(POP_TYPE), INTENT(INOUT) :: POP INTEGER(i4b), INTENT(IN) :: g, t - INTEGER(i4b), PARAMETER :: SPAN = NYEAR_WINDOW + INTEGER(i4b) :: SPAN REAL(dp) :: x(NYEAR_SMOOTH), y(NYEAR_SMOOTH), a, b, r REAL(dp) :: sumflux, sumsmooth, flux(NYEAR_HISTORY), smoothed_flux REAL(dp) :: dbuf INTEGER(i4b) :: t0, tt, n, k ! update cat_mortality_history + SPAN = NYEAR_WINDOW IF (t.gt.NYEAR_HISTORY) THEN DO k = 1, NYEAR_HISTORY-1 @@ -3741,15 +4069,115 @@ SUBROUTINE alloc_POP(POP, arraysize) TYPE(POP_TYPE), INTENT(INOUT) :: POP INTEGER, INTENT(IN) :: arraysize + ! Iterators for allocation + INTEGER :: LS, P, L ! Iterate over landscapes, patches, layers + IF (.NOT. ALLOCATED(POP%POP_Grid)) ALLOCATE(POP%POP_Grid(arraysize)) IF (.NOT. ALLOCATED(POP%Iwood)) ALLOCATE(POP%Iwood(arraysize)) ! IF (.NOT. ALLOCATED(POP%LU)) ALLOCATE (POP%LU(arraysize)) IF (.NOT. ALLOCATED(POP%it_pop)) ALLOCATE(POP%it_pop(arraysize)) + ! We need to do some more work here now that the POP type is configurable + ! Are the allocated checks necessary? Why would they ever be already + ! allocated? + ! Need to loop over each subtype + AllocateLandscapes: DO LS = 1, arraysize + ALLOCATE(POP%POP_grid(LS)%patch(NPATCH2D)) + ALLOCATE(POP%POP_grid(LS)%freq(NPATCH2D)) + ALLOCATE(POP%POP_grid(LS)%freq_old(NPATCH2D)) + ALLOCATE(POP%POP_grid(LS)%fire_freq(NPATCH2D)) + ALLOCATE(POP%POP_grid(LS)%fire_freq_old(NPATCH2D)) + ALLOCATE(POP%POP_grid(LS)%cat_freq(NPATCH2D)) + ALLOCATE(POP%POP_grid(LS)%cat_freq_old(NPATCH2D)) + ALLOCATE(POP%POP_grid(LS)%freq_ranked_age_unique(NPATCH2D, NDISTURB)) + ALLOCATE(POP%POP_grid(LS)%ranked_age_unique(NPATCH2D, NDISTURB)) + ALLOCATE(POP%POP_grid(LS)%n_age(NDISTURB)) + ALLOCATE(POP%POP_grid(LS)%cmass_stem_bin(HEIGHT_BINS)) + ALLOCATE(POP%POP_grid(LS)%densindiv_bin(HEIGHT_BINS)) + ALLOCATE(POP%POP_grid(LS)%height_bin(HEIGHT_BINS)) + ALLOCATE(POP%POP_grid(LS)%diameter_bin(HEIGHT_BINS)) + ALLOCATE(POP%POP_grid(LS)%bin_labels(HEIGHT_BINS)) + ALLOCATE(POP%POP_grid(LS)%fire_mortality_history(NYEAR_HISTORY)) + ALLOCATE(POP%POP_grid(LS)%cat_mortality_history(NYEAR_HISTORY)) + ALLOCATE(POP%POP_grid(LS)%freq_age(AGEMAX)) + ALLOCATE(POP%POP_grid(LS)%biomass_age(AGEMAX)) + AllocatePatches: DO P = 1, NPATCH2D + ALLOCATE(POP%POP_grid(LS)%patch(P)%disturbance_interval(NDISTURB)) + ALLOCATE(POP%POP_grid(LS)%patch(P)%first_disturbance_year(NDISTURB)) + ALLOCATE(POP%POP_grid(LS)%patch(P)%age(NDISTURB)) + ALLOCATE(POP%POP_grid(LS)%patch(P)%layer(NLayer)) + AllocateLayers: DO L = 1, NLAYER + ALLOCATE(POP%POP_grid(LS)%patch(P)%layer(L)%Cohort(NCOHORT_MAX)) + END DO AllocateLayers + END DO AllocatePatches + END DO AllocateLandscapes + + END SUBROUTINE alloc_POP !******************************************************************************* + SUBROUTINE update_POP_parameters(disturbance_interval, disturbance_intensity) + !*## Purpose + ! + ! Update any POP parameters defined by the user in the namelist. + ! + !## Method + ! + ! Read the pop.nml namelist and use the values in said namelist to write + ! over the default values defined in the POP_Constants module. + + ! The POP_Constants contains all the values which parameterise the POP + ! algorithm. At the moment, all the values from said module are imported + ! into the POP module via USE, which we want to move away from in the + ! future. For now, given they're defined within the scope of this module, + ! we can just write over them with minimal effort. + ! There are two veg parameters that have come into play in the POP routines, + ! which are veg%disturbance_interval and veg%disturbance_intensity. Going to + ! import the veg_parameter_type we pass around to make those parameters + ! modifiable here, but might not be a good idea in the long term. + + INTEGER, DIMENSION(:,:), POINTER, INTENT(INOUT) :: disturbance_interval + REAL(dp), DIMENSION(:,:), POINTER, INTENT(INOUT) :: disturbance_intensity + + ! Status checker for namelist + INTEGER :: nmlUnit, ios + CHARACTER(LEN=200) :: ioMessage + + ! This is just in the order they appear in POP_Constants module + NAMELIST /popnml/ fulton_alpha, densindiv_max, densindiv_min, kbiometric,& + wd, growth_efficiency_min, pmort, mort_max,& + theta_recruit, cmass_stem_init, powerbiomass,& + powergrowthefficiency, crowdingfactor, alpha_cpc,& + k_allom1, k_rp, ksapwood, shootfrac, CtoNw, CtoNl, CtoNr,& + ncohort_max, ndisturb, patch_reps, patch_reps1,& + patch_reps2, height_bins, bin_power, nlayer,& + allom_switch, max_height_switch, resource_switch,& + interp_switch, smooth_switch, nyear_window,& + agemax, disturbance_interval, disturbance_intensity + + OPEN(NEWUNIT=nmlUnit, FILE='pop.nml', STATUS='OLD', ACTION='read',& + IOSTAT=ios) + + IF (ios == 29) THEN + ! Until this change is reflected in the common configuration, leave the + ! option to continue if the namelist doesn't exist. File not found returns + ! an IOSTAT=29. + CONTINUE + ELSE + ! In all other instances, handle the namelist reading as usual + READ(nmlUnit, NML=popnml, IOSTAT=ios, IOMSG=ioMessage) + !CALL handle_iostat(ios, ioMessage) + END IF + + ! Any parameters that are defined strictly in relation to other parameters + Npatch = patch_reps1 * patch_reps2 + Npatch2D = Npatch + Nyear_smooth = 2 * Nyear_window + 1 + Nyear_history = nyear_smooth - nyear_window + +END SUBROUTINE update_POP_parameters + END MODULE POPModule !******************************************************************************* diff --git a/core/biogeochem/POPLUC.F90 b/core/biogeochem/POPLUC.F90 index e59c88916..9614cb041 100644 --- a/core/biogeochem/POPLUC.F90 +++ b/core/biogeochem/POPLUC.F90 @@ -1651,6 +1651,8 @@ SUBROUTINE POPLUC_Init(POPLUC, LUC_EXPT, casapool, casaflux, casabiome, veg, POP ENDDO ENDIF + ! Set the disturbance interval from the veg type + disturbance_interval = veg%disturbance_interval END SUBROUTINE POPLUC_Init diff --git a/core/biogeochem/pop_io.F90 b/core/biogeochem/pop_io.F90 index 7b2d6bb4d..a3779180c 100755 --- a/core/biogeochem/pop_io.F90 +++ b/core/biogeochem/pop_io.F90 @@ -1329,13 +1329,13 @@ SUBROUTINE POP_IO(POP, casamet, YEAR, ACTION, CF) DO m = 1, mp SELECT CASE ( i ) CASE( 1) - POP%pop_grid(m)%biomass = R2(m,:) + POP%pop_grid(m)%biomass = R2(m,nlayer) CASE( 2) - POP%pop_grid(m)%density = R2(m,:) + POP%pop_grid(m)%density = R2(m,nlayer) CASE( 3) - POP%pop_grid(m)%hmean = R2(m,:) + POP%pop_grid(m)%hmean = R2(m,nlayer) CASE( 4) - POP%pop_grid(m)%hmax = R2(m,:) + POP%pop_grid(m)%hmax = R2(m,nlayer) CASE default write(*,*) "Parameter not assigned in pop_io.F90!" #ifdef __MPI__ diff --git a/offline/cable_input.F90 b/offline/cable_input.F90 index 6ec79fec0..8b69ac130 100644 --- a/offline/cable_input.F90 +++ b/offline/cable_input.F90 @@ -2431,7 +2431,7 @@ SUBROUTINE load_parameters(met, air, ssnow, veg, bgc, soil, canopy, rough, rad, ! landpt%type - via cable_IO_vars_module (nap,cstart,cend,ilon,ilat) ! max_vegpatches - via cable_IO_vars_module !! vh_js !! - USE POPmodule, ONLY: POP_INIT + USE POPmodule, ONLY: POP_INIT, adjust_POP_parameters USE POPLUC_module, ONLY: POPLUC_INIT USE CABLE_LUC_EXPT, ONLY: LUC_EXPT_TYPE use casaparm, only: initcasa @@ -2580,6 +2580,17 @@ SUBROUTINE load_parameters(met, air, ssnow, veg, bgc, soil, canopy, rough, rad, ENDIF ENDDO + ! The parameters veg%disturbance_interval and veg%disturbance_intensity + ! are intended to part of the tunable parameters for POP. We don't want + ! to expose the veg type to the POP routines directly, but still have + ! them modifiable via namelist, so we call the namelist reading + ! routine here. Unfortunately, all the POP parameters exist in a POP + ! module, so we're going to do a bit of modifying module data and a bit + ! of modifying argument data. This is far from a perfect solution, but + ! an acceptable temporary one. + CALL update_POP_parameters(veg%disturbance_interval,& + veg%disturbance_intensity) + CALL POP_init(POP, veg%disturbance_interval(Iwood,:), mp_POP, Iwood) IF (.NOT. (spinup .OR. CABLE_USER%POP_fromZero)) & CALL POP_IO(POP, casamet, cable_user%YearStart, "READ_RST", .TRUE.)