Skip to content

Commit be39875

Browse files
authored
Merge pull request #1 from ForModLabUHel/master
Before Euca_PARA
2 parents f71496d + 0790c55 commit be39875

11 files changed

Lines changed: 1092 additions & 148 deletions

R/clcut.r

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
ClCutD_Pine <- function(ETSmean,ETSthres,siteType){
1+
ClCutD_Pine <- function(ETSmean,ETSthres,siteType){ #TEST DELETE THIS COMMENT LAuri
22
if(siteType<=3 & ETSmean>=ETSthres) inDclct <- ClCut_pine[1,1]
33
if(siteType==4 & ETSmean>=ETSthres) inDclct <- ClCut_pine[2,1]
44
if(siteType>=5 & ETSmean>=ETSthres) inDclct <- ClCut_pine[3,1]

R/multiSitePrebas.r

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -62,7 +62,7 @@ InitMultiSite <- function(nYearsMS,
6262
nVar <- length(varNam)
6363

6464
nClimID <- length(unique(climIDs))
65-
if(!all((1:nClimID) %in% climIDs) | length(climIDs) != nSites) return("check consistency between weather inputs and climIDs")
65+
if(!all((1:nClimID) %in% climIDs) | length(climIDs) != nSites) warning("check consistency between weather inputs and climIDs")
6666
if(nClimID == 1){
6767
nClimID = 2
6868
climIDs[1] <- 2
@@ -148,12 +148,12 @@ InitMultiSite <- function(nYearsMS,
148148
c(ClCutD_Pine(ETSmean[climIDs[i]],ETSthres,siteInfo[i,3]),
149149
ClCutD_Spruce(ETSmean[climIDs[i]],ETSthres,siteInfo[i,3]),
150150
ClCutD_Birch(ETSmean[climIDs[i]],ETSthres,siteInfo[i,3]),
151-
NA,NA,NA,NA) ###"fasy","pipi","eugl","rops")
151+
NA,NA,NA,NA,NA) ###"fasy","pipi","eugl","rops","popu")
152152
if(ClCut[i]==1 & all(is.na(inAclct[i,]))) inAclct[i,] <-
153153
c(ClCutA_Pine(ETSmean[climIDs[i]],ETSthres,siteInfo[i,3]),
154154
ClCutA_Spruce(ETSmean[climIDs[i]],ETSthres,siteInfo[i,3]),
155155
ClCutA_Birch(ETSmean[climIDs[i]],ETSthres,siteInfo[i,3]),
156-
80,50,13,30) ###"fasy","pipi","eugl","rops" )
156+
80,50,13,30,50) ###"fasy","pipi","eugl","rops","popu" )
157157
if(any(!is.na(inDclct[i,]))) inDclct[i,is.na(inDclct[i,])] <- max(inDclct[i,],na.rm=T)
158158
if(all(is.na(inDclct[i,]))) inDclct[i,] <- 9999999.99
159159
if(any(!is.na(inAclct[i,]))) inAclct[i,is.na(inAclct[i,])] <- max(inAclct[i,],na.rm=T)
@@ -231,7 +231,8 @@ InitMultiSite <- function(nYearsMS,
231231
multiInitVar <- aaply(multiInitVar,1,findHcNAs,pHcMod)
232232
}
233233

234-
234+
###age cannot be lower than 1 year
235+
multiInitVar[,2,][which(multiInitVar[,2,]<1)] <- 1
235236

236237
####compute A
237238
for(ikj in 1:maxNlayers){

R/parameters.r

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -116,7 +116,7 @@ pHcM[1:3,7] <- c(0.04237,-0.13308,0.31382)
116116
pHcM[1:3,8] <- c(0.04237,-0.13308,0.31382)
117117

118118
litterSizeDef <- matrix(0.,3,8,dimnames = list(NULL,speciesNam))
119-
litterSizeDef[1,] <- c(30,30,10,30,30,20,20,20)
119+
litterSizeDef[1,] <- c(10,10,5,10,10,7,7,7)
120120
litterSizeDef[2,] <- 2
121121

122122
ClCut_birch <- matrix(NA,2,4)

R/prebas.r

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,3 @@
1-
21
prebas <- function(nYears,
32
pCROBAS = pCROB,
43
pHcMod = pHcM,

R/utilStaff.r

Lines changed: 192 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -92,3 +92,195 @@ varNames <- c('siteID','gammaC','sitetype','species','ETS' ,'P0','age', 'DeadWo
9292
aP = colSums(matrix(Precip,365,nYears))
9393
return(aP)
9494
}
95+
96+
97+
98+
### Function for the calculation of monthly mean values from a daily model output variable GPP
99+
monthlyGPP <- function(GPPx,func="sum"){
100+
monthsDays <- c(rep(1,31),rep(2,28),rep(3,31),rep(4,30),rep(5,31),rep(6,30), # Number of days in each month
101+
rep(7,31),rep(8,31),rep(9,30),rep(10,31),rep(11,30),rep(12,31))
102+
GPPbyYear <- matrix(GPPx,365,length(GPPx)/365) # Divide GPP data into a matrix so that rows are days and columns are years
103+
GPPm = unlist(apply(GPPbyYear, 2, function(x) # Applies aggregation on the columns of input matrix of GPPbyYear.
104+
aggregate(x,by=list(monthsDays),FUN=func)$x)) # Aggregation is done by groups in the variable monthsDays, as there is no timestamp on datapoints of GPP
105+
GPPm <- as.vector(GPPm)
106+
return(GPPm)
107+
}
108+
109+
monthlyWeather <- function(varx,func){ # This function works similarly as the monthlyGPP above but with different data set
110+
monthsDays <- c(rep(1,31),rep(2,28),rep(3,31),rep(4,30),rep(5,31),rep(6,30),
111+
rep(7,31),rep(8,31),rep(9,30),rep(10,31),rep(11,30),rep(12,31))
112+
# varByYear <- matrix(varx,365,length(varx)/365)
113+
varM = unlist(apply(varx, 1, function(x)
114+
aggregate(x,by=list(monthsDays),FUN=func)$x))
115+
varM <- as.vector(varM)
116+
return(varM)
117+
}
118+
119+
### A post-process function for modifying prebas model output for meeting the Digital Twin Earth project requirements
120+
monthlyFluxes <- function(modOut,weatherOption=1){
121+
cueGV <- 0.5 ###carbon use efficiency (NPP/GPP) of ground vegetation
122+
123+
# if(is.matrix(mGPP)){
124+
##aggregating total GPP (tree layers + GV) at monthly time step
125+
mGPPtot <- t(apply(modOut$dailyPRELES[,,1],1,monthlyGPP,"sum")) # Daily GPP is converted into daily GPP using the function monthlyGPP, which is defined above
126+
# $dailyPRELES is a output matrix from the prebas model
127+
### Calculate CUE for different layers
128+
cueTrees <- modOut$multiOut[,,18,,1]/modOut$multiOut[,,44,,1] # CUE = calculate npp / GPPspecies, but include data only from the standing trees
129+
cueTrees[which(is.na(cueTrees))] <- 0. # Replace NAs with a value
130+
cueAll <- pGPP <- array(NA,dim=c(modOut$nSites,modOut$maxYears,(modOut$maxNlayers+1))) # Create a pair of 3D arrays: x=1:7, y=1:150, z=1:4
131+
cueAll[,,1] <- cueGV # When layer is 1, then CUE is cueGV?
132+
cueAll[,,2:(modOut$maxNlayers+1)] <- cueTrees # ... and otherwise its cueTrees?
133+
134+
### Calculate total GPP (totGPP) and ratio of GPP at a layer/total GPP (pGPP)
135+
totGPP <- apply(modOut$multiOut[,,44,,1],1:2,sum) + modOut$GVout[,,3] # Sum of GPP for each year and site in a 2d array
136+
pGPP[,,1] <- modOut$GVout[,,3]/totGPP # Ratio between CUE of vegetation and CUE total?
137+
138+
for(i in 1:modOut$maxNlayers){ # Ratio between CUE of vegetation and pine or spruce or mixed
139+
pGPP[,,(i+1)] <- modOut$multiOut[,,44,i,1] / totGPP
140+
}
141+
### Create arrays for monthly GPP and NPP values
142+
mGPP <- mNPP <- array(NA,dim=c(modOut$nSites,modOut$maxYears*12,(modOut$maxNlayers+1))) # Create another 3D array, similar as before.
143+
for(i in 1:dim(pGPP)[3]){ # Loop goes through forest layers 1-4 with i
144+
for(ij in 1:modOut$nSites){ # ... and for every layers it goes through different sites ij
145+
mGPP[ij,,i] <- rep(pGPP[ij,,i],each=12) * mGPPtot[ij,] # Calculate monthly GPP for different layers by multiplying layer's share of GPP with monthly total GPP
146+
mNPP[ij,,i] <- mGPP[ij,,i] * rep(cueAll[ij,,i],each=12) # Calculate monthly NPP for different layers by multiplying monthly total GPP with layers yearly CUE
147+
}
148+
}
149+
150+
151+
###litterfall calculations lit is splitted to August and September
152+
# Woody litter
153+
# These 5 variables below are different litter categories, have 3d structure: site,litter variable,layer
154+
Lf <- aperm(apply(modOut$multiOut[,,26,,1],c(1,3),mLit,months=8:9),c(2,1,3)) # Runs model output variable Litter_fol through the 'mLit' function and transforms the containing array by switching places between rows and columns
155+
Lfr <- aperm(apply(modOut$multiOut[,,27,,1],c(1,3),mLit,months=8:9),c(2,1,3)) # The mLit function moves litter fall to the August and September time period in the model output
156+
Lnw <- Lf + Lfr # Non-woody litter = The sum of Foliage litter (var. Litter_fol) and Fine root litter (var. Litter_fr) prebas output variables
157+
Lfw <- aperm(apply(modOut$multiOut[,,28,,1],c(1,3),mLit,months=8:9),c(2,1,3)) # fine woody litter
158+
Lw <- aperm(apply(modOut$multiOut[,,29,,1],c(1,3),mLit,months=8:9),c(2,1,3)) # Coarse woody litter
159+
160+
### Create a input array for the Yasso model with the litter fall information
161+
nYears <- dim(modOut$multiOut)[2]
162+
nMonths <- nYears * 12
163+
nSites <- modOut$nSites
164+
nLayers <- dim(modOut$multiOut)[4]
165+
litter <- array(0.,dim=c(nSites,nMonths,nLayers,3))
166+
litter[,,,1] <- Lnw # Non-woody litter
167+
litter[,,,2] <- Lfw # Branch litter
168+
litter[,,,3] <- Lw # Woody litter
169+
# litter <- litter*12
170+
### Prepare also other initialization information for Yasso
171+
species <- modOut$multiOut[,1,4,,1]
172+
nSp <- max(species)
173+
litterSize <- modOut$litterSize[,1:nSp]
174+
soilC <- array(0.,dim=c(nSites,(nMonths+1),5,3,nLayers))
175+
soilC[,1,,,] <- modOut$soilC[,1,,,] ###initialize with soilC from simulations
176+
nClimID <- modOut$nClimID
177+
climIDs <- modOut$siteInfo[,2]
178+
weatherYasso <- array(0,dim=c(nClimID,nMonths,3))
179+
if(weatherOption==1){
180+
weatherYasso[,,1] <- t(apply(modOut$weather[,,,2],1,monthlyWeather,"mean")) ###monthly tempearture
181+
weatherYasso[,,2] <- t(apply(modOut$weather[,,,4],1,monthlyWeather,"sum")) *12 ###monthly precipitation ###*12 rescales to annual mm
182+
}else if(weatherOption==2){
183+
weatherYasso[,,1] <- t(apply(modOut$weather[,,,2],1,monthlyWeather,"mean")) ###monthly tempearture
184+
weatherYasso[,,2] <- t(apply(modOut$weather[,,,4],1,monthlyWeather,"sum")) *12 ###monthly precipitation ###*12 rescales to annual mm
185+
weatherYasso[,,3] <- (t(apply(modOut$weather[,,,2],1,monthlyWeather,"max")) -
186+
t(apply(modOut$weather[,,,2],1,monthlyWeather,"min"))) /2
187+
}else if(weatherOption==3){
188+
for(i in 1:nSites){
189+
for(j in 1:nYears){
190+
weatherYasso[i,(((j-1)*12)+1):(j*12),1] <- modOut$weatherYasso[i,j,1]
191+
weatherYasso[i,(((j-1)*12)+1):(j*12),2] <- modOut$weatherYasso[i,j,2]
192+
weatherYasso[i,(((j-1)*12)+1):(j*12),3] <- modOut$weatherYasso[i,j,3]
193+
}
194+
}
195+
}
196+
### Run the Yasso model, which is a function in the src/A_routines.90 file
197+
soilCtrees <- .Fortran("runYassoMonthly",litter=as.array(litter*12), ###*12 rescale monthly litter to annual units
198+
litterSize=as.array(litterSize),
199+
nMonths=as.integer(nMonths),
200+
nLayers=as.integer(nLayers),
201+
nSites=as.integer(nSites),
202+
nSp=as.integer(nSp),
203+
species=as.matrix(species),
204+
nClimID=as.integer(nClimID),
205+
climIDs=as.integer(climIDs),
206+
pAWEN=as.matrix(parsAWEN),
207+
pYasso=as.double(pYAS),
208+
climate=as.matrix(weatherYasso),
209+
soilC=as.array(soilC))
210+
211+
###calculate soil C for gv
212+
fAPAR <- modOut$fAPAR # fAPAR from the preles model within prebas. 3d array, value for each year, but no layer dimension
213+
fAPAR[which(is.na(modOut$fAPAR),arr.ind = T)] <- 0. # Convert NAs to values
214+
AWENgv <- array(NA,dim=c(dim(modOut$fAPAR),4)) # Add the layer dimension on the fAPAR array
215+
p0 = modOut$multiOut[,,6,1,1] # Annual potential photosynthesis, 3d array: site,year,value
216+
ETSy = modOut$multiOut[,,5,1,1]
217+
218+
for(ij in 1:nYears){
219+
AWENgv[,ij,] <- t(sapply(1:nrow(fAPAR), function(i) .Fortran("fAPARgv",fAPAR[i,ij],
220+
ETSy[i,ij],modOut$siteInfo[i,2],
221+
0,0,p0[i,ij],rep(0,4))[[7]]))
222+
}
223+
224+
mAWEN <- array(0,dim=c(nSites,nMonths,5))
225+
mAWEN[,,1:4] <- aperm(apply(AWENgv,c(1,3),mLit,months=6:9),c(2,1,3))*12 ###*12 rescale monthly litter to annual units
226+
mGVit <- t(apply(modOut$GVout[,,2],1,mLit,months=6:9))
227+
###calculate steady state soil C per GV
228+
# ststGV <- matrix(NA,nSites,5)
229+
soilGV <- array(0,dim=c(nSites,(nMonths+1),5))
230+
231+
232+
soilCgv <- .Fortran("runYassoAWENinMonthly",
233+
mAWEN=as.array(mAWEN),
234+
nMonths=as.integer(nMonths),
235+
nSites=as.integer(nSites),
236+
litSize=0.,
237+
nClimID=as.integer(nClimID),
238+
climIDs=as.integer(climIDs),
239+
pYasso=as.double(pYAS),
240+
climate=as.matrix(weatherYasso),
241+
soilCgv=as.array(soilGV))
242+
####add gvsoilc to first layer foliage soilC
243+
# check in normal runs where ground vegetation soilC is calculated
244+
soilCtot <- apply(soilCtrees$soilC,1:2,sum) + apply(soilCgv$soilCgv,1:2,sum)
245+
mRhTot <- soilCtot[,1:nMonths]/10 - soilCtot[,2:(nMonths+1)]/10 +
246+
apply(litter,1:2,sum)/10 + mGVit/10
247+
mNPPtot <- apply(mNPP,1:2,sum)
248+
mRaTot <- mGPPtot - mNPPtot
249+
mNEPtot <- mNPPtot - mRhTot
250+
return(list(mGPP=mGPPtot,mNPP=mNPPtot,mRa=mRaTot,mRh=mRhTot,
251+
mNEP=mNEPtot,soilC=soilCtot))
252+
}
253+
254+
255+
256+
### This function is used for adding yearling litter fall to the fall period on a timeseries data set, used in the monthlyFluxes function
257+
mLit <- function(aLit,months=8:9){
258+
nYears <- length(aLit) # Gets the length of the array
259+
nMonths <- length(months)
260+
lit <- matrix(0,12,nYears) # Creates a matrix by 12 rows, nYears columns and fills it with zeros
261+
lit[months,] <- matrix(aLit/nMonths,nMonths,nYears,byrow = T) # This appears to create a new matrix within a matrix, on rows 8-9. Is this just used to fill these rows with data?
262+
lit <- as.vector(lit) # Matrix is converted intoa vector, which removes the table structure and appears to create a timeseries data set
263+
return(lit)
264+
}
265+
266+
267+
###Function to calculate basal weighted mean of a PREBAS output
268+
### modOut = multisite PREBAS output
269+
### varX = index of the variable for which the basal area weighted mean needs to be calculated
270+
baWmean <- function(modOut,varX){
271+
###calculates basal area weighted mean for single site runs
272+
if(class(modOut)=="prebas"){
273+
weightXs <- as.matrix(apply(modOut$output[,13,,1],1,FUN=function(vec)vec/sum(vec)))
274+
weightXs <- t(weightXs)
275+
weightXs[which(is.na(weightXs))] <- 0.
276+
weigthedMean <- rowSums(modOut$output[,varX,,1]*weightXs)
277+
}else{
278+
###calculates basal area weighted mean for multi site runs
279+
weightXs <- apply(modOut$multiOut[,,13,,1],1:2,FUN=function(vec)vec/sum(vec))
280+
weightXs <- aperm(weightXs,c(2:3,1))
281+
weightXs[which(is.na(weightXs))] <- 0.
282+
weigthedMean <- apply(modOut$multiOut[,,varX,,1]*weightXs,1:2,sum)
283+
}
284+
return(weigthedMean)
285+
}
286+

0 commit comments

Comments
 (0)