Smoth growth

R/multiSitePrebas.r

@@ -10,8 +10,9 @@ InitMultiSite <- function(nYearsMS,
                           multiThin = NA,
                           multiNthin = NA,
                           multiInitClearCut = NA,
-                          fixBAinitClarcut = 1.,
-                          initCLcutRatio = NA,
+                          fixBAinitClarcut = 1.,  ###if 1 when clearcut occur the species inital biomass is fixed at replanting using the values in initCLcutRatio else at replanting the replanting follows species relBa at last year 
+                          initCLcutRatio = NA,  ###BA ratio per each species/layer (default is the ba ratio at the begginning of the simulations)
+                          areas = NA,
                           PAR,
                           TAir,
                           VPD,
@@ -27,19 +28,18 @@ InitMultiSite <- function(nYearsMS,
                           inDclct = NA,
                           inAclct = NA,
                           yassoRun = 0,
-                          lukeRuns){
+                          lukeRuns,
+                          smoothP0 = 1,
+                          smoothETS = 1,
+                          smoothYear=5){
 
   nSites <- length(nYearsMS)
+  if(all(is.na(areas))) areas <- rep(1.,nSites) ###each site is 1 ha (used to scale regional harvest)
   if(all(is.na(siteInfo))){
     siteInfo = matrix(c(1,1,3,160,0,0,20,3,3),nSites,9,byrow = T) ###default values for nspecies and site type = 3
     siteInfo[,1] <- 1:nSites
   }
   nLayers <- siteInfo[,8]
-  if(length(fixBAinitClarcut)==1) fixBAinitClarcut=rep(fixBAinitClarcut,nSites)
-  if(all(is.na(initCLcutRatio))){
-    initCLcutRatio <- matrix(0.,nSites,max(nLayers))
-    for(iz in 1:nSites) initCLcutRatio[iz,1:nLayers[iz]] <- rep(1/nLayers[iz],nLayers[iz])
-  }
   # nSp <- siteInfo[,9]
   climIDs <- siteInfo[,2]
   # if(all(is.na(multiInitVar)) & all(is.na(nSp)) nSp <- rep(3,nSites)
@@ -96,7 +96,10 @@ InitMultiSite <- function(nYearsMS,
                                                 which(is.na(multiInitClearCut[sitesClimID,5])),round(Ainit))
   }
   ETSthres <- 1000; ETSmean <- rowMeans(multiETS)
-
+  if(smoothETS==1. & maxYears > 1){
+    for(i in 2:maxYears) multiETS[,i] <- multiETS[,(i-1)] + (multiETS[,i]-multiETS[,(i-1)])/min(i,smoothYear)
+  } 
+
   ####process clearcut
   for(i in 1: nSites){
     if(ClCut[i]==1 & all(is.na(inDclct[i,]))) inDclct[i,] <-
@@ -140,16 +143,24 @@ InitMultiSite <- function(nYearsMS,
   ### compute P0
   ###if P0 is not provided use preles to compute P0
   if(all(is.na(multiP0))){
-    multiP0 <- matrix(NA,nClimID,maxYears)
+    multiP0 <- array(NA,dim=c(nClimID,maxYears,2))
     for(climID in 1:nClimID){
       nYearsX <- max(nYearsMS[which(climIDs==climID)])
       P0 <- PRELES(DOY=rep(1:365,nYearsX),PAR=PAR[climID,1:(365*nYearsX)],
                    TAir=TAir[climID,1:(365*nYearsX)],VPD=VPD[climID,1:(365*nYearsX)],
                    Precip=Precip[climID,1:(365*nYearsX)],CO2=rep(380,(365*nYearsX)),
                    fAPAR=rep(1,(365*nYearsX)),LOGFLAG=0,p=pPRELES)$GPP
       P0 <- matrix(P0,365,nYearsX)
-      multiP0[climID,(1:nYearsX)] <- colSums(P0)
-    }}
+      multiP0[climID,(1:nYearsX),1] <- colSums(P0)
+    }
+    if(smoothP0==1 & maxYears > 1){
+      multiP0[,1,2] <- multiP0[,1,1]
+      for(i in 2:maxYears) multiP0[,i,2] <- multiP0[,(i-1),2] + (multiP0[,i,1]-multiP0[,(i-1),2])/min(i,smoothYear)
+      # multiP0[,,2] <- matrix(rowMeans(multiP0[,,1]),nClimID,maxYears,byrow = F)
+    } else{
+    multiP0[,,2] <- multiP0[,,1]
+    }
+  }
 
   if (all(is.na(multiInitVar))){
     multiInitVar <- array(NA,dim=c(nSites,6,maxNlayers))
@@ -158,6 +169,12 @@ InitMultiSite <- function(nYearsMS,
     multiInitVar[,5,] <- initClearcut[3]/maxNlayers; multiInitVar[,6,] <- initClearcut[4]
     multiInitVar[,2,] <- matrix(multiInitClearCut[,5],nSites,maxNlayers)
   }
+  if(length(fixBAinitClarcut)==1) fixBAinitClarcut=rep(fixBAinitClarcut,nSites)
+  if(all(is.na(initCLcutRatio))){
+    initCLcutRatio <- multiInitVar[,5,]/rowSums(multiInitVar[,5,])
+    # initCLcutRatio <- matrix(0.,nSites,max(nLayers))
+    # for(iz in 1:nSites) initCLcutRatio[iz,1:nLayers[iz]] <- rep(1/nLayers[iz],nLayers[iz])
+  }
   if(all(is.na(litterSize))){
     litterSize <- matrix(0,3,allSp)
     litterSize[2,] <- 2
@@ -173,6 +190,7 @@ InitMultiSite <- function(nYearsMS,
     maxYears = maxYears,
     maxThin = maxThin,
     nYears = nYearsMS,
+    areas = areas,
     thinning = multiThin,
     pCROBAS = pCROBAS,
     allSp = allSp,
@@ -203,7 +221,9 @@ InitMultiSite <- function(nYearsMS,
     dailyPRELES = array(-999,dim=c(nSites,(maxYears*365),3)),
     yassoRun = yassoRun,
     lukeRuns = lukeRuns,
-    PREBASversion = PREBASversion)
+    PREBASversion = PREBASversion,
+    smoothP0 = smoothP0,
+    smoothETS = smoothETS)
   return(multiSiteInit)
 }
 
@@ -227,7 +247,7 @@ multiPrebas <- function(multiSiteInit){
                      fixBAinitClearcut = as.double(multiSiteInit$fixBAinitClarcut),
                      initCLcutRatio = as.matrix(multiSiteInit$initCLcutRatio),
                      ETSy=as.matrix(multiSiteInit$ETSy),
-                     P0y=as.matrix(multiSiteInit$P0y),
+                     P0y=as.array(multiSiteInit$P0y),
                      multiInitVar=as.array(multiSiteInit$multiInitVar),
                      weather=as.array(multiSiteInit$weather),
                      DOY= as.integer(multiSiteInit$DOY),
@@ -268,6 +288,7 @@ regionPrebas <- function(multiSiteInit,
                    minDharv = as.double(minDharv),
                    multiOut = as.array(multiSiteInit$multiOut),
                    nSites = as.integer(multiSiteInit$nSites),
+                   areas = as.double(multiSiteInit$areas),
                    nClimID = as.integer(multiSiteInit$nClimID),
                    nLayers = as.integer(multiSiteInit$nLayers),######
                    maxYears = as.integer(multiSiteInit$maxYears),
@@ -284,7 +305,7 @@ regionPrebas <- function(multiSiteInit,
                    fixBAinitClarcut = as.double(multiSiteInit$fixBAinitClarcut),
                    initCLcutRatio = as.matrix(multiSiteInit$initCLcutRatio),
                    ETSy=as.matrix(multiSiteInit$ETSy),
-                   P0y=as.matrix(multiSiteInit$P0y),
+                   P0y=as.array(multiSiteInit$P0y),
                    multiInitVar=as.array(multiSiteInit$multiInitVar),
                    weather=as.array(multiSiteInit$weather),
                    DOY= as.integer(multiSiteInit$DOY),
@@ -306,9 +327,10 @@ regionPrebas <- function(multiSiteInit,
                    lukeRuns=as.double(multiSiteInit$lukeRuns))
 class(prebas) <- "regionPrebas"
 if(prebas$maxNlayers>1){
-    prebas$totHarv <- apply(prebas$multiOut[,,37,,1],2,sum)
+    rescalVbyArea <- prebas$multiOut[,,37,,1] * prebas$areas
+    prebas$totHarv <- apply(rescalVbyArea,2,sum)
   }else{
-    prebas$totHarv <- prebas$multiOut[,,37,,1]
+    prebas$totHarv <- colSums(prebas$multiOut[,,37,1,1]*prebas$areas)
   }
 return(prebas)
 }

R/prebas.r

@@ -22,7 +22,10 @@ prebas <- function(nYears,
                    ClCut = 1.,
                    inDclct = NA,
                    inAclct = NA,
-                   yassoRun = 0){
+                   yassoRun = 0,
+                   smoothP0 = 1,
+                   smoothETS = 1,
+                   smoothYear=5){
 
   ###process weather###
   if(length(PAR) >= (nYears*365)){
@@ -65,6 +68,7 @@ prebas <- function(nYears,
   Temp <- TAir[1:(365*nYears)]-5
   ETS <- pmax(0,Temp,na.rm=T)
   ETS <- matrix(ETS,365,nYears); ETS <- colSums(ETS)
+  if(smoothETS==1.) for(i in 2:nYears) ETS[i] <- ETS[(i-1)] + (ETS[i]-ETS[(i-1)])/min(i,smoothYear)
 
   ###if P0 is not provided use preles to compute P0
   if(is.na(P0)){
@@ -73,6 +77,11 @@ prebas <- function(nYears,
                  fAPAR=rep(1,length(PAR)),LOGFLAG=0,p=pPRELES)$GPP
     P0 <- matrix(P0,365,nYears);P0 <- colSums(P0)
   }
+  P0 <- matrix(P0,nYears,2)
+  if(smoothP0==1.){
+    P0[1,2] <- P0[1,1]
+    for(i in 2:nYears) P0[i,2] <- P0[(i-1),2] + (P0[i,1]-P0[(i-1),2])/min(i,smoothYear)
+  } 
 
   ETSthres <- 1000; ETSmean <- mean(ETS)
 
@@ -147,7 +156,7 @@ prebas <- function(nYears,
                      fixBAinitClarcut=as.numeric(fixBAinitClarcut),
                      initCLcutRatio = as.double(initCLcutRatio),
                      ETS = as.numeric(ETS),
-                     P0 = as.numeric(P0),
+                     P0 = as.matrix(P0),
                      weather=as.array(weatherPreles),
                      DOY= as.integer(1:365),
                      pPRELES=as.numeric(pPRELES),

src/B_prebas_v0.f90

@@ -27,7 +27,7 @@ subroutine prebas_v0(nYears,nLayers,nSp,siteInfo,pCrobas,initVar,thinning,output
  integer, intent(inout) :: nThinning
  real (kind=8), intent(out) :: fAPAR(nYears)
  real (kind=8), intent(inout) :: dailyPRELES((nYears*365),3)
- real (kind=8), intent(in) :: initVar(6,nLayers),P0y(nYears),ETSy(nYears),initCLcutRatio(nLayers)!
+ real (kind=8), intent(in) :: initVar(6,nLayers),P0y(nYears,2),ETSy(nYears),initCLcutRatio(nLayers)!
  real (kind=8), intent(inout) :: siteInfo(7)
  real (kind=8), intent(out) :: output(nYears,nVar,nLayers,2)
  real (kind=8), intent(inout) :: soilCinOut(nYears,5,3,nLayers),soilCtotInOut(nYears) !dimensions = nyears,AWENH,treeOrgans(woody,fineWoody,Foliage),species
@@ -41,7 +41,7 @@ subroutine prebas_v0(nYears,nLayers,nSp,siteInfo,pCrobas,initVar,thinning,output
  real (kind=8) :: STAND(nVar),STAND_tot(nVar),param(npar)!, output(nYear,nSites,nVar)
  integer :: i, ij, ijj,species,layer,nSpec,ll! tree species 1,2,3 = scots pine, norway spruce, birch
 
- real (kind=8) :: p0_ref, ETS_ref
+ real (kind=8) :: p0_ref, ETS_ref,P0yX(nYears,2)
  integer :: time, ki, year,yearX,Ainit, countThinning,domSp(1)
  real (kind=8) :: step, totBA
 
@@ -69,7 +69,7 @@ subroutine prebas_v0(nYears,nLayers,nSp,siteInfo,pCrobas,initVar,thinning,output
  real (kind=8) :: hb, A, B2,beta0, beta1,beta2, betas, betab
  real (kind=8) :: c,dHc,dH,dLc,g0,g1,g2,g3,g4,g5
  real (kind=8) :: npp, p_eff_all
- real (kind=8) :: p_eff, par_alfar,p
+ real (kind=8) :: p_eff, par_alfar,p,gpp_sp
  real (kind=8) :: s0,par_s0scale
  real (kind=8) :: weight, dNp,dNb,dNs
  real (kind=8) :: W_wsap, respi_m, respi_tot, V_scrown, V_bole, V,Vold
@@ -79,12 +79,9 @@ subroutine prebas_v0(nYears,nLayers,nSp,siteInfo,pCrobas,initVar,thinning,output
 
 !fix parameters
  real (kind=8) :: qcTOT0,Atot,fAPARprel(365)
-
+!v1 version definitions
  real (kind=8) :: theta
 
-
- ! open(2,file="test.txt")
- ! write(2,*) "site = ",siteInfo(1)
 !###initialize model###!
 fbAWENH = 0.
 folAWENH = 0.
@@ -99,14 +96,14 @@ subroutine prebas_v0(nYears,nLayers,nSp,siteInfo,pCrobas,initVar,thinning,output
 pars(25) = siteInfo(5)!CWinit
 pars(26) = siteInfo(6) !SOGinit
 pars(27) = siteInfo(7) !Sinit
-
+P0yX = P0y
 
  do i = 1,nLayers
   modOut(:,4,i,1) = initVar(1,i)  ! assign species
   modOut(:,7,i,1) = initVar(2,i) ! assign initAge !age can be made species specific assigning different ages to different species
   modOut(1,39,i,1) = sum(soilC(1,:,:,i)) !assign initial soilC
-  modOut(:,5,i,1) = ETSy	! assign ETS
-  modOut(:,6,i,1) = P0y		! assign P0
+  modOut(:,5,i,1) = ETSy! assign ETS
+  modOut(:,6,i,1) = P0yX(:,2)	! assign P0
  enddo
  modOut(:,1,:,1) = siteInfo(1); modOut(:,2,:,1) = siteInfo(2)	!! assign siteID and climID
  modOut(1,11,:,1) = initVar(3,:)
@@ -357,8 +354,8 @@ subroutine prebas_v0(nYears,nLayers,nSp,siteInfo,pCrobas,initVar,thinning,output
    pars(26) = prelesOut(4); siteInfo(6) = prelesOut(4) !SOGinit
    pars(27) = prelesOut(14); siteInfo(7) = prelesOut(14) !Sinit
 
-   STAND_all(10,:) = prelesOut(1)/1000.! Photosynthesis in g C m-2 (converted to kg C m-2)
-		endif
+   STAND_all(10,:) = prelesOut(1)/1000. ! Photosynthesis in g C m-2 (converted to kg C m-2)
+
 endif
 !enddo !! end site loop
 
@@ -426,7 +423,7 @@ subroutine prebas_v0(nYears,nLayers,nSp,siteInfo,pCrobas,initVar,thinning,output
   leff = STAND(19)
   keff = STAND(20)
   lproj = STAND(21)
-  p_eff_all = STAND(10) !!##!!2
+  p_eff_all = STAND(10)*P0yX(year,2)/P0yX(year,1) !!##!!2    smoothing PHOTOSYNTHESIS
   weight = STAND(23)
 
   rc = Lc / (H-1.3) !crown ratio
@@ -456,7 +453,8 @@ subroutine prebas_v0(nYears,nLayers,nSp,siteInfo,pCrobas,initVar,thinning,output
    par_alfar = par_alfar5
   end if
 
-!relate metabolic and structural parameters to site conditions
+!!relate metabolic and structural parameters to site conditions
+
 !  par_mf = par_mf0 * p0 / p0_ref
 !  par_mr = par_mr0 * p0 / p0_ref
 !  par_mw = par_mw0 * p0 / p0_ref
@@ -477,6 +475,7 @@ subroutine prebas_v0(nYears,nLayers,nSp,siteInfo,pCrobas,initVar,thinning,output
     !GPP all STAND$species   UNITS: g C  /  m2
     ! -------------------------------------
         p_eff = weight * p_eff_all
+		gpp_sp = weight * STAND(10)
 
     if(wf_STKG > 0.) then
         s0 = min(par_s0scale * P0 * par_k * par_sla, P_eff / wf_STKG * 10000.)
@@ -595,8 +594,10 @@ subroutine prebas_v0(nYears,nLayers,nSp,siteInfo,pCrobas,initVar,thinning,output
       W_wsap = N * par_rhow * A * (beta1 * H + beta2 * Hc)
       Respi_m = (par_mf + par_alfar*par_mr)* wf_STKG + par_mw * W_wsap
 ! note changes in the equations below AM 15.5.2015
-      npp = (weight * p_eff_all - Respi_m / 10000.) / (1.+par_c)
-      Respi_tot = weight * p_eff_all - npp
+      ! npp = (weight * p_eff_all - Respi_m / 10000.) / (1.+par_c)
+      npp = (gpp_sp - Respi_m / 10000.) / (1.+par_c)
+	  ! Respi_tot = weight * p_eff_all - npp
+	  Respi_tot = gpp_sp - npp
       V_scrown =  A * (par_betas*Lc)
 ! note that this equation has changed AM 15.5.2015
       V_bole = (A+B+sqrt(A*B)) * Hc /2.9
@@ -652,7 +653,7 @@ subroutine prebas_v0(nYears,nLayers,nSp,siteInfo,pCrobas,initVar,thinning,output
   STAND(35) = B
   STAND(36) = Light
   STAND(42) = Vold* min(1.,-dN*step/Nold)
-  STAND(44) = p_eff
+  STAND(44) = gpp_sp
 else
   STAND(8:21) = 0. !#!#
   STAND(23:37) = 0. !#!#
@@ -1085,7 +1086,8 @@ subroutine prebas_v0(nYears,nLayers,nSp,siteInfo,pCrobas,initVar,thinning,output
 
  ! write(2,*) "here4"
 
-modOut(:,46,:,1) = modOut(:,44,:,1) - modOut(:,9,:,1) - modOut(:,45,:,1)
+modOut(:,46,:,1) = modOut(:,44,:,1) - modOut(:,9,:,1) - modOut(:,45,:,1) !!Gpp is not smoothed
+!modOut(:,46,:,1) = modOut(:,18,:,1) - modOut(:,45,:,1) !!!everything smoothed
 
 ! write(2,*) "here5"
 
@@ -1100,5 +1102,3 @@ subroutine prebas_v0(nYears,nLayers,nSp,siteInfo,pCrobas,initVar,thinning,output
 end subroutine
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-

src/C_regionPrebas.f90

@@ -2,7 +2,7 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !subroutine bridging
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-subroutine regionPrebas(siteOrder,HarvLim,minDharv,multiOut,nSites,nClimID,nLayers,maxYears,maxThin, &
+subroutine regionPrebas(siteOrder,HarvLim,minDharv,multiOut,nSites,areas,nClimID,nLayers,maxYears,maxThin, &
 		nYears,thinning,pCrobas,allSP,siteInfo, maxNlayers, &
 		nThinning,fAPAR,initClearcut,fixBAinitClarcut,initCLcutRatio,ETSy,P0y, initVar,&
 		weatherPRELES,DOY,pPRELES,etmodel, soilCinOut,pYasso,&
@@ -18,7 +18,7 @@ subroutine regionPrebas(siteOrder,HarvLim,minDharv,multiOut,nSites,nClimID,nLaye
 real (kind=8), intent(in) :: weatherPRELES(nClimID,maxYears,365,5),HarvLim(maxYears),minDharv
  integer, intent(in) :: DOY(365),etmodel
  real (kind=8), intent(in) :: pPRELES(30),pCrobas(npar,allSP)
- real (kind=8), intent(inout) :: siteInfo(nSites,7)
+ real (kind=8), intent(inout) :: siteInfo(nSites,7), areas(nSites)
  real (kind=8), intent(in) :: thinning(nSites,maxThin,8),pAWEN(12,allSP)
  real (kind=8), intent(inout) :: dailyPRELES(nSites,(maxYears*365),3)
  real (kind=8), intent(inout) :: initClearcut(nSites,5),fixBAinitClarcut(nSites),initCLcutRatio(nSites,maxNlayers)	!initial stand conditions after clear cut. (H,D,totBA,Hc,Ainit)
@@ -28,7 +28,7 @@ subroutine regionPrebas(siteOrder,HarvLim,minDharv,multiOut,nSites,nClimID,nLaye
 ! integer, intent(in) :: siteThinning(nSites)
  integer, intent(inout) :: nThinning(nSites)
  real (kind=8), intent(out) :: fAPAR(nSites,maxYears)
- real (kind=8), intent(inout) :: initVar(nSites,6,maxNlayers),P0y(nClimID,maxYears),ETSy(nClimID,maxYears)!,par_common
+ real (kind=8), intent(inout) :: initVar(nSites,6,maxNlayers),P0y(nClimID,maxYears,2),ETSy(nClimID,maxYears)!,par_common
  real (kind=8), intent(inout) :: multiOut(nSites,maxYears,nVar,maxNlayers,2)
  real (kind=8), intent(inout) :: soilCinOut(nSites,maxYears,5,3,maxNlayers),soilCtotInOut(nSites,maxYears) !dimensions = nyears,AWENH,treeOrgans(woody,fineWoody,Foliage),species
  real (kind=8) :: soilC(nSites,maxYears,5,3,maxNlayers),soilCtot(nSites,maxYears) !dimensions = nyears,AWENH,treeOrgans(woody,fineWoody,Foliage),species
@@ -106,15 +106,15 @@ subroutine regionPrebas(siteOrder,HarvLim,minDharv,multiOut,nSites,nClimID,nLaye
 	if(prebasVersion(i)==0.) then
 	  call prebas_v0(1,nLayers(i),allSP,siteInfo(i,:),pCrobas,initVar(i,:,1:nLayers(i)),&
 		thinningX(1:az,:),output(1,:,1:nLayers(i),:),az,maxYearSite,fAPAR(i,ij),initClearcut(i,:),&
-		fixBAinitClarcut(i),initCLcutRatio(i,1:nLayers(i)),ETSy(climID,ij),P0y(climID,ij),&
+		fixBAinitClarcut(i),initCLcutRatio(i,1:nLayers(i)),ETSy(climID,ij),P0y(climID,ij,:),&
 		weatherPRELES(climID,ij,:,:),DOY,pPRELES,etmodel, &
 		soilC(i,ij,:,:,1:nLayers(i)),pYasso,pAWEN,weatherYasso(climID,ij,:),&
 		litterSize,soilCtot(i,ij),&
 		defaultThinX,ClCutX,inDclct(i,:),inAclct(i,:),dailyPRELES(i,(((ij-1)*365)+1):(ij*365),:),yassoRun(i))
 	elseif(prebasVersion(i)==1.) then
 	  call prebas_v1(1,nLayers(i),allSP,siteInfo(i,:),pCrobas,initVar(i,:,1:nLayers(i)),&
 		thinningX(1:az,:),output(1,:,1:nLayers(i),:),az,maxYearSite,fAPAR(i,ij),initClearcut(i,:),&
-		fixBAinitClarcut(i),initCLcutRatio(i,1:nLayers(i)),ETSy(climID,ij),P0y(climID,ij),&
+		fixBAinitClarcut(i),initCLcutRatio(i,1:nLayers(i)),ETSy(climID,ij),P0y(climID,ij,:),&
 		weatherPRELES(climID,ij,:,:),DOY,pPRELES,etmodel, &
 		soilC(i,ij,:,:,1:nLayers(i)),pYasso,pAWEN,weatherYasso(climID,ij,:),&
 		litterSize,soilCtot(i,ij),&
@@ -165,7 +165,7 @@ subroutine regionPrebas(siteOrder,HarvLim,minDharv,multiOut,nSites,nClimID,nLaye
 	initVar(i,1,1:nLayers(i)) = output(1,4,1:nLayers(i),1)
 	initVar(i,2,1:nLayers(i)) = output(1,7,1:nLayers(i),1)
 	initVar(i,3:6,1:nLayers(i)) = output(1,11:14,1:nLayers(i),1)
-	HarvArea = HarvArea + sum(output(1,37,1:nLayers(i),1))
+	HarvArea = HarvArea + sum(output(1,37,1:nLayers(i),1))* areas(i)
  end do !iz i
 
 
@@ -186,7 +186,7 @@ subroutine regionPrebas(siteOrder,HarvLim,minDharv,multiOut,nSites,nClimID,nLaye
 if(maxState(siteX)>minDharv .and. ClCut(siteX) > 0.) then
   ! close(10)
 !!   !!clearcut!!
-   HarvArea = HarvArea + sum(multiOut(siteX,ij,30,1:nLayers(siteX),1))
+   HarvArea = HarvArea + sum(multiOut(siteX,ij,30,1:nLayers(siteX),1))*areas(i)
    multiOut(siteX,ij,37,:,1) = multiOut(siteX,ij,37,1:nLayers(siteX),1) + multiOut(siteX,ij,30,1:nLayers(siteX),1)
    do ijj = 1, nLayers(siteX)
     multiOut(siteX,ij,6:nVar,ijj,2) = multiOut(siteX,ij,6:nVar,ijj,1)