diff --git a/R/multiSitePrebas.r b/R/multiSitePrebas.r index 5f4f1e73..ae5722e7 100644 --- a/R/multiSitePrebas.r +++ b/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,7 +143,7 @@ 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)], @@ -148,8 +151,16 @@ InitMultiSite <- function(nYearsMS, 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) } diff --git a/R/prebas.r b/R/prebas.r index a04be980..daced420 100644 --- a/R/prebas.r +++ b/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), diff --git a/src/B_prebas_v0.f90 b/src/B_prebas_v0.f90 index 43671145..3a08b7a5 100644 --- a/src/B_prebas_v0.f90 +++ b/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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - diff --git a/src/C_regionPrebas.f90 b/src/C_regionPrebas.f90 index 862b611d..196ab71c 100644 --- a/src/C_regionPrebas.f90 +++ b/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,7 +106,7 @@ 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),& @@ -114,7 +114,7 @@ subroutine regionPrebas(siteOrder,HarvLim,minDharv,multiOut,nSites,nClimID,nLaye 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)