! This file contains Fortran90 code for the POP model, ! a stand-alone tree demography and landscape structure module for Earth system models ! 17-01-2014 ! Written by Vanessa Haverd, Ben Smith and Lars Nieradzik ! Report Bugs to Vanessa.Haverd@csiro.au !CITATION !-------------------------------------------------------- !When referring to this code in publications, please cite: ! Haverd, V., Smith, B., Cook, G., Briggs, P.R., Nieradzik, L., Roxburgh, S.R., Liedloff, A., ! Meyer, C.P. and Canadell, J.G., 2013. ! A stand-alone tree demography and landscape structure module for Earth system models. ! Geophysical Research Letters, 40: 1-6. !DISCLAIMER, COPYRIGHT AND LICENCE !-------------------------------------------------------- ! Use of this code is subject to the Legal Notice and Disclaimer at ! http://www.csiro.au/org/LegalNoticeAndDisclaimer.html ! This code is Copyright, CSIRO, 2014. ! This code is made available under the conditions of the Creative Commons ! Attribution-Share Alike 3.0 License: ! http://creativecommons.org/licenses/by-sa/3.0/ !=============================================================================== MODULE POPModule !------------------------------------------------------------------------------- ! * This module contains all subroutines for POP calcs at a single time step. !------------------------------------------------------------------------------- USE TYPEdef, ONLY: sp, i4b USE POP_Types USE POP_Constants CONTAINS !============================================================================= SUBROUTINE ZeroPOP(POP,n) TYPE(POP_TYPE), INTENT(INOUT) :: POP INTEGER, OPTIONAL, INTENT(IN) ::n INTEGER:: g,k,l,c, np,a,b IF ( .NOT. ALLOCATED(pop%pop_grid) ) THEN WRITE(*,*)" POP not allocated! Abort in ZeroPOP." STOP -1 ENDIF np = SIZE(pop%pop_grid) ! optional interger n intended for zeroing secondary forest tiles IF (PRESENT(n)) THEN a = n b = n !pop%LU(n) = 2 POP%pop_grid(n)%LU = 2 ELSE a = 1 b = np !pop%LU = 1 POP%pop_grid%LU = 1 ENDIF DO g=a,b POP%pop_grid(g)%freq = 0 ! patch weighting POP%pop_grid(g)%freq_old = 0 ! patch weighting POP%pop_grid(g)%fire_freq = 0 POP%pop_grid(g)%fire_freq_old = 0 POP%pop_grid(g)%cat_freq = 0 POP%pop_grid(g)%cat_freq_old = 0 POP%pop_grid(g)%biomass = 0 ! landscape stem biomass (weighted mean over patches) POP%pop_grid(g)%density= 0 ! landscape tree density (weighted mean over patches) POP%pop_grid(g)%hmean= 0 ! landscape mean treen height (weighted mean over patches) POP%pop_grid(g)%hmax= 0 ! landscape max tree height POP%pop_grid(g)%cmass_stem_bin= 0 ! biomass by height bin POP%pop_grid(g)%densindiv_bin= 0 ! density by height bin POP%pop_grid(g)%height_bin= 0 ! mean height in each bin POP%pop_grid(g)%diameter_bin= 0 ! mean diameter in each bin POP%pop_grid(g)%bin_labels= ' ' ! text strings for bin bounds POP%pop_grid(g)%cmass_sum= 0 ! landscape biomass POP%pop_grid(g)%csapwood_sum= 0 ! landscape biomass POP%pop_grid(g)%cmass_sum_old= 0 ! landscape biomass POP%pop_grid(g)%csapwood_sum_old= 0 ! landscape biomass POP%pop_grid(g)%cheartwood_sum= 0 ! landscape biomass POP%pop_grid(g)%densindiv= 0 ! landscape density of individuals POP%pop_grid(g)%height_mean= 0 POP%pop_grid(g)%height_max= 0 POP%pop_grid(g)%basal_area= 0 POP%pop_grid(g)%sapwood_loss = 0 POP%pop_grid(g)%sapwood_area_loss = 0 POP%pop_grid(g)%crowding_mortality = 0 ! (kg C m-2 y-1) POP%pop_grid(g)%stress_mortality = 0 ! (kg C m-2 y-1) POP%pop_grid(g)%fire_mortality = 0 ! (kg C m-2 y-1) POP%pop_grid(g)%cat_mortality = 0 ! (kg C m-2 y-1) POP%pop_grid(g)%res_mortality = 0 ! (kg C m-2 y-1) POP%pop_grid(g)%growth= 0 POP%pop_grid(g)%area_growth= 0 POP%pop_grid(g)%crown_cover = 0 POP%pop_grid(g)%crown_area = 0 POP%pop_grid(g)%crown_volume = 0 POP%pop_grid(g)%sapwood_area = 0 POP%pop_grid(g)%sapwood_area_old = 0 POP%pop_grid(g)%Kclump = 1 POP%pop_grid(g)%fire_mortality_smoothed = 0 POP%pop_grid(g)%cat_mortality_smoothed = 0 POP%pop_grid(g)%smoothing_buffer = 0 POP%pop_grid(g)%smoothing_buffer_cat = 0 POP%pop_grid(g)%fire_mortality_history = 0 POP%pop_grid(g)%cat_mortality_history = 0 POP%pop_grid(g)%freq_age = 0 IF (PRESENT(n)) THEN POP%pop_grid(g)%freq_age(1) = 1 ENDIF POP%pop_grid(g)%biomass_age = 0 DO k=1,NPATCH2D POP%pop_grid(g)%patch(k)%factor_recruit= 0 POP%pop_grid(g)%patch(k)%pgap= 0 POP%pop_grid(g)%patch(k)%lai= 0 POP%pop_grid(g)%patch(k)%biomass= 0 ! total biomass in patch POP%pop_grid(g)%patch(k)%biomass_old= 0 POP%pop_grid(g)%patch(k)%sapwood = 0 ! total biomass in patch (sapwood) POP%pop_grid(g)%patch(k)%heartwood = 0 ! total biomass in patch (heartwood) POP%pop_grid(g)%patch(k)%sapwood_old= 0 POP%pop_grid(g)%patch(k)%sapwood_area = 0 POP%pop_grid(g)%patch(k)%sapwood_area_old= 0 POP%pop_grid(g)%patch(k)%stress_mortality = 0 ! biomass lost in each patch due to stress POP%pop_grid(g)%patch(k)%fire_mortality= 0 ! biomass lost in each patch due to fire partial dist POP%pop_grid(g)%patch(k)%cat_mortality= 0 ! biomass lost in each patch due to fire partial dist POP%pop_grid(g)%patch(k)%crowding_mortality = 0 POP%pop_grid(g)%patch(k)%cpc = 0 POP%pop_grid(g)%patch(k)%mortality = 0 POP%pop_grid(g)%patch(k)%sapwood_loss = 0 POP%pop_grid(g)%patch(k)%sapwood_area_loss = 0 POP%pop_grid(g)%patch(k)%growth= 0 ! biomass growth in each patch due stem increment POP%pop_grid(g)%patch(k)%area_growth= 0 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 POP%pop_grid(g)%patch(k)%frac_respiration = 0 POP%pop_grid(g)%patch(k)%frac_light_uptake = 0 DO l=1,NLAYER POP%pop_grid(g)%patch(k)%Layer(L)%ncohort = 0 ! number of cohorts with density >0 POP%pop_grid(g)%patch(k)%Layer(L)%biomass = 0 ! layer biomass POP%pop_grid(g)%patch(k)%Layer(L)%density= 0 ! layer tree density POP%pop_grid(g)%patch(k)%Layer(L)%hmean= 0 ! layer mean tree height (weighted mean over patches) POP%pop_grid(g)%patch(k)%Layer(L)%hmax= 0 ! 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 ! cohort biomass POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%density = 0 ! landscape tree density (weighted mean over patches) POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%frac_resource_uptake = 0 POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%frac_light_uptake = 0 POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%frac_interception = 0 POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%frac_respiration = 0 POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%frac_NPP = 0 POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%respiration_scalar = 0 POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%crown_area = 0 POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%Pgap = 0 POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%height = 0 POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%diameter = 0 POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%sapwood = 0 ! cohort sapwood POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%heartwood = 0 ! cohort heartwood POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%sapwood_area = 0 POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%basal_area = 0 POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%LAI = 0 POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%Cleaf = 0 POP%pop_grid(g)%patch(k)%Layer(L)%cohort(c)%Croot = 0 ENDDO ENDDO ENDDO ENDDO 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 IMPLICIT NONE TYPE(POP_TYPE), INTENT(INOUT) :: POP INTEGER(i4b), INTENT(IN) :: mean_disturbance_interval(:,:) INTEGER(i4b), INTENT(IN), OPTIONAL :: m INTEGER(i4b) :: j, k, g, ipatch, idist, p, c, n, i INTEGER(i4b) :: disturbance_interval INTEGER(i4b):: patch_disturbance_interval_idist(NDISTURB,NPATCH2D) INTEGER(i4b):: patch_first_disturbance_year_idist(NDISTURB,NPATCH2D), & patch_first_disturbance_year_unique(NDISTURB,NPATCH2D) INTEGER(i4b):: Poisson_age(1000),Poisson_freq(1000) REAL(dp):: Poisson_weight(1000), CumPoisson_weight(1000) INTEGER(i4b):: disturbances_per_timebase, timebase INTEGER:: i_min, i_max, age_sample(2,NAGE_MAX**2), tmp(NAGE_MAX**2) INTEGER:: age_tmp, tmp_unique(NAGE_MAX**2), n_age, np REAL(dp):: disturbance_freq, tmp1 INTEGER:: tmp2(PATCH_REPS1) , n_first_disturbance_year(NDISTURB), tmp3(PATCH_REPS2) INTEGER:: a,b np = SIZE(POP%pop_grid) a = 1 b = np IF (PRESENT(m)) THEN a = m b = m ENDIF DO g=a,b ! calculate Poisson weights for each of the 2 mean disturbance intervals IF (NPATCH.GT.1) THEN DO idist=1,NDISTURB disturbance_freq=1.0/REAL(mean_disturbance_interval(g,idist)) DO p = 1,1000 Poisson_age(p) = p Poisson_weight(p) = Exponential(disturbance_freq,p) CumPoisson_weight(p) = CumExponential(disturbance_freq,REAL(p,dp)) ENDDO ! set max age to correspond to 95% percentile of cum pdf DO k =1,NPATCH2D i_max = MAXLOC(Poisson_age,1,CumPoisson_weight.LE.0.95) POP%pop_grid(g)%patch(k)%disturbance_interval(idist) = Poisson_age(i_max) POP%pop_grid(g)%patch(k)%id = k POP%pop_grid(g)%patch(k)%age = 0 ENDDO ENDDO ! DO idist =1,ndisturb ! set first disturbance year for first dist interval class IF (idist .EQ. 1) THEN disturbance_interval = POP%pop_grid(g)%patch(1)%disturbance_interval(idist) DO c = 1,PATCH_REPS1 IF (c==1) THEN tmp2(1) = 1 ELSE tmp2(1) = MAX(disturbance_interval*(c-1)/(PATCH_REPS1),1)+1 ENDIF tmp2(2) = MAX(disturbance_interval*c/(PATCH_REPS1),1) tmp2(3) = tmp2(1) ! write(*,*) 'tmp2', c, disturbance_interval, tmp2(1),tmp2(2) DO j = 1,PATCH_REPS2 ipatch = (c-1)*PATCH_REPS2 + j POP%pop_grid(g)%patch(ipatch)%first_disturbance_year(idist) = tmp2(3) tmp2(3)=tmp2(3)+1 IF (tmp2(3)>tmp2(2)) THEN tmp2(3) = tmp2(1) ENDIF ENDDO ENDDO !$ ! set first disturbance year for first dist interval class !$ idist = 1 !$ disturbance_interval = POP%pop_grid(g)%patch(1)%disturbance_interval(idist) !$ DO c = 1,PATCH_REPS1 !$ tmp2(c) = max(disturbance_interval*c/(PATCH_REPS1),1) !$ ENDDO !$ DO c = 1,PATCH_REPS1 !$ i = 0 !$ DO j = 1,PATCH_REPS2 !$ ipatch = (j-1)*PATCH_REPS1 + c !$ i = i+1 !$ IF (i.gt.PATCH_REPS1) then !$ i = 1 !$ ENDIF !$ do while ((tmp2(i+1).eq. tmp2(i)).and.(i.lt.PATCH_REPS1)) !$ i = i+1 !$ IF (i.gt.PATCH_REPS1) then !$ i = 1 !$ ENDIF !$ !$ ENDDO !$ !$ !$ write(*,*) i, tmp2(i) !$ POP%pop_grid(g)%patch(ipatch)%first_disturbance_year(idist) = tmp2(i) !$ !$ ENDDO !$ ENDDO ! set first disturbance year for first 2nd interval class ELSEIF (idist.EQ.2) THEN disturbance_interval = POP%pop_grid(g)%patch(1)%disturbance_interval(idist) DO c = 1,PATCH_REPS2 tmp3(c) = MAX(disturbance_interval*(c-1)/(PATCH_REPS2),1) ENDDO DO c = 1,PATCH_REPS2 i = 0 DO j = 1,PATCH_REPS1 ipatch = (j-1)*PATCH_REPS2 + c POP%pop_grid(g)%patch(ipatch)%first_disturbance_year(idist) = & tmp3(c) +(j-1)*MAX((tmp3(idist)-tmp3(1))/PATCH_REPS1,1) ! i = i+1 ! if (i.gt.(tmp3(2)-tmp3(1))) i = 0 ENDDO ENDDO ENDIF ENDDO ELSE ! NPATCH =1 (single patch mode) k = 1 DO idist=1,NDISTURB POP%pop_grid(g)%patch(k)%disturbance_interval(idist) = mean_disturbance_interval(g,idist) POP%pop_grid(g)%patch(k)%first_disturbance_year(idist) = 113 POP%pop_grid(g)%patch(k)%age = 0 POP%pop_grid(g)%patch(k)%id = k ENDDO ENDIF POP%pop_grid(g)%npatch_active = NPATCH ENDDO END SUBROUTINE InitPOP2D_Poisson ! !============================================================================= SUBROUTINE POPStep(POP, StemNPP, disturbance_interval, disturbance_intensity,LAI,Cleaf,Croot, & NPPtoGPP, StemNPP_av,frac_intensity1,precip ) IMPLICIT NONE TYPE(POP_TYPE), INTENT(INOUT) :: POP REAL(dp), INTENT(IN) :: StemNPP(:,:) REAL(dp), INTENT(IN) :: disturbance_intensity(:,:) INTEGER(i4b), INTENT(IN) :: disturbance_interval(:,:) REAL(dp), INTENT(IN) :: LAI(:) REAL(dp), INTENT(IN) :: Cleaf(:) REAL(dp), INTENT(IN) :: Croot(:) REAL(dp), INTENT(IN) :: NPPtoGPP(:) REAL(dp), INTENT(IN), OPTIONAL :: frac_intensity1(:), precip(:) REAL(dp), INTENT(IN), OPTIONAL :: StemNPP_av(:) INTEGER(i4b) :: idisturb,np,g INTEGER(i4b), ALLOCATABLE :: it(:) REAL(dp):: dallocW pop%it_pop = pop%it_pop + 1 !it = pop%it_pop(1) np = SIZE(POP%POP_grid) ALLOCATE(it(np)) 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, it) !$ ENDDO !$ !$ CALL GetPatchFrequencies(POP,it) IF (PRESENT(precip)) THEN IF(PRESENT(StemNPP_av)) THEN CALL PatchAnnualDynamics(POP, StemNPP,NPPtoGPP,disturbance_interval, it, precip=precip,StemNPP_av=StemNPP_av) ELSE CALL PatchAnnualDynamics(POP, StemNPP,NPPtoGPP,disturbance_interval, it, precip=precip) ENDIF ELSE IF(PRESENT(StemNPP_av)) THEN CALL PatchAnnualDynamics(POP, StemNPP,NPPtoGPP,disturbance_interval, it,StemNPP_av=StemNPP_av) ELSE CALL PatchAnnualDynamics(POP, StemNPP,NPPtoGPP,disturbance_interval, it) ENDIF ENDIF IF (NDISTURB.EQ.1) THEN IF (PRESENT(precip)) THEN ! CALL Patch_disturb(POP,it,1,precip) CALL Patch_partial_disturb2(POP,1,precip) ELSE CALL Patch_disturb(POP,1) ! CALL Patch_partial_disturb2(POP,it,1) ENDIF ELSEIF (NDISTURB.EQ.2) THEN IF (PRESENT(frac_intensity1)) THEN IF (PRESENT(precip)) THEN CALL Patch_partial_disturb(POP,1,disturbance_intensity,precip,frac_intensity1=frac_intensity1) ELSE CALL Patch_partial_disturb(POP,1,disturbance_intensity,frac_intensity1=frac_intensity1) ENDIF ELSE IF (PRESENT(precip)) THEN CALL Patch_partial_disturb(POP,1,disturbance_intensity,precip=precip) ELSE CALL Patch_partial_disturb(POP,1,disturbance_intensity) ENDIF ENDIF IF (PRESENT(precip)) THEN !CALL Patch_partial_disturb2(POP,it,2,precip) CALL Patch_disturb(POP,2,precip) ELSE ! CALL Patch_partial_disturb2(POP,it,2) CALL Patch_disturb(POP,2) ENDIF ENDIF DO idisturb = 1,NDISTURB CALL GetUniqueAgeFrequencies(POP, disturbance_interval, idisturb, it) ENDDO CALL GetPatchFrequencies(POP,it) !PRINT*,"Get Diags" IF (PRESENT(precip)) THEN CALL GetDiagnostics(pop, LAI,Cleaf,Croot, disturbance_interval, it,precip) ELSE CALL GetDiagnostics(pop, LAI,Cleaf,Croot, disturbance_interval, it) ENDIF END SUBROUTINE POPStep !============================================================================= SUBROUTINE PatchAnnualDynamics(pop, StemNPP,NPPtoGPP, disturbance_interval, it, precip,StemNPP_av) IMPLICIT NONE TYPE( POP_TYPE ), INTENT(INOUT) :: pop REAL(dp), INTENT(IN) :: StemNPP(:,:) REAL(dp), INTENT(IN) :: NPPtoGPP(:) INTEGER(i4b), INTENT(IN) :: disturbance_interval(:,:) REAL(dp), INTENT(IN), OPTIONAL :: precip(:) REAL(dp), OPTIONAL, INTENT(IN) :: StemNPP_av(:) INTEGER(i4b), INTENT(IN) :: it(:) REAL(dp) :: f, mu, densindiv, cmass REAL(dp) :: tmp,tmp_light,tmp_respiration,tmp_fracnpp, cmass_stem_sum,cmass_stem_inc INTEGER(i4b) :: j, k,c, ncohort, idist INTEGER(i4b) :: ivec(NCOHORT_MAX), nc, tmp1(NPATCH2D), np, idisturb REAL(dp) :: growth_efficiency,cmass_stem REAL(dp) :: mort, mort_bg, fire_mort REAL(dp) :: s2, cpc, crown_area REAL(dp) :: mort_cpc REAL(dp) :: basal, ht, diam, area_growth_grid , basal_grid, basal_new, basal_old REAL(dp) :: tmp2(NCOHORT_MAX), freq,tmp3(NPATCH2D),tmp4(NPATCH2D),tmp5(NPATCH2D) idisturb = 1 np = SIZE(POP%POP_grid) ! growth ! Distributes layer biomass increment among cohorts and increments age ! calculate fractional resource uptake by each cohort DO j=1,np basal_grid = 0.0 area_growth_grid = 0.0 pop%pop_grid(j)%sapwood_area_old = pop%pop_grid(j)%sapwood_area pop%pop_grid(j)%freq_old = pop%pop_grid(j)%freq pop%pop_grid(j)%fire_freq_old = pop%pop_grid(j)%fire_freq pop%pop_grid(j)%cat_freq_old = pop%pop_grid(j)%cat_freq ! Get fraction allocation for each patch tmp = 0.0 tmp_light = 0.0 tmp_respiration = 0.0 tmp_fracNPP = 0.0 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 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_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 ENDDO tmp_light = tmp_light + freq*(1. - pop%pop_grid(j)%patch(k)%Pgap) ENDDO IF (tmp_respiration .GT. 1.e-8 .AND. tmp_light .GT. 1.e-8) THEN DO k=1,NPATCH2D ! fraction respiration and un-normalised NPP for each patch 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 ! frac gpp pop%pop_grid(j)%patch(k)%frac_light_uptake = & (1. - pop%pop_grid(j)%patch(k)%pgap) /tmp_light ! frac npp pop%pop_grid(j)%patch(k)%frac_NPP = & MAX(pop%pop_grid(j)%patch(k)%frac_light_uptake*(1/NPPtoGPP(j)) - & pop%pop_grid(j)%patch(k)%frac_respiration*(1/NPPtoGPP(j)-1.0),0.0_dp) tmp_fracNPP = tmp_fracNPP + freq*pop%pop_grid(j)%patch(k)%frac_NPP ENDDO ! normalised fraction NPP DO k=1,NPATCH2D pop%pop_grid(j)%patch(k)%frac_NPP = & pop%pop_grid(j)%patch(k)%frac_NPP/tmp_fracNPP ENDDO ELSE pop%pop_grid(j)%patch(:)%frac_NPP = 1.0 pop%pop_grid(j)%patch(:)%frac_respiration = 1.0 pop%pop_grid(j)%patch(:)%frac_light_uptake = 1.0 ENDIF ELSE pop%pop_grid(j)%patch(:)%frac_NPP = 1.0 pop%pop_grid(j)%patch(:)%frac_respiration = 1.0 pop%pop_grid(j)%patch(:)%frac_light_uptake = 1.0 ENDIF ! End Get fraction allocation for each patch ! Get fraction allocation for each cohort in each patch DO k=1,NPATCH2D tmp = 0.0 tmp_light = 0.0 tmp_respiration = 0.0 tmp_fracNPP = 0.0 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 densindiv = pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%density ! get initial basal area IF ( PRESENT(precip) ) THEN CALL GET_ALLOMETRY( ALLOM_SWITCH, cmass_stem, densindiv, ht, diam, basal_old, precip(j)) ELSE CALL GET_ALLOMETRY( ALLOM_SWITCH, cmass_stem, densindiv, ht, diam, basal_old ) ENDIF IF (ALLOM_SWITCH.EQ.1) THEN ! assumes crown radius (m) = 0.1492 * dbh (cm) (from G. Cook, pers. comm.) crown_area = densindiv*PI*(diam*100.*0.1492)**2 ELSE crown_area = densindiv*PI*(((k_allom1 * diam ** k_rp )/PI)**0.5)**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_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 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 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 ELSE pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_interception = 1.0 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 !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) 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 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/NPPtoGPP(j)) - & pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_respiration*(1/NPPtoGPP(j)-1.0),0.0_dp) 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 ENDDO ! fraction net resource uptake 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 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)%frac_NPP ENDIF ENDDO ELSE c = 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 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 ENDIF ENDIF ENDDO tmp = 0 DO k=1,NPATCH2D pop%pop_grid(j)%patch(k)%sapwood_loss = 0.0 pop%pop_grid(j)%patch(k)%sapwood_area_loss = 0.0 pop%pop_grid(j)%patch(k)%sapwood_old = pop%pop_grid(j)%patch(k)%sapwood pop%pop_grid(j)%patch(k)%sapwood_area_old = pop%pop_grid(j)%patch(k)%sapwood_area pop%pop_grid(j)%patch(k)%biomass_old = pop%pop_grid(j)%patch(k)%biomass pop%pop_grid(j)%patch(k)%growth = 0.0 pop%pop_grid(j)%patch(k)%area_growth = 0.0 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 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 IF ( PRESENT(precip) ) THEN CALL GET_ALLOMETRY( ALLOM_SWITCH, cmass_stem, densindiv, ht, diam, basal_old, precip(j)) ELSE CALL GET_ALLOMETRY( ALLOM_SWITCH, cmass_stem, densindiv, ht, diam, basal_old ) 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) 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 CALL GET_ALLOMETRY( ALLOM_SWITCH, cmass_stem, densindiv, ht, diam, basal_new, precip(j)) ELSE CALL GET_ALLOMETRY( ALLOM_SWITCH, cmass_stem, densindiv, ht, diam, basal_new ) ENDIF ! 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) ! 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 ! 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 ! decrease sapwood in cohort (accounting for loss to heartwood) pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood = & (1. - 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 ! 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) ! patch sapwood area increment pop%pop_grid(j)%patch(k)%area_growth = pop%pop_grid(j)%patch(k)%area_growth + & basal_new-basal_old ! 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 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) ENDDO ENDDO ! Mortality !Implements resource stress mortality and crowding mortality for all cohorts in layer DO j=1,np DO k=1,NPATCH2D nc = 0 ivec = 0 pop%pop_grid(j)%patch(k)%stress_mortality = 0.0 pop%pop_grid(j)%patch(k)%fire_mortality = 0.0 pop%pop_grid(j)%patch(k)%crowding_mortality = 0.0 pop%pop_grid(j)%patch(k)%mortality = 0.0 crown_area = 0.0 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_av)) THEN growth_efficiency=pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%frac_resource_uptake* & StemNPP_av(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)) mort=MORT_MAX/(1.0+(growth_efficiency/GROWTH_EFFICIENCY_MIN)**Pmort) ! mort = 0 ! test pop%pop_grid(j)%patch(k)%stress_mortality = pop%pop_grid(j)%patch(k)%stress_mortality & + mort*cmass_stem IF (pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%diameter*100.GT.1.) 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* & PI*(pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%diameter*100.*0.14)**2 ELSE 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 crown_area = crown_area + 0.5*pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%LAI ENDIF cpc = 1. - EXP(-crown_area) pop%pop_grid(j)%patch(k)%cpc = cpc IF (cpc.GT.1e-3 .AND. alpha_cpc * (1. - 1./cpc).GT.-50.0) THEN mort_cpc = EXP(alpha_cpc * (1. - 1./cpc)) ELSE mort_cpc = 0. ENDIF !mort_cpc = 0 ! test pop%pop_grid(j)%patch(k)%crowding_mortality = pop%pop_grid(j)%patch(k)%crowding_mortality + & MIN((mort_cpc*CrowdingFactor),cmass_stem_inc/cmass_stem)*cmass_stem mort = mort + MIN((mort_cpc*CrowdingFactor),cmass_stem_inc/cmass_stem) pop%pop_grid(j)%patch(k)%mortality = pop%pop_grid(j)%patch(k)%mortality + 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 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 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%biomass = cmass_stem*(1.-mort) pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%heartwood = & pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%heartwood *(1.-mort) pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood = & pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood *(1.-mort) pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%density = & pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%density*(1.-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 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)%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)%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 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%heartwood = 0.0 ELSE 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 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 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%sapwood = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%sapwood_area = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%basal_area = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%heartwood = 0.0 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) ENDDO ENDDO ! recruitment IF (PRESENT(precip)) THEN CALL layer_recruitment(pop, precip) ELSE CALL layer_recruitment(pop) ENDIF ! Update time since last patch disturbance DO j=1,np DO k=1,NPATCH2D DO idist =1, NDISTURB pop%pop_grid(j)%patch(k)%age(idist) = pop%pop_grid(j)%patch(k)%age(idist) + 1 ENDDO ENDDO ENDDO END SUBROUTINE PatchAnnualDynamics !============================================================================= SUBROUTINE GetUniqueAgeFrequencies(pop, disturbance_interval, idisturb, it) IMPLICIT NONE TYPE(POP_TYPE), INTENT(INOUT) :: POP INTEGER(i4b), INTENT(IN) :: disturbance_interval(:,:), idisturb, it(:) INTEGER(i4b) :: g, i,j,k,ct,lastct,agecopy,idcopy REAL(dp), ALLOCATABLE :: midpoint(:) INTEGER(i4b), ALLOCATABLE :: ranked_age(:), ranked_age_init(:) INTEGER(i4b) :: tmp_count, tmp_i, age_tmp INTEGER(i4b), ALLOCATABLE :: ranked_age_unique_id(:), ranked_age_id(:), counter(:) REAL(dp), ALLOCATABLE :: tmp(:), freq_tmp(:), freq_tmp1(:) REAL(dp) :: p,cump,lastcump, freq, tmp1 INTEGER(i4b) :: n_age ! number of unique ages INTEGER(i4b) :: npatch_active ! number of active patches REAL(dp):: disturbance_freq INTEGER(i4b) :: i_max, age_max, Poisson_age(1000), np REAL(dp):: Poisson_weight(1000), CumPoisson_weight(1000) INTEGER(i4b), ALLOCATABLE :: bound(:,:), unique_age(:) !Fills array freq with weights (frequencies across landscape) for each unique age ! given specified mean disturbance interval 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)) ! 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 ranked_age_unique_id = 0 freq_tmp = 0.0 freq = 0 pop%pop_grid(g)%freq_ranked_age_unique(:, idisturb) = 0.0 midpoint = 0.0 DO i = 1, npatch_active -1 DO j = i+1, npatch_active IF (ranked_age(i).GT.ranked_age(j)) THEN agecopy = ranked_age(i) idcopy = ranked_age_id(i) ranked_age(i) = ranked_age(j) ranked_age_id(i) = ranked_age_id(j) ranked_age(j) = agecopy ranked_age_id(j) = idcopy ENDIF ENDDO ENDDO ! subset to unique ages k=0 age_tmp = -1 DO i = 1, npatch_active 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) age_tmp = ranked_age(i) n_age = k ENDDO disturbance_freq=1.0/REAL(disturbance_interval(g,idisturb)) DO i =1,1000 Poisson_age(i) = i CumPoisson_weight(i) = CumExponential(disturbance_freq,REAL(i,dp)) ENDDO ! construct upper and lower bounds for each unique age: these set the range of ages to be ! represented by an unique age ALLOCATE(bound(n_age,2)) ALLOCATE (unique_age(n_age)) bound = 0 unique_age = pop%pop_grid(g)%ranked_age_unique(1:n_age,idisturb) DO i=1,n_age IF (unique_age(i).EQ.0) THEN bound(i,1) = 0 bound(i,2) = 0 ELSEIF ((i.EQ.1).AND.(unique_age(i).GT.0)) THEN bound(i,1) = 0 bound(i,2) = unique_age(i) ELSEIF ((unique_age(i).GT.0).AND.(i.GT.1).AND.(unique_age(i-1).EQ.unique_age(i)-1)) THEN bound(i,1) = unique_age(i) IF (i.LT.n_age) THEN bound(i,2) = unique_age(i) ELSE i_max = MAXLOC(Poisson_age, 1, CumPoisson_weight.LE.0.99) bound(i, 2) = Poisson_age(i_max) ENDIF ELSEIF ((unique_age(i).GT.0).AND.(i.GT.1).AND.(unique_age(i-1).NE.unique_age(i)-1)) THEN bound(i,1) = bound(i-1,2)+1 IF (i.LT.n_age) THEN bound(i,2) = (unique_age(i)+ unique_age(i+1))/2 ELSE i_max = MAXLOC(Poisson_age, 1, CumPoisson_weight.LE.0.99) bound(i, 2) = Poisson_age(i_max) ENDIF ENDIF ENDDO ! calculate weighting for each unique age DO i=1,n_age DO j = bound(i,1),bound(i,2) !IF (pop%LU(g)==2) THEN ! secondary forest IF (POP%pop_grid(g)%LU ==2) THEN freq_tmp(i) = freq_tmp(i) + pop%pop_grid(g)%freq_age(j+1) ELSE freq_tmp(i) = freq_tmp(i) + REALExponential(disturbance_freq,REAL(j,dp)) ENDIF ENDDO ENDDO pop%pop_grid(g)%freq_ranked_age_unique(1:npatch_active,idisturb) = freq_tmp pop%pop_grid(g)%n_age(idisturb) = n_age DEALLOCATE (bound) DEALLOCATE (unique_age) ENDDO END SUBROUTINE GetUniqueAgeFrequencies !============================================================================= SUBROUTINE GetPatchFrequencies(pop,it) IMPLICIT NONE TYPE(POP_TYPE), INTENT(INOUT) :: POP INTEGER(i4b), INTENT(IN) :: it(:) INTEGER(i4b) :: n1, n2, g, REPCOUNT, tmp1(NPATCH1D), np, idist REAL(dp) :: tmp2(NPATCH1D), tmp3(NPATCH1D), sum_freq np = SIZE(Pop%pop_grid) DO g=1,np pop%pop_grid(g)%freq = 0.0 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 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 ENDDO ENDDO ENDIF ENDDO ! end loop over idist sum_freq = SUM(pop%pop_grid(g)%freq) IF (sum_freq.GT.0.0) pop%pop_grid(g)%freq = pop%pop_grid(g)%freq/sum_freq ENDDO END SUBROUTINE GetPatchFrequencies !============================================================================= SUBROUTINE GetDiagnostics(pop,LAI,Cleaf,Croot,disturbance_interval, it, precip) ! Gets diagnostic data for current landscape structure IMPLICIT NONE TYPE(POP_TYPE), INTENT(INOUT) :: POP REAL(dp), INTENT(IN) :: LAI(:) REAL(dp), INTENT(IN) :: Cleaf(:) REAL(dp), INTENT(IN) :: Croot(:) INTEGER(i4b), INTENT(IN) :: disturbance_interval(:,:) REAL(dp), INTENT(IN), OPTIONAL :: precip(:) INTEGER(i4b), INTENT(IN) :: it(:) INTEGER(i4b) :: P, g,i,j,ct, ct_highres REAL(dp) :: limits(HEIGHT_BINS+1) REAL(dp) :: ht, htmax, cmass_stem,densindiv, freq, freq_old CHARACTER(len=12) :: string1, string2 CHARACTER(len=9) :: fmt INTEGER(i4b) :: npatch_active ! number of active patches INTEGER(i4b) :: np, nc, i_height, id REAL(dp) :: diam,basal, cump REAL(dp) :: patch_crown_area(NPATCH2D), patch_crown_cover(NPATCH2D) REAL(dp), ALLOCATABLE :: height_list(:), height_list_weight(:) REAL(dp) :: height_copy, weight_copy, Pwc, FAVD INTEGER(i4b), PARAMETER :: HEIGHT_BINS_highres=100 ! bins for assessing height_max REAL(dp), ALLOCATABLE :: limits_highres(:), DENSINDIV_HIGHRES(:) REAL(dp) :: a, b, C, res_flux, fire_flux, cat_flux, vol, vol1, tmp, tmp2 INTEGER :: arg1 fmt = '(f5.1)' limits(1) = 0. IF(.NOT.ALLOCATED(limits_highres)) ALLOCATE(limits_highres(HEIGHT_BINS_highres+1)) IF(.NOT.ALLOCATED(DENSINDIV_HIGHRES)) ALLOCATE(DENSINDIV_HIGHRES(HEIGHT_BINS_highres)) limits_highres(1) = 0. np = SIZE(Pop%pop_grid) DO g=1,np npatch_active = NPATCH2D IF (MAX_HEIGHT_SWITCH.EQ.1) THEN ALLOCATE(height_list(NPATCH2D*NCOHORT_MAX)) ALLOCATE(height_list_weight(NPATCH2D*NCOHORT_MAX)) ENDIF ! IF(.NOT.ALLOCATED(MASK)) ALLOCATE(MASK(POP%pop_grid%npatch_active)) DO i=1,HEIGHT_BINS limits(i+1) = BIN_POWER**REAL(i) WRITE(string1,fmt) (limits(i)) WRITE(string2,fmt) (limits(i+1)) pop%pop_grid(g)%bin_labels(i) = 'Height_'//TRIM(ADJUSTL(string1))//'-'//TRIM(ADJUSTL(string2))//'m' pop%pop_grid(g)%cmass_stem_bin(i) = 0.0 pop%pop_grid(g)%densindiv_bin(i) = 0.0 pop%pop_grid(g)%cmass_stem_bin(i) = 0.0 pop%pop_grid(g)%height_bin(i) = REAL(limits(i)+limits(i+1))/2. pop%pop_grid(g)%diameter_bin(i) = ((REAL(limits(i))/Kbiometric)**(3/2)+(REAL(limits(i+1))/Kbiometric)**(3/2))/2. ENDDO DO i=1,HEIGHT_BINS_highres limits_highres(i+1) = REAL(i) ENDDO IF (MAX_HEIGHT_SWITCH.EQ.1) THEN height_list = 0.0 height_list_weight = 0.0 ENDIF i_height = 0 pop%pop_grid(g)%cmass_sum_old = pop%pop_grid(g)%cmass_sum pop%pop_grid(g)%csapwood_sum_old = pop%pop_grid(g)%csapwood_sum pop%pop_grid(g)%cmass_sum = 0.0 pop%pop_grid(g)%csapwood_sum = 0.0 pop%pop_grid(g)%cheartwood_sum = 0.0 pop%pop_grid(g)%height_mean = 0.0 pop%pop_grid(g)%fire_mortality = 0.0 pop%pop_grid(g)%cat_mortality = 0.0 pop%pop_grid(g)%res_mortality = 0.0 pop%pop_grid(g)%stress_mortality = 0.0 pop%pop_grid(g)%crowding_mortality = 0.0 pop%pop_grid(g)%sapwood_loss = 0.0 pop%pop_grid(g)%sapwood_area_loss = 0.0 pop%pop_grid(g)%growth = 0.0 pop%pop_grid(g)%area_growth = 0.0 pop%pop_grid(g)%basal_area = 0.0 pop%pop_grid(g)%densindiv = 0.0 pop%pop_grid(g)%height_max = 0.0 pop%pop_grid(g)%crown_cover = 0.0 pop%pop_grid(g)%crown_area = 0.0 pop%pop_grid(g)%sapwood_area = 0.0 pop%pop_grid(g)%crown_volume = 0.0 densindiv_highres = 0.0 ! loop through patches DO P = 1, npatch_active pop%pop_grid(g)%patch(p)%biomass = 0.0 pop%pop_grid(g)%patch(p)%sapwood = 0.0 pop%pop_grid(g)%patch(p)%sapwood_area = 0.0 pop%pop_grid(g)%patch(p)%heartwood = 0.0 pop%pop_grid(g)%patch(p)%layer(1)%biomass = 0.0 pop%pop_grid(g)%patch(p)%layer(1)%density = 0.0 patch_crown_area(p) = 0.0 patch_crown_cover(p) = 0.0 tmp2 = SUM(pop%pop_grid(g)%patch(p)%layer(1)%cohort(1:pop%pop_grid(g)%patch(p)%layer(1)%ncohort)%sapwood_area) freq = pop%pop_grid(g)%freq(pop%pop_grid(g)%patch(p)%id) freq_old = pop%pop_grid(g)%freq_old(pop%pop_grid(g)%patch(p)%id) ! loop through cohorts DO i = 1, pop%pop_grid(g)%patch(p)%layer(1)%ncohort cmass_stem = pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%biomass densindiv = pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%density IF ( PRESENT(precip) ) THEN CALL GET_ALLOMETRY( ALLOM_SWITCH, cmass_stem, densindiv, ht, diam, basal, precip(g)) ELSE CALL GET_ALLOMETRY( ALLOM_SWITCH, cmass_stem, densindiv, ht, diam, basal ) ENDIF pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%height = ht pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%diameter = diam ! basal area in each cohort pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%basal_area = basal ! sapwood area in each cohort pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%sapwood_area = basal - & ! m2 ha-1 pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%heartwood/(ht*WD)*1.e4 pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%sapwood_area = & MAX(pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%sapwood_area, 0.0) ! get bin ct = 1 DO j=1,HEIGHT_BINS IF (ht.GT.limits(j)) ct = j ENDDO ! bins ! get high res bin ct_highres = 1 DO j=1,HEIGHT_BINS_highres IF (ht.GT.limits_highres(j)) ct_highres = j ENDDO ! bins pop%pop_grid(g)%patch(p)%layer(1)%biomass = pop%pop_grid(g)%patch(p)%layer(1)%biomass + cmass_stem pop%pop_grid(g)%patch(p)%layer(1)%density = pop%pop_grid(g)%patch(p)%layer(1)%density + densindiv IF (diam*100.0.GT.1.) THEN patch_crown_area(p) = patch_crown_area(p) + densindiv*PI*(diam*100.*0.1492)**2 ! uses GC relationship pop%pop_grid(g)%crown_volume = pop%pop_grid(g)%crown_volume + & freq*densindiv*(4./3.)*PI*(diam*100.*0.1492)**2*(1.5*(diam*100.*0.1492)) ! assumes vertical radius = 1.5 * horizontal radius ENDIF IF (diam*100..GT.5.) THEN IF (ALLOM_SWITCH.EQ.1) THEN ! assumes crown radius (m) = 0.1492 * dbh (cm) (from G. Cook, pers. comm.) ! assumes vertical radius = 1.5 * horizontal radius pop%pop_grid(g)%crown_volume = pop%pop_grid(g)%crown_volume + & freq*densindiv*(4./3.)*PI*(diam*100.*0.1492)**2*(1.5*(diam*100.*0.1492)) ELSE ! global allometry ! assumes vertical radius = 1.5 * horizontal radius pop%pop_grid(g)%crown_volume = pop%pop_grid(g)%crown_volume + & freq*densindiv*(4./3.)*PI*1.5*((k_allom1 * diam ** k_rp )/PI)**1.5 ENDIF ENDIF pop%pop_grid(g)%patch(p)%sapwood = pop%pop_grid(g)%patch(p)%sapwood + & pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%sapwood pop%pop_grid(g)%patch(p)%sapwood_area = pop%pop_grid(g)%patch(p)%sapwood_area + & pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%sapwood_area pop%pop_grid(g)%patch(p)%heartwood = pop%pop_grid(g)%patch(p)%heartwood + & pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%heartwood pop%pop_grid(g)%patch(p)%biomass = pop%pop_grid(g)%patch(p)%biomass + cmass_stem pop%pop_grid(g)%cmass_stem_bin(ct) = pop%pop_grid(g)%cmass_stem_bin(ct) + freq*cmass_stem pop%pop_grid(g)%densindiv_bin(ct) = pop%pop_grid(g)%densindiv_bin(ct) + freq*densindiv pop%pop_grid(g)%cmass_sum = pop%pop_grid(g)%cmass_sum + freq*cmass_stem pop%pop_grid(g)%csapwood_sum = pop%pop_grid(g)%csapwood_sum + & freq*pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%sapwood pop%pop_grid(g)%cheartwood_sum = pop%pop_grid(g)%cheartwood_sum + & freq*pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%heartwood pop%pop_grid(g)%densindiv = pop%pop_grid(g)%densindiv + freq*densindiv pop%pop_grid(g)%height_mean = pop%pop_grid(g)%height_mean + ht*freq*densindiv pop%pop_grid(g)%basal_area = pop%pop_grid(g)%basal_area +basal*freq densindiv_highres(ct_highres) = densindiv_highres(ct_highres) + freq*densindiv IF (MAX_HEIGHT_SWITCH.EQ.1) THEN i_height = i_height+1 height_list(i_height) = ht height_list_weight(i_height) = densindiv*freq ENDIF ENDDO ! cohorts pop%pop_grid(g)%stress_mortality = pop%pop_grid(g)%stress_mortality + & freq*pop%pop_grid(g)%patch(p)%stress_mortality pop%pop_grid(g)%crowding_mortality = pop%pop_grid(g)%crowding_mortality + & freq*pop%pop_grid(g)%patch(p)%crowding_mortality pop%pop_grid(g)%sapwood_area = pop%pop_grid(g)%sapwood_area + & freq*SUM(pop%pop_grid(g)%patch(p)%layer(1)%cohort(1:pop%pop_grid(g)%patch(p)%layer(1)%ncohort)%sapwood_area) pop%pop_grid(g)%sapwood_loss = pop%pop_grid(g)%sapwood_loss + & freq*pop%pop_grid(g)%patch(p)%sapwood_loss pop%pop_grid(g)%sapwood_loss = pop%pop_grid(g)%sapwood_loss + & pop%pop_grid(g)%patch(p)%sapwood_old*(freq_old-freq) pop%pop_grid(g)%sapwood_area_loss = pop%pop_grid(g)%sapwood_area_loss + & freq*pop%pop_grid(g)%patch(p)%sapwood_area_loss pop%pop_grid(g)%sapwood_area_loss = pop%pop_grid(g)%sapwood_area_loss + & pop%pop_grid(g)%patch(p)%sapwood_area_old*(freq_old-freq) pop%pop_grid(g)%cat_mortality = pop%pop_grid(g)%cat_mortality + & freq*pop%pop_grid(g)%patch(p)%cat_mortality pop%pop_grid(g)%res_mortality = pop%pop_grid(g)%res_mortality + & pop%pop_grid(g)%patch(p)%biomass_old*(freq_old-freq) pop%pop_grid(g)%fire_mortality = pop%pop_grid(g)%fire_mortality + & freq*pop%pop_grid(g)%patch(p)%fire_mortality pop%pop_grid(g)%growth = pop%pop_grid(g)%growth + freq_old*pop%pop_grid(g)%patch(p)%growth !write(*,*) 'freq_old',freq_old pop%pop_grid(g)%area_growth = pop%pop_grid(g)%area_growth + & freq*pop%pop_grid(g)%patch(p)%area_growth ENDDO ! patches IF (INTERP_SWITCH==1.AND.NDISTURB.EQ.2) THEN !CALL INTERPOLATE_BIOMASS_2D(pop, disturbance_interval,it) CALL INTERPOLATE_BIOMASS_2D(pop, disturbance_interval,it(g),g) ELSEIF (INTERP_SWITCH==1.AND.NDISTURB.EQ.1) THEN CALL INTERPOLATE_BIOMASS_1D(pop, disturbance_interval,it(g),g) ENDIF arg1 = NYEAR_SMOOTH-(NYEAR_SMOOTH/2) IF (SMOOTH_SWITCH==1) THEN IF (it(g).LE.NYEAR_SMOOTH-NYEAR_SMOOTH/2) THEN CALL SMOOTH_FLUX(POP,g,it(g)) ELSE CALL SMOOTH_FLUX(POP,g,INT(arg1,i4b)) ENDIF ENDIF ! leaf area index in each cohort DO P = 1, npatch_active freq = pop%pop_grid(g)%freq(pop%pop_grid(g)%patch(p)%id) ! loop through cohorts DO i = 1, pop%pop_grid(g)%patch(p)%layer(1)%ncohort cmass_stem = pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%biomass densindiv = pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%density basal=PI*(pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%diameter/2.0)* & (pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%diameter/2.0)*densindiv*1.e4 ! leaf area index in each cohort pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%LAI = LAI(g) * & MIN(pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%sapwood_area & /MAX(pop%pop_grid(g)%sapwood_area,1e-3), 10.0_dp) pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%Cleaf = Cleaf(g) * & MIN(pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%sapwood_area & /MAX(pop%pop_grid(g)%sapwood_area,1.e-3), 10.0_dp) pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%Croot = Croot(g) * & MIN(pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%sapwood_area & /MAX(pop%pop_grid(g)%sapwood_area,1e-3), 10.0_dp) ENDDO ! cohorts pop%pop_grid(g)%patch(p)%LAI = SUM(pop%pop_grid(g)%patch(p)%layer(1)% & cohort(1:pop%pop_grid(g)%patch(p)%layer(1)%ncohort)%LAI) ENDDO ! patches ! PGap = (1-fcover) calculation IF (pop%pop_grid(g)%crown_volume>0.0) THEN FAVD = LAI(g)/pop%pop_grid(g)%crown_volume ! foliage area volume density ELSE FAVD = 0.0 ENDIF DO P = 1, npatch_active freq = pop%pop_grid(g)%freq(pop%pop_grid(g)%patch(p)%id) nc = pop%pop_grid(g)%patch(p)%layer(1)%ncohort ! loop through cohorts DO i = 1, nc cmass_stem = pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%biomass densindiv = pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%density IF ( PRESENT(precip) ) THEN CALL GET_ALLOMETRY( ALLOM_SWITCH, cmass_stem, densindiv, ht, diam, basal, precip(g)) ELSE CALL GET_ALLOMETRY( ALLOM_SWITCH, cmass_stem, densindiv, ht, diam, basal ) ENDIF IF (diam*100.GT.1.) THEN IF (ALLOM_SWITCH.EQ.1) THEN ! assumes crown radius (m) = 0.1492 * dbh (cm) (from G. Cook, pers. comm.) pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%crown_area = densindiv*PI*(diam*100.*0.1492)**2 Pwc = EXP(-0.5 * pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%LAI/ & pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%crown_area) pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%crown_area = densindiv*PI*(diam*100.*0.1492)**2*(1.-Pwc) ELSE ! global allometry pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%crown_area = & densindiv*PI*(((k_allom1 * diam ** k_rp )/PI)**0.5)**2 Pwc = EXP(MAX(-0.5 * pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%LAI/ & MAX(pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%crown_area,1.e-3),-20.0)) pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%crown_area = & pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%crown_area*(1.-Pwc) ! *1.4142 ENDIF ELSE pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%crown_area = & 0.5*pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%LAI ! *1.4142 ENDIF pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%crown_area= & MAX(pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%crown_area,0.01) pop%pop_grid(g)%crown_area = pop%pop_grid(g)%crown_area + & freq*pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%crown_area IF (i.EQ.1) THEN ! top cohort pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%Pgap = & EXP(-pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%crown_area) pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%frac_interception = & 1- EXP(-pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%crown_area) ELSE pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%Pgap = & pop%pop_grid(g)%patch(p)%layer(1)%cohort(i-1)%Pgap* & EXP(-pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%crown_area) pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%frac_interception = & pop%pop_grid(g)%patch(p)%layer(1)%cohort(i-1)%Pgap - & pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%Pgap ENDIF pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%respiration_scalar = & pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%sapwood/shootfrac/CtoNw + & pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%Cleaf/CtoNl + & pop%pop_grid(g)%patch(p)%layer(1)%cohort(i)%Croot/CtoNr ENDDO ! cohorts IF (nc>0) THEN pop%pop_grid(g)%patch(p)%pgap = & pop%pop_grid(g)%patch(p)%layer(1)%cohort(nc)%Pgap ELSE pop%pop_grid(g)%patch(p)%pgap = 1 ENDIF ENDDO ! patches pop%pop_grid(g)%Kclump = MAX(pop%pop_grid(g)%crown_area/(0.5*LAI(g)),0.1_dp) pop%pop_grid(g)%crown_cover = 1.-EXP(-pop%pop_grid(g)%crown_area) pop%pop_grid(g)%height_mean = pop%pop_grid(g)%height_mean/MAX(pop%pop_grid(g)%densindiv,1.0e-5) ! Height Diagnostics IF (MAX_HEIGHT_SWITCH.EQ.0) THEN ! Set landscape maximum height to centre of bin with <5% of trees in a bin of higher size classes cump = 0. j = 1 DO WHILE (cump.LT.0.95) cump = cump + pop%pop_grid(g)%densindiv_bin(j)/MAX(pop%pop_grid(g)%densindiv,1.0e-5) pop%pop_grid(g)%height_max = pop%pop_grid(g)%height_bin(j) j = j+1 ENDDO ELSEIF (MAX_HEIGHT_SWITCH.EQ.1) THEN ! sort height list DO i = 1, i_height -1 DO j = i+1, i_height IF (height_list(i).GT.height_list(j)) THEN height_copy = height_list(i) weight_copy = height_list_weight(i) height_list(i) = height_list(j) height_list_weight(i) = height_list_weight(j) height_list(j) = height_copy height_list_weight(j) = weight_copy ENDIF ENDDO ENDDO ! end sort height list ! normailse height list weights height_list_weight=height_list_weight/SUM(height_list_weight(1:i_height)) cump = 0. j = 1 DO WHILE (cump.LT.0.95) cump = cump + height_list_weight(j) pop%pop_grid(g)%height_max = height_list(j) j = j+1 ENDDO DEALLOCATE(height_list) DEALLOCATE(height_list_weight) ELSEIF (MAX_HEIGHT_SWITCH.EQ.2) THEN cump = 0. j = 1 densindiv_highres= densindiv_highres/MAX(SUM(densindiv_highres),1.0e-5) DO WHILE ((cump.LT.0.95).AND.(j.LE.HEIGHT_BINS_highres)) cump = cump + densindiv_highres(j) pop%pop_grid(g)%height_max = (limits_highres(j+1) + limits_highres(j))/2. j = j+1 ENDDO ENDIF !deallocate(MASK) ENDDO ! end loop over grid cells END SUBROUTINE GetDiagnostics !============================================================================= SUBROUTINE Patch_partial_disturb(pop,idisturb,intensity,precip,frac_intensity1) IMPLICIT NONE TYPE(POP_TYPE), INTENT(INOUT) :: POP INTEGER(i4b), INTENT(IN) :: idisturb REAL(dp), INTENT(IN) :: intensity(:,:) REAL(dp), INTENT(IN), OPTIONAL :: precip(:), frac_intensity1(:) INTEGER(i4b) :: j, k, i, g, c, nc, np INTEGER(i4b) :: ivec(NCOHORT_MAX) REAL(dp) :: ht, diam REAL(dp) :: Psurvival_l, Psurvival_s, Psurvival, char_height np = SIZE(Pop%pop_grid) !Print*,"CLN Hier PROBLEM PSurv1 unten l1345 wegen PRESENT() " ! Kills a fraction of biomass in patch when prescribed disturbance interval is reached DO j=1,np DO k=1,NPATCH pop%pop_grid(j)%patch(k)%fire_mortality = 0.0 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. & (pop%pop_grid(j)%patch(k)%disturbance_interval(idisturb).EQ.pop%pop_grid(j)%patch(k)%age(idisturb))) THEN ! loop through cohorts ivec = 0 nc = 0 DO c = 1, pop%pop_grid(j)%patch(k)%layer(1)%ncohort ! kill fraction of each cohort char_height = 3.7*(1.-EXP(-0.19*Intensity(j,1))) ht = pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%height diam = pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%diameter*100. ! diameter in cm IF ((ht.GT.8.5).AND.(ht.GT.char_height)) THEN Psurvival_s =(-0.0011*Intensity(j,1) -0.00002)*ht & +(0.0075*Intensity(j,1)+1.) ELSEIF ((ht.LE.8.5).AND.(ht.GT.char_height)) THEN Psurvival_s =(0.0178*Intensity(j,1) + 0.0144)*ht & + (-0.1174*Intensity(j,1)+0.9158) ELSE Psurvival_s = 0.0 ENDIF Psurvival_s = MIN(Psurvival_s,1.0_dp) Psurvival_s = MAX(Psurvival_s,1e-3_dp) Psurvival = Psurvival_s IF (PRESENT(frac_intensity1)) THEN char_height = 3.7*(1.-EXP(-0.19*Intensity(j,2))) ht = pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%height IF ((ht.GT.8.5).AND.(ht.GT.char_height)) THEN Psurvival_s =(-0.0011*Intensity(j,2) -0.00002)*ht & +(0.0075*Intensity(j,2)+1.) ELSEIF ((ht.LE.8.5).AND.(ht.GT.char_height)) THEN Psurvival_s =(0.0178*Intensity(j,2) + 0.0144)*ht & + (-0.1174*Intensity(j,2)+0.9158) ELSE Psurvival_s = 0.0 ENDIF Psurvival_s = MIN(Psurvival_s,1.0_dp) Psurvival_s = MAX(Psurvival_s,1e-3_dp) Psurvival = Psurvival_s*(1.-frac_intensity1(j)) + Psurvival*frac_intensity1(j) ENDIF !Psurvival = 1 ! test pop%pop_grid(j)%patch(k)%fire_mortality = pop%pop_grid(j)%patch(k)%fire_mortality + & (1.-Psurvival)*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 + & (1.-Psurvival)*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 + & (1.-Psurvival)*pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood_area pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass = & Psurvival*pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood = & Psurvival*pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%heartwood = & 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 ! 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 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)%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 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 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%heartwood = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood_area = 0.0 ELSE 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 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 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%sapwood = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%sapwood_area = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%heartwood = 0.0 ENDIF pop%pop_grid(j)%patch(k)%age(idisturb) = 0 pop%pop_grid(j)%patch(k)%first_disturbance_year(idisturb) = 0 ENDIF ENDDO ENDDO END SUBROUTINE Patch_partial_disturb !============================================================================= SUBROUTINE Patch_partial_disturb2(pop,idisturb,precip) IMPLICIT NONE TYPE(POP_TYPE), INTENT(INOUT) :: POP INTEGER(i4b), INTENT(IN) :: idisturb REAL(dp), INTENT(IN), OPTIONAL :: precip(:) INTEGER(i4b) :: j, k, i, g, c, nc, np INTEGER(i4b) :: ivec(NCOHORT_MAX) REAL(dp) :: ht, diam REAL(dp) :: Psurvival_l, Psurvival_s, Psurvival, char_height, frac_mort, Pmort np = SIZE(Pop%pop_grid) ! Kills a fraction (80%) biomass in patch when prescribed disturbance interval is reached DO j=1,np DO k=1,NPATCH2D pop%pop_grid(j)%patch(k)%cat_mortality = 0. ! 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. & (pop%pop_grid(j)%patch(k)%disturbance_interval(idisturb).EQ.pop%pop_grid(j)%patch(k)%age(idisturb))) THEN ! loop through cohorts ivec = 0 nc = 0 frac_mort = 0.0 pop%pop_grid(j)%patch(k)%cat_mortality = 0.0 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 * pop%pop_grid(j)%patch(k)%Layer(1)%biomass ) THEN Pmort = MIN( (0.8*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 ENDIF Psurvival = 1- Pmort pop%pop_grid(j)%patch(k)%cat_mortality = pop%pop_grid(j)%patch(k)%cat_mortality + & (1.-Psurvival)*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 + & (1.-Psurvival)*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 + & (1.-Psurvival)*pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood_area pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass = & Psurvival*pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%biomass pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood = & Psurvival*pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%sapwood pop%pop_grid(j)%patch(k)%layer(1)%cohort(c)%heartwood = & 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 ! 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 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)%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 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 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%heartwood = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(c)%sapwood_area = 0.0 ELSE 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 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 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%sapwood = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%sapwood_area = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(nc+1:NCOHORT_MAX)%heartwood = 0.0 ENDIF pop%pop_grid(j)%patch(k)%age(idisturb) = 0 pop%pop_grid(j)%patch(k)%first_disturbance_year(idisturb) = 0 ENDIF ENDDO ENDDO END SUBROUTINE Patch_partial_disturb2 !============================================================================= SUBROUTINE Patch_disturb(pop,idisturb,precip) IMPLICIT NONE TYPE(POP_TYPE), INTENT(INOUT) :: POP REAL(dp), INTENT(IN), OPTIONAL :: precip(:) !INTEGER(i4b), INTENT(IN) :: it(:),idisturb INTEGER(i4b), INTENT(IN) :: idisturb INTEGER(i4b) :: j, k, np, nc np = SIZE(Pop%pop_grid) ! Kills all biomass in patch when prescribed disturbance interval is reached ! Should be called after accounting for this year DO j=1,np DO k=1,NPATCH2D pop%pop_grid(j)%patch(k)%cat_mortality = 0. IF (pop%pop_grid(j)%patch(k)%first_disturbance_year(idisturb).NE.0) THEN IF ((pop%pop_grid(j)%patch(k)%first_disturbance_year(idisturb).EQ.pop%pop_grid(j)%patch(k)%age(idisturb)).OR. & (pop%pop_grid(j)%patch(k)%disturbance_interval(idisturb).EQ.pop%pop_grid(j)%patch(k)%age(idisturb)) ) THEN ! kill entire layer nc = pop%pop_grid(j)%patch(k)%layer(1)%ncohort ! pop%pop_grid(j)%patch(k)%fire_mortality = SUM(pop%pop_grid(j)%patch(k)%layer(1)%cohort(1:nc)%biomass) 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)%sapwood_loss = pop%pop_grid(j)%patch(k)%sapwood_loss + & 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 pop%pop_grid(j)%patch(k)%layer(1:NLayer)%density = 0.0 pop%pop_grid(j)%patch(k)%layer(1:NLayer)%hmean = 0.0 pop%pop_grid(j)%patch(k)%layer(1:NLayer)%hmax = 0.0 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 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 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%sapwood = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%sapwood_area = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%heartwood = 0.0 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 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%height = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%pgap = 1.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_resource_uptake = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_light_uptake = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_interception = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_respiration = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_NPP = 0.0 pop%pop_grid(j)%patch(k)%area_growth = 0.0 pop%pop_grid(j)%patch(k)%pgap = 1.0 ! understorey recruitment IF (PRESENT(precip)) THEN CALL layer_recruitment_single_patch(pop,k,j,precip) ELSE CALL layer_recruitment_single_patch(pop,k,j) ENDIF ENDIF ELSEIF (pop%pop_grid(j)%patch(k)%disturbance_interval(idisturb).EQ.pop%pop_grid(j)%patch(k)%age(idisturb)) THEN ! kill entire layer nc = pop%pop_grid(j)%patch(k)%layer(1)%ncohort pop%pop_grid(j)%patch(k)%sapwood_loss = pop%pop_grid(j)%patch(k)%sapwood_loss + & 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)%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 pop%pop_grid(j)%patch(k)%layer(1:NLayer)%density = 0.0 pop%pop_grid(j)%patch(k)%layer(1:NLayer)%hmean = 0.0 pop%pop_grid(j)%patch(k)%layer(1:NLayer)%hmax = 0.0 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 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 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%sapwood = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%sapwood_area = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%heartwood = 0.0 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 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%height = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%pgap = 1.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_resource_uptake = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_light_uptake = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_interception = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_respiration = 0.0 pop%pop_grid(j)%patch(k)%Layer(1)%cohort(1:NCOHORT_MAX)%frac_NPP = 0.0 pop%pop_grid(j)%patch(k)%area_growth = 0.0 pop%pop_grid(j)%patch(k)%pgap = 1.0 ! understorey recruitment IF (PRESENT(precip)) THEN CALL layer_recruitment_single_patch(pop,k,j,precip) ELSE CALL layer_recruitment_single_patch(pop,k,j) ENDIF ENDIF ENDDO ENDDO END SUBROUTINE Patch_disturb !============================================================================= SUBROUTINE layer_recruitment(pop,precip) IMPLICIT NONE TYPE(POP_TYPE), INTENT(INOUT) :: POP REAL(dp), INTENT(IN), OPTIONAL :: precip(:) REAL(dp) :: f, mu, densindiv, cmass, ht REAL(dp) :: tmp, cmass_stem_sum,cmass_stem_inc INTEGER(i4b) :: j, k,c, ncohort, np REAL(dp) :: diam,basal np = SIZE(Pop%pop_grid) DO j=1,np DO k=1,NPATCH2D IF (RECRUIT_SWITCH==0) THEN pop%pop_grid(j)%patch(k)%factor_recruit = EXP(-0.6*((pop%pop_grid(j)%patch(k)%Layer(1)%biomass)**(0.6667))) ELSEIF (RECRUIT_SWITCH==1) THEN pop%pop_grid(j)%patch(k)%factor_recruit = MAX(pop%pop_grid(j)%patch(k)%pgap,1e-3) ENDIF f = pop%pop_grid(j)%patch(k)%factor_recruit mu=EXP(MAX(FULTON_ALPHA*(1.0-2*THETA_recruit/(f+1-SQRT((f+1)*(f+1)-4*THETA_recruit*f))),-50.0)); densindiv=DENSINDIV_MAX*mu; cmass=CMASS_STEM_INIT*densindiv/DENSINDIV_MAX; !COMMLN below: should not be cohort +1 or .LE. ! IF (cmass>EPS*10..AND.densindiv>DENSINDIV_MIN.AND.(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 IF ( PRESENT(precip) ) THEN CALL GET_ALLOMETRY( ALLOM_SWITCH, cmass, densindiv, ht, diam, basal, precip(j)) ELSE CALL GET_ALLOMETRY( ALLOM_SWITCH, cmass, densindiv, ht, diam, basal ) ENDIF pop%pop_grid(j)%patch(k)%layer(1)%cohort(ncohort)%height = ht pop%pop_grid(j)%patch(k)%layer(1)%cohort(ncohort)%diameter = diam ENDIF ENDDO ENDDO END SUBROUTINE layer_recruitment !============================================================================= SUBROUTINE layer_recruitment_single_patch(pop, index, grid_index,precip) IMPLICIT NONE TYPE(POP_TYPE), INTENT(INOUT) :: POP REAL(dp), INTENT(IN), OPTIONAL :: precip(:) INTEGER(i4b), INTENT(IN) :: index, grid_index REAL(dp) :: f, mu, densindiv, cmass, ht REAL(dp) :: tmp, cmass_stem_sum,cmass_stem_inc INTEGER(i4b) :: j, k,c, ncohort, np REAL(dp) :: diam,basal np = SIZE(Pop%pop_grid) 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*((pop%pop_grid(j)%patch(k)%Layer(1)%biomass)**(0.6667))) 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 ENDIF f = pop%pop_grid(j)%patch(k)%factor_recruit mu=EXP(FULTON_ALPHA*(1.0-2*THETA_recruit/(f+1-SQRT((f+1)*(f+1)-4*THETA_recruit*f)))); densindiv=DENSINDIV_MAX*mu; cmass=CMASS_STEM_INIT*densindiv/DENSINDIV_MAX; ! write(*,*) 'layer_recruitment_single_patch', densindiv, pop%pop_grid(j)%patch(k)%Layer(1)%biomass IF (cmass>EPS*10..AND.densindiv>DENSINDIV_MIN.AND.(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 IF ( PRESENT(precip) ) THEN CALL GET_ALLOMETRY( ALLOM_SWITCH, cmass, densindiv, ht, diam, basal, precip(j)) ELSE CALL GET_ALLOMETRY( ALLOM_SWITCH, cmass, densindiv, ht, diam, basal ) ENDIF pop%pop_grid(j)%patch(k)%layer(1)%cohort(ncohort)%height = ht pop%pop_grid(j)%patch(k)%layer(1)%cohort(ncohort)%diameter = diam ENDIF ENDDO ENDDO END SUBROUTINE layer_recruitment_single_patch !============================================================================= ! Exponential distribution ! Returns probability of a given time-between-events (x) ! Given a Poisson process with expected frequency (events per unit time) lambda ! Reference: http://en.wikipedia.org/wiki/Exponential_distribution ! Use to determine average age (x, years) of patches with a given random disturbance ! frequency lambda (disturbances per year) REAL(dp) FUNCTION Exponential(lambda, x) IMPLICIT NONE INTEGER(i4b), INTENT(IN) :: x REAL(dp), INTENT(IN) :: lambda IF (x.LT.0) THEN ! Shouldn't happen but ... Exponential=0.0 ELSE Exponential=lambda*EXP(-lambda*x) ENDIF END FUNCTION Exponential !============================================================================= ! Exponential distribution ! Returns probability of a given time-between-events (x) ! Given a Poisson process with expected frequency (events per unit time) lambda ! Reference: http://en.wikipedia.org/wiki/Exponential_distribution ! Use to determine average age (x, years) of patches with a given random disturbance ! frequency lambda (disturbances per year) REAL(dp) FUNCTION REALExponential(lambda, x) IMPLICIT NONE REAL(dp), INTENT(IN) :: x REAL(dp), INTENT(IN) :: lambda IF (x.LT.0) THEN ! Shouldn't happen but ... REALExponential=0.0 ELSE REALExponential=lambda*EXP(-lambda*x) ENDIF END FUNCTION REALExponential !============================================================================= REAL(dp) FUNCTION CumExponential(lambda, x) IMPLICIT NONE REAL(dp), INTENT(IN) :: x REAL(dp), INTENT(IN) :: lambda IF (x.LT.0) THEN ! Shouldn't happen but ... CumExponential=0.0 ELSE CumExponential=1.-EXP(-lambda*x) ENDIF END FUNCTION CumExponential !============================================================================= REAL(dp) FUNCTION Factorial(n) IMPLICIT NONE INTEGER , INTENT(IN) :: n INTEGER :: i REAL(dp) :: Ans Ans = 1. DO i = 1, n Ans = Ans *REAL(i,dp) END DO Factorial = Ans END FUNCTION Factorial !============================================================================= ! ALLOMETRY !============================================================================= SUBROUTINE GET_ALLOMETRY( ALLOM_SWITCH, biomass, density, ht, diam, basal, precip ) IMPLICIT NONE INTEGER(i4b), INTENT(IN) :: ALLOM_SWITCH REAL(dp), INTENT(IN) :: biomass REAL(dp), INTENT(IN) :: density REAL(dp), INTENT(IN), OPTIONAL :: precip REAL(dp), INTENT(OUT):: ht, diam, basal ! Standard Allometry IF (ALLOM_SWITCH.EQ.0) THEN ht = (Kbiometric**(3.0/4.0))*(4.*biomass/(MAX(density,1e-5)*WD*PI))**(1.0/4.0) diam = (ht/Kbiometric)**(1.5) basal= PI * (diam/2.0) * (diam/2.0) * density * 1.e4 ! Top-End Allometry following G.Cook ELSEIF (ALLOM_SWITCH.EQ.1.AND.PRESENT(precip)) THEN ht =GetHeight(precip,biomass,density) CALL Allometry(ht,biomass,density,diam,basal) ! Allometry following Williams 2005, Model 5b ELSEIF ( ALLOM_SWITCH.EQ.2 ) THEN CALL Williams_Allometry(biomass,density,ht,diam,basal) ELSE WRITE(*,*)"Invalid Allometry settings in POP!" WRITE(*,*)"ALLOM_SWITCH = ",ALLOM_SWITCH WRITE(*,*)"Precip present = ",PRESENT(precip) STOP -1 ENDIF END SUBROUTINE GET_ALLOMETRY !============================================================================= ! TOP-END ALLOMETRY STARTS HERE !============================================================================= ! Tree height based on precipitation and Gary Cook Top-End allometry ! Bisection solution for tree height (m) based on modified height-DBH relationship ! from Garry Cook (pers. comm. 15/4/2013) ! "I have been using H=0.054xExp(0.0014xRF)xD + 4.05*exp(-0.00032*rf) with Rf in mm, D in cm and height in m." ! Since the above expression does not go to zero at diameter=0 it is linked to a simple linear equation ! for initial height growth (H=50*D with D in m) using a non-rectangular hyperbola to smooth between the two ! Mathematical derivation: see POP documentation ! Arguments: ! precip = annual precipitation (mm) ! biomass = tree stem C biomass across patch (kgC/m2) ! density = tree density (indiv/m2) REAL(dp) FUNCTION GetHeight(precip,biomass,density) IMPLICIT NONE REAL(dp), INTENT(IN) :: precip REAL(dp), INTENT(IN) :: biomass REAL(dp), INTENT(IN) :: density REAL(dp),PARAMETER:: THETA=0.99 ! Shape parameter, should be slightly <1 REAL(dp),PARAMETER:: HMIN=0.001 ! min bound for tree height REAL(dp),PARAMETER:: HMAX=100 ! max bound for tree height REAL(dp),PARAMETER:: EPS=0.01 ! precision of the root INTEGER(i4b), PARAMETER :: MAXTRIES=25 REAL(dp) :: alpha,beta,delta,rh,st,x1,x2,rtbis,dx,fmid,xmid,lhs,rhs INTEGER(i4b) :: b alpha=4.05*EXP(-0.00032*precip) beta=5.4*EXP(0.0014*precip) delta=2.0*SQRT(biomass/density/WD/PI) x1=HMIN x2=HMAX rtbis=x1 dx=x2-x1 b=0 fmid=EPS+1.0 DO WHILE (ABS(dx).GT.EPS.AND.b.LE.MAXTRIES) b=b+1 dx=dx*0.5 xmid=rtbis+dx ! Evaluate LHS-RHS at height=xmid ! LHS-RHS should increase with increasing height lhs=xmid rh=1.0/SQRT(xmid) st=alpha+beta*delta*rh+100*delta*rh rhs=1.0/2.0/THETA* & (st-SQRT(st*st-400*THETA*alpha*delta*rh- & 400*THETA*beta*delta*delta/xmid)) fmid=lhs-rhs IF (fmid.LT.0.0) rtbis=xmid ENDDO GetHeight=xmid END FUNCTION GetHeight !============================================================================= SUBROUTINE INTERPOLATE_BIOMASS_1D(pop, disturbance_interval,it,g) IMPLICIT NONE TYPE(POP_TYPE), INTENT(INOUT) :: POP INTEGER(i4b), INTENT(IN) :: disturbance_interval(:,:) INTEGER(i4b), INTENT(IN) :: it,g INTEGER(i4b) :: nage,iage, i_min, i_max, tmp_array(NPATCH2D) INTEGER(i4b) :: i_min_growth, i_max_growth REAL(dp) :: disturbance_freq,tmp_min,tmp_max, tmp1_min, tmp1_max REAL(dp) :: tmp2_min, tmp2_max REAL(dp) :: tmp3_min, tmp3_max REAL(dp) :: tmp4_min, tmp4_max LOGICAL :: MASK(NPATCH2D) INTEGER(i4b) :: age_min, age_max INTEGER(i4b) :: age_min_growth, age_max_growth INTEGER(i4b), ALLOCATABLE :: age(:) REAL(dp), ALLOCATABLE ::cmass_age(:), stress_mort_age(:), crowd_mort_age(:) REAL(dp), ALLOCATABLE ::csapwood_age(:), sapwood_area_age(:), growth_age(:) REAL(dp), ALLOCATABLE ::freq_age(:) ! get interpolated biomass,sapwood, stress mortality, crowding mortality, disturbance mortality POP%pop_grid(g)%cmass_sum= 0 POP%pop_grid(g)%stress_mortality = 0 POP%pop_grid(g)%cat_mortality = 0 pop%pop_grid(g)%crowding_mortality = 0 pop%pop_grid(g)%csapwood_sum = 0 pop%pop_grid(g)%sapwood_area = 0 tmp_array = 0 nage = MIN(POP%pop_grid(g)%patch(1)%disturbance_interval(1),it)+1 ! maximum age !nage = maxval(pop%pop_grid(g)%patch(:)%age(1)) disturbance_freq=1.0/REAL(disturbance_interval(g,1)) IF(.NOT.ALLOCATED(age)) ALLOCATE(age(nage)) IF(.NOT.ALLOCATED(freq_age)) ALLOCATE(freq_age(nage)) IF(.NOT.ALLOCATED(cmass_age)) ALLOCATE(cmass_age(nage)) IF(.NOT.ALLOCATED(growth_age)) ALLOCATE(growth_age(nage)) IF(.NOT.ALLOCATED(csapwood_age)) ALLOCATE(csapwood_age(nage)) IF(.NOT.ALLOCATED(sapwood_area_age)) ALLOCATE(sapwood_area_age(nage)) IF(.NOT.ALLOCATED(stress_mort_age)) ALLOCATE(stress_mort_age(nage)) IF(.NOT.ALLOCATED(crowd_mort_age)) ALLOCATE(crowd_mort_age(nage)) !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 tmp_max = 0.0 pop%pop_grid(g)%biomass_age = 0.0 !IF (pop%LU(g)==2) THEN ! secondary forest IF (POP%pop_grid(g)%LU==2) THEN ! secondary forest DO iage = 1, nage age(iage) = iage-1 freq_age(iage) = pop%pop_grid(g)%freq_age(iage) ENDDO ELSE DO iage = 1, nage age(iage) = iage-1 freq_age(iage) = REALExponential(disturbance_freq,REAL(age(iage),dp)) pop%pop_grid(g)%freq_age(iage) = freq_age(iage) END DO ENDIF IF (SUM(freq_age)>0.0) freq_age = freq_age/SUM(freq_age) 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 age_min_growth = age_min age_max_growth = age_max i_min_growth = i_min i_max_growth = i_max !write(*,*) 'I0',it, age(iage), age_min, age_max, i_min, i_max 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 ELSEWHERE tmp_array=0 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 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) ENDIF ! interpolate or extrapolate (growth) 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 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) 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) 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 ELSEWHERE tmp_array=0 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 tmp1_min = 0.0 tmp2_min = 0.0 tmp3_min = 0.0 tmp4_min = 0.0 ENDIF MASK = pop%pop_grid(g)%patch(:)%age(1).EQ.age_max WHERE (MASK) tmp_array= 1 ELSEWHERE tmp_array=0 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)* & REAL(age(iage)-age_min) stress_mort_age(iage) = tmp1_min + (tmp1_max-tmp1_min)/REAL(age_max-age_min)* & REAL(age(iage)-age_min) crowd_mort_age(iage) = tmp2_min + (tmp2_max-tmp2_min)/REAL(age_max-age_min)* & REAL(age(iage)-age_min) csapwood_age(iage) = tmp3_min + (tmp3_max-tmp3_min)/REAL(age_max-age_min)* & REAL(age(iage)-age_min) sapwood_area_age(iage) = tmp4_min + (tmp4_max-tmp4_min)/REAL(age_max-age_min)* & REAL(age(iage)-age_min) 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 ELSEWHERE tmp_array=0 endwhere tmp_min = SUM(pop%pop_grid(g)%patch(:)%growth,MASK)/SUM(tmp_array) ELSE tmp_min = 0.0 ENDIF MASK = pop%pop_grid(g)%patch(:)%age(1).EQ.age_max_growth WHERE (MASK) tmp_array= 1 ELSEWHERE tmp_array=0 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)* & REAL(age(iage)-age_min_growth) !cmass_age(iage) = cmass_age(iage) + growth_age(iage) - stress_mort_age(iage) - crowd_mort_age(iage) ENDIF ! write(*,*) 'age_min, age_max 1:', age(iage), age_min, age_max, i_min, i_max !if (it.eq.403) then ! witer(58,*) age(iage),age_min, age_max, tmp1_min, tmp1_max, stress_mort_age(iage), i_min !endif POP%pop_grid(g)%cmass_sum = POP%pop_grid(g)%cmass_sum + & freq_age(iage)*cmass_age(iage) POP%pop_grid(g)%csapwood_sum = POP%pop_grid(g)%csapwood_sum + & freq_age(iage)*csapwood_age(iage) POP%pop_grid(g)%sapwood_area = POP%pop_grid(g)%sapwood_area + & freq_age(iage)*sapwood_area_age(iage) !$if (g==2) then !$write(71, "(2i4, 350e16.6)") it, iage, freq_age(iage), cmass_age(iage), growth_age(iage), stress_mort_age(iage), & !$ crowd_mort_age(iage), tmp_min, tmp_max, real(age_max_growth), real(age_min_growth) !$endif !$if (g==2) write(72, "(2i4, 350e16.6)") it, iage, freq_age(iage), cmass_age(iage) !$if (g==1) write(71, "(2i4, 350e16.6)") it, iage, freq_age(iage), cmass_age(iage) POP%pop_grid(g)%stress_mortality = POP%pop_grid(g)%stress_mortality + & freq_age(iage)*stress_mort_age(iage) POP%pop_grid(g)%crowding_mortality = POP%pop_grid(g)%crowding_mortality + & freq_age(iage)*crowd_mort_age(iage) POP%pop_grid(g)%cat_mortality = POP%pop_grid(g)%growth - & POP%pop_grid(g)%stress_mortality - & POP%pop_grid(g)%crowding_mortality - & ( POP%pop_grid(g)%cmass_sum- POP%pop_grid(g)%cmass_sum_old) POP%pop_grid(g)%fire_mortality = 0.0 POP%pop_grid(g)%res_mortality = 0.0 pop%pop_grid(g)%biomass_age(iage) = cmass_age(iage) ENDDO !$if (it.gt.400) then !$ write(*,*) 'it, nage', it, nage !$ write(591, "(350e16.6)") freq_age !$ write(601,"(350e16.6)") cmass_age !$ write(602,"(350e16.6)") stress_mort_age !$ write(603,"(350e16.6)") real(age) !$endif DEALLOCATE(age) DEALLOCATE(freq_age) DEALLOCATE(cmass_age) DEALLOCATE(stress_mort_age) END SUBROUTINE INTERPOLATE_BIOMASS_1D !============================================================================= SUBROUTINE INTERPOLATE_BIOMASS_2D(pop, disturbance_interval,it,g) IMPLICIT NONE TYPE(POP_TYPE), INTENT(INOUT) :: POP INTEGER(i4b), INTENT(IN) :: disturbance_interval(:,:) INTEGER(i4b), INTENT(IN) :: it,g INTEGER(i4b), ALLOCATABLE :: A1(:), A2(:) ! interpolated ages INTEGER(i4b), ALLOCATABLE :: xobs(:), yobs(:) ! observed ages REAL(dp), ALLOCATABLE :: z1obs(:), z2obs(:), z3obs(:) ! observed biomass, stress_mort, crowd_mort REAL(dp), ALLOCATABLE :: z1interp(:), z2interp(:), z3interp(:) ! interpolated biomass, stress mortality, crowding mortality REAL(dp), ALLOCATABLE :: freq_interp(:) ! weightings for interpolated age pairs REAL(dp), ALLOCATABLE :: zp(:) ! euclidean distance from interpolated age pair to observed age pairs INTEGER(i4b) :: age_max(2), nrep(NPATCH2D+1) INTEGER(i4b) :: tmp1, tmp2, I1, I2,I3,I4 INTEGER(i4b) :: x1, x2, x3, x4, y1, y2, y3, y4 INTEGER(i4b) :: p, j, k, n, np, nobs, count_extrap,l, ct LOGICAL :: flag REAL(dp) :: biomass(NPATCH2D+1), stress_mort(NPATCH2D+1), crowd_mort(NPATCH2D+1) INTEGER(i4b) :: age1(NPATCH2D+1), age2(NPATCH2D+1) REAL(dp) :: zmin INTEGER(i4b), ALLOCATABLE :: interp_case(:), tmp_array(:), tmp(:) REAL(dp) :: area(4,4), x(4), y(4), disturbance_freq1, disturbance_freq2 INTEGER(i4b) :: triangle_points(4,3), I_inside_triangle, Ineighbour(8) LOGICAL :: MASK_INSIDE_TRIANGLE(4), IS_NEIGHBOUR(8), tmp_logical LOGICAL, ALLOCATABLE :: MASK2(:), MASK3(:), MASK4(:) INTEGER(i4b), ALLOCATABLE :: address(:,:) POP%pop_grid(g)%cmass_sum = 0.0 POP%pop_grid(g)%stress_mortality = 0.0 POP%pop_grid(g)%crowding_mortality = 0.0 ! Construct Age Interpolating Grid age_max(1) = MIN(POP%pop_grid(g)%patch(1)%disturbance_interval(1),it)+1 ! maximum age age_max(2) = MIN(POP%pop_grid(g)%patch(1)%disturbance_interval(NDISTURB),it)+1 ! maximum age p=0; DO j=0,age_max(1) DO k=0, age_max(2) IF (k>j) THEN p=p+1 ENDIF ENDDO ENDDO np = p ALLOCATE(A1(np)) ALLOCATE(A2(np)) ALLOCATE(z1interp(np)) ALLOCATE(z2interp(np)) ALLOCATE(z3interp(np)) ALLOCATE(freq_interp(np)) ALLOCATE(interp_case(np)) ALLOCATE(tmp_array(np)) ALLOCATE(tmp(np)) ALLOCATE(address(age_max(1)+1,age_max(2)+1)) p=0; address = -9999 DO j=0,age_max(1) DO k=0, age_max(2) IF (k>j) THEN p=p+1 A1(p) = j A2(p) = k address(j+1,k+1) = p ENDIF ENDDO ENDDO ! Construct Age observations nrep = 1 flag = .FALSE. j = 1 ! create "observed" age pair (0,0) with zero biomass stress_mort(j) = 0.0 biomass(j) = 0.0 crowd_mort(j) = 0.0 age1(j) = 0 age2(j) = 0 j = j+1 DO k=1,NPATCH2D tmp1 = pop%pop_grid(g)%patch(k)%age(1) tmp2 = pop%pop_grid(g)%patch(k)%age(NDISTURB) IF (j.GT.1) THEN DO n=1,j-1 IF (tmp1 == age1(n) .AND. tmp2==age2(n)) THEN flag = .TRUE. nrep(n) = nrep(n) + 1 biomass(n) = biomass(n) + pop%pop_grid(g)%patch(k)%layer(1)%biomass stress_mort(n) = stress_mort(n) + pop%pop_grid(g)%patch(k)%stress_mortality crowd_mort(n) = crowd_mort(n) + pop%pop_grid(g)%patch(k)%crowding_mortality ENDIF ENDDO ENDIF IF (flag .EQV. .FALSE.) THEN age1(j) = tmp1 age2(j) = tmp2 biomass(j) = pop%pop_grid(g)%patch(k)%layer(1)%biomass stress_mort(j) = pop%pop_grid(g)%patch(k)%stress_mortality crowd_mort(j) = pop%pop_grid(g)%patch(k)%crowding_mortality IF (k.NE.NPATCH2D) j = j+1 ELSE flag = .FALSE. ENDIF ENDDO nobs = j biomass(1:nobs) = biomass(1:nobs)/nrep(1:nobs) stress_mort(1:nobs) = stress_mort(1:nobs)/nrep(1:nobs) crowd_mort(1:nobs) = crowd_mort(1:nobs)/nrep(1:nobs) ALLOCATE(xobs(nobs)) ALLOCATE(yobs(nobs)) ALLOCATE(z1obs(nobs)) ALLOCATE(z2obs(nobs)) ALLOCATE(z3obs(nobs)) ALLOCATE(zp(nobs)) ALLOCATE(MASK2(nobs)) ALLOCATE(MASK3(nobs)) ALLOCATE(MASK4(nobs)) xobs = age1(1:nobs) yobs = age2(1:nobs) z1obs = biomass(1:nobs) z2obs = stress_mort(1:nobs) z3obs = crowd_mort(1:nobs) IF (it.EQ.500) THEN DO k=1,nobs WRITE(502, "(2i5,2e16.6,i5)") xobs(k), yobs(k), z1obs(k), z3obs(k), nrep(k) ENDDO ENDIF ! get weightings for each interpolated age pair DO k = 1, np disturbance_freq1=1.0/REAL(disturbance_interval(g,1)) disturbance_freq2=1.0/REAL(disturbance_interval(g,2)) freq_interp(k) = REALExponential(disturbance_freq1,REAL(A1(k),dp)) * & REALExponential(disturbance_freq2,REAL(A2(k),dp)) ENDDO freq_interp = freq_interp/SUM(freq_interp) ! interpolate DO p=1,np ! loop over interpolated age pairs ! get distance to all observations DO j=1,nobs zp(j) = SQRT((REAL(A1(p))-REAL(xobs(j)))**2+(REAL(A2(p))-REAL(yobs(j)))**2) ENDDO ! get closest point zmin = MINVAL(zp) I1 = MINLOC(zp,1) x1 = xobs(I1) y1 = yobs(I1) ! check for obs locations forming a quadrangle around interpolating point MASK2 = (SIGN(1,A1(p)-xobs)== -SIGN(1,A1(p)-x1)).AND.(SIGN(1,A2(p)-yobs)== SIGN(1,A2(p)-y1).AND.A1(p).NE.xobs) MASK3 = (SIGN(1,A1(p)-xobs)== SIGN(1,A1(p)-x1)).AND.(SIGN(1,A2(p)-yobs)== -SIGN(1,A2(p)-y1).AND.A2(p).NE.yobs) MASK4 = (SIGN(1,A1(p)-xobs)== -SIGN(1,A1(p)-x1)).AND.(SIGN(1,A2(p)-yobs)== -SIGN(1,A2(p)-y1) & .AND.A1(p).NE.xobs.AND.A2(p).NE.yobs) IF ((ANY(MASK2)).AND.(ANY(MASK3)).AND.(ANY(MASK4))) THEN ! get nearest point with opposing sign of x displacement I2 = MINLOC(zp,1, MASK2) x2 = xobs(I2) y2 = yobs(I2) ! get nearest point with opposing sign of y displacement I3 = MINLOC(zp,1, MASK3) x3 = xobs(I3) y3 = yobs(I3) ! get nearest point with opposing sign of x & y displacements I4 = MINLOC(zp,1, MASK4) x4 = xobs(I4) y4 = yobs(I4) tmp_logical = .NOT.( (x2.EQ.0.AND.y2.EQ.0).OR.(x3.EQ.0.AND.y3.EQ.0).OR.(x4.EQ.0.AND.y4.EQ.0)) ENDIF IF ((A1(p)==x1).AND.(A2(p)==y1)) THEN interp_case(p)=1 ELSEIF ((ANY(MASK2)).AND.(ANY(MASK3)).AND.(ANY(MASK4)) .AND.tmp_logical) THEN ! quadrangle (without (0,0)) exists interp_case(p) = 2 ELSE interp_case(p) = 3 ENDIF!j=1 SELECT CASE (interp_case(p)) CASE(1) ! interpolated point is the same as observation z1interp(p) = z1obs(I1) z2interp(p) = z2obs(I1) z3interp(p) = z3obs(I1) CASE(3) ! extrapolation required: set to value of nearest observation z1interp(p) = z1obs(I1) z2interp(p) = z2obs(I1) z3interp(p) = z3obs(I1) CASE(2) ! quadrangle ! get nearest point with opposing sign of x displacement I2 = MINLOC(zp,1, MASK2) x2 = xobs(I2) y2 = yobs(I2) ! get nearest point with opposing sign of y displacement I3 = MINLOC(zp,1, MASK3) x3 = xobs(I3) y3 = yobs(I3) ! get nearest point with opposing sign of x & y displacements I4 = MINLOC(zp,1, MASK4) x4 = xobs(I4) y4 = yobs(I4) x=REAL((/x1, x2, x3, x4/),dp) y =REAL((/y1, y2, y3, y4/),dp) ! get area of four possible triangles from 4 vertices, and corresponding 3 partial ! triangles, with observation as one vertex area = 0.0 triangle_points = 0 triangle_points(1,:) = (/I1, I2, I3/) area(1,1) = area_triangle(x(1), y(1),x(2), y(2),x(3), y(3)) area(1,2) = area_triangle(REAL(A1(p),dp), REAL(A2(p),dp), x(2), y(2), x(3), y(3)) area(1,3) = area_triangle(REAL(A1(p),dp), REAL(A2(p),dp), x(1), y(1), x(3), y(3)) area(1,4) = area_triangle(REAL(A1(p),dp), REAL(A2(p),dp), x(1), y(1), x(2), y(2)) triangle_points(2,:) = (/I1, I2, I4/) area(2,1) = area_triangle(x(1), y(1), x(2), y(2), x(4), y(4)) area(2,2) = area_triangle(REAL(A1(p),dp), REAL(A2(p),dp), x(2), y(2), x(4), y(4)) area(2,3) = area_triangle(REAL(A1(p),dp), REAL(A2(p),dp), x(1), y(1), x(4), y(4)) area(2,4) = area_triangle(REAL(A1(p),dp), REAL(A2(p),dp), x(1), y(1), x(2), y(2)) triangle_points(3,:) = (/I1, I4, I3/) area(3,1) = area_triangle(x(1) , y(1), x(3), y(3),x(4), y(4)) area(3,2) = area_triangle(REAL(A1(p),dp), REAL(A2(p),dp), x(4), y(4), x(3), y(3)) area(3,3) = area_triangle(REAL(A1(p),dp), REAL(A2(p),dp), x(1), y(1), x(3), y(3)) area(3,4) = area_triangle(REAL(A1(p),dp), REAL(A2(p),dp), x(1), y(1), x(4), y(4)) triangle_points(4,:) = (/I2, I4, I3/) area(4,1) = area_triangle(x(2), y(2), x(3), y(3), x(4), y(4)); area(4,2) = area_triangle(REAL(A1(p),dp), REAL(A2(p),dp), x(4), y(4), x(3), y(3)) area(4,3) = area_triangle(REAL(A1(p),dp), REAL(A2(p),dp), x(2), y(2), x(3), y(3)) area(4,4) = area_triangle(REAL(A1(p),dp), REAL(A2(p),dp), x(2), y(2), x(4), y(4)) MASK_INSIDE_TRIANGLE=(area(:,1)==SUM(area(:,2:4),2)); I_inside_triangle = MINLOC(area(:,1),1,MASK_INSIDE_TRIANGLE) z1interp(p) =SUM( z1obs(triangle_points(I_inside_triangle,:))* & area(I_inside_triangle,2:4))/SUM(area(I_inside_triangle,2:4)) z2interp(p) =SUM( z2obs(triangle_points(I_inside_triangle,:))* & area(I_inside_triangle,2:4))/SUM(area(I_inside_triangle,2:4)) z3interp(p) =SUM( z3obs(triangle_points(I_inside_triangle,:))* & area(I_inside_triangle,2:4))/SUM(area(I_inside_triangle,2:4)) CASE default WRITE(*,*) " illegal interpolation case." STOP END SELECT ! interpolation case ENDDO ! loop over interpolated age pairs p = 0 ! counter for while loop count_extrap = 0 ! Extrapolation IF (ANY(interp_case.NE.3)) THEN DO WHILE (ANY(interp_case==3)) DO ct = 3,3 DO p=1,np ! find number of neighbours for each extrapolable point is_neighbour = .FALSE. Ineighbour = 1 tmp(p) = 0 IF (interp_case(p)==3) THEN DO k=1,8 IF(k==1) THEN tmp1=A1(p)-1 tmp2=A2(p) ELSEIF(k==2) THEN tmp1 = A1(p)+1 tmp2 = A2(p) ELSEIF (k==3) THEN tmp1 = A1(p) tmp2 = A2(p)-1 ELSEIF (k==4) THEN tmp1 = A1(p) tmp2 = A2(p)+1 ELSEIF (k==5) THEN tmp1=A1(p)-1 tmp2=A2(p)-1 ELSEIF (k==6) THEN tmp1=A1(p)+1 tmp2=A2(p)+1 ELSEIF (k==7) THEN tmp1=A1(p)+1 tmp2=A2(p)-1 ELSEIF (k==8) THEN tmp1=A1(p)-1 tmp2=A2(p)+1 ENDIF IF(tmp1.GE.0.AND.tmp1.LE.age_max(1).AND.tmp2.GE.0.AND.tmp2.LE.age_max(2).AND.tmp2.GT.tmp1) THEN is_neighbour(k) = (interp_case(address(tmp1+1,tmp2+1)).NE.3) Ineighbour(k) = address(tmp1+1,tmp2+1) ENDIF ENDDO tmp(p) = COUNT(is_neighbour) IF (tmp(p).GE.ct) THEN z1interp(p) = SUM(z1interp(Ineighbour) ,is_neighbour)/REAL(tmp(p)) z2interp(p) = SUM(z2interp(Ineighbour) ,is_neighbour)/REAL(tmp(p)) z3interp(p) = SUM(z3interp(Ineighbour) ,is_neighbour)/REAL(tmp(p)) interp_case(p) = 2 ENDIF ENDIF ENDDO !1, np ENDDO IF (ALL(tmp.LT.3)) THEN ! maximum of one neighbour p = MAXLOC(tmp,1) is_neighbour = .FALSE. Ineighbour = 1 tmp(p) = 0 DO k=1,8 IF(k==1) THEN tmp1=A1(p)-1 tmp2=A2(p) ELSEIF(k==2) THEN tmp1 = A1(p)+1 tmp2 = A2(p) ELSEIF (k==3) THEN tmp1 = A1(p) tmp2 = A2(p)-1 ELSEIF (k==4) THEN tmp1 = A1(p) tmp2 = A2(p)+1 ELSEIF (k==5) THEN tmp1=A1(p)-1 tmp2=A2(p)-1 ELSEIF (k==6) THEN tmp1=A1(p)+1 tmp2=A2(p)+1 ELSEIF (k==7) THEN tmp1=A1(p)+1 tmp2=A2(p)-1 ELSEIF (k==8) THEN tmp1=A1(p)-1 tmp2=A2(p)+1 ENDIF IF(tmp1.GE.0.AND.tmp1.LE.age_max(1).AND.tmp2.GE.0.AND.tmp2.LE.age_max(2).AND.tmp2.GT.tmp1) THEN is_neighbour(k) = (interp_case(address(tmp1+1,tmp2+1)).NE.3) Ineighbour(k) = address(tmp1+1,tmp2+1) ENDIF ENDDO tmp(p) = COUNT(is_neighbour) z1interp(p) = SUM(z1interp(Ineighbour) ,is_neighbour)/REAL(tmp(p)) z2interp(p) = SUM(z2interp(Ineighbour) ,is_neighbour)/REAL(tmp(p)) z3interp(p) = SUM(z3interp(Ineighbour) ,is_neighbour)/REAL(tmp(p)) interp_case(p) = 2 ENDIF ! ALL(tmp.lt.2) ENDDO ! do while ENDIF DO p = 1,np POP%pop_grid(g)%cmass_sum = POP%pop_grid(g)%cmass_sum + & freq_interp(p)*z1interp(p) POP%pop_grid(g)%stress_mortality = POP%pop_grid(g)%stress_mortality + & freq_interp(p)*z2interp(p) POP%pop_grid(g)%crowding_mortality = POP%pop_grid(g)%crowding_mortality + & freq_interp(p)*z3interp(p) ENDDO POP%pop_grid(g)%cat_mortality = POP%pop_grid(g)%cmass_sum * disturbance_freq2 POP%pop_grid(g)%fire_mortality = POP%pop_grid(g)%growth - & POP%pop_grid(g)%cat_mortality - & POP%pop_grid(g)%stress_mortality - & POP%pop_grid(g)%crowding_mortality - & ( POP%pop_grid(g)%cmass_sum- POP%pop_grid(g)%cmass_sum_old) DEALLOCATE(xobs) DEALLOCATE(yobs) DEALLOCATE(z1obs) DEALLOCATE(z2obs) DEALLOCATE(z3obs) DEALLOCATE(zp) DEALLOCATE(MASK2) DEALLOCATE(MASK3) DEALLOCATE(MASK4) DEALLOCATE(A1) DEALLOCATE(A2) DEALLOCATE(z1interp) DEALLOCATE(z2interp) DEALLOCATE(z3interp) DEALLOCATE(freq_interp) DEALLOCATE(interp_case) DEALLOCATE(tmp) DEALLOCATE(tmp_array) DEALLOCATE(address) END SUBROUTINE INTERPOLATE_BIOMASS_2D !============================================================================= SUBROUTINE SMOOTH_FLUX(POP,g,t) IMPLICIT NONE TYPE(POP_TYPE), INTENT(INOUT) :: POP INTEGER(i4b), INTENT(IN) :: g, t INTEGER(i4b), PARAMETER :: SPAN = NYEAR_SMOOTH/2 REAL(dp) :: x(SPAN+1), y(SPAN+1), 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 IF (t.GT.NYEAR_HISTORY) THEN DO k = 1, NYEAR_HISTORY-1 POP%pop_grid(g)%fire_mortality_history(k) = POP%pop_grid(g)%fire_mortality_history(k+1) ENDDO ENDIF POP%pop_grid(g)%fire_mortality_history(t) = POP%pop_grid(g)%fire_mortality flux = POP%pop_grid(g)%fire_mortality_history t0 = t-SPAN n = 0 sumflux = 0.0 sumsmooth = 0.0 DO tt = 1,NYEAR_SMOOTH IF ((t0+tt).GE.1 .AND. (t0+tt).LE.t+1) THEN sumflux = sumflux + flux(t0+tt) y(tt) = flux(t0+tt) x(tt) = tt n = n+1 IF ((t0+tt).EQ.t+1) THEN CALL regress(x,y,n,a,b,r) ENDIF ELSE sumflux = sumflux + (a + b*tt) n = n+ 1 ENDIF ENDDO dbuf =POP%pop_grid(g)%smoothing_buffer/(REAL(NYEAR_SMOOTH)/2.0) !MC smoothed_flux=max(sumflux/real(n)+dbuf, 0.0) smoothed_flux=MAX(sumflux/REAL(n)+dbuf, 0.0_dp) POP%pop_grid(g)%smoothing_buffer = POP%pop_grid(g)%smoothing_buffer + flux(t) - smoothed_flux POP%pop_grid(g)%fire_mortality_smoothed = smoothed_flux END SUBROUTINE SMOOTH_FLUX !============================================================================= SUBROUTINE SMOOTH_FLUX_cat(POP,g,t) IMPLICIT NONE TYPE(POP_TYPE), INTENT(INOUT) :: POP INTEGER(i4b), INTENT(IN) :: g, t INTEGER(i4b), PARAMETER :: SPAN = NYEAR_SMOOTH/2 REAL(dp) :: x(SPAN+1), y(SPAN+1), 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 IF (t.GT.NYEAR_HISTORY) THEN DO k = 1, NYEAR_HISTORY-1 POP%pop_grid(g)%cat_mortality_history(k) = POP%pop_grid(g)%cat_mortality_history(k+1) ENDDO ENDIF POP%pop_grid(g)%cat_mortality_history(t) = POP%pop_grid(g)%cat_mortality flux = POP%pop_grid(g)%cat_mortality_history t0 = t-SPAN n = 0 sumflux = 0.0 sumsmooth = 0.0 DO tt = 1,NYEAR_SMOOTH IF ((t0+tt).GE.1 .AND. (t0+tt).LE.t+1) THEN sumflux = sumflux + flux(t0+tt-1) y(tt) = flux(t0+tt-1) x(tt) = tt n = n+1 IF ((t0+tt).EQ.t+1) THEN CALL regress(x,y,n,a,b,r) ENDIF ELSE sumflux = sumflux + (a + b*tt) n = n+ 1 ENDIF ENDDO dbuf =POP%pop_grid(g)%smoothing_buffer_cat/(REAL(NYEAR_SMOOTH)/2.0) !MC smoothed_flux=max(sumflux/real(n)+dbuf, 0.0) smoothed_flux=MAX(sumflux/REAL(n)+dbuf, 0.0_dp) POP%pop_grid(g)%smoothing_buffer_cat = POP%pop_grid(g)%smoothing_buffer_cat + flux(t) - smoothed_flux POP%pop_grid(g)%cat_mortality_smoothed = smoothed_flux END SUBROUTINE SMOOTH_FLUX_cat !============================================================================= SUBROUTINE REGRESS(x, y, n, a, b, r) IMPLICIT NONE REAL(dp), INTENT(IN) :: x(:), y(:) REAL(dp), INTENT(OUT) :: a, b, r INTEGER(i4b), INTENT(IN) :: n REAL(dp):: sx,sy,sxx,sxy,delta,meanx,meany,sdx,sdy INTEGER(i4b) :: i !Performs a linear regression of array y on array x (n values) ! returning parameters a and b in the fitted model: y=a+bx ! Source: Press et al 1986, Sect 14.2 ! also returns Pearson r sx=0.0 sy=0.0 sxx=0.0 sxy=0.0 DO i=1,n sx = sx + x(i) sy = sy + y(i) sxx = sxx + x(i) * y(i) sxy = sxy + x(i)*y(i) ENDDO delta = REAL(n)*sxx - sx*sx a=(sxx*sy-sx*