Estou tendo dificuldades em acertar o formato que devo organizar os outputs de um looping que fiz. 

Tenho um banco de dados com datas e algumas variáveis climáticas. Dessas datas, preciso simular alguns cálculos para dias específicos. São 36 simulações por ano, durante um período de 34 anos, totalizando 1224 simulações.

Para alcançar esse objetivo estou trabalhando com 3 comandos for, um dentro do outro.

O primeiro seleciona as datas/posições que devem ser feitos os cálculos.

O segundo executa os devidos cálculos para as datas selecionadas e o ultimo integra os valores calculados, gerando a informação de interesse.

Acredito que a logica de cálculos dentro de cada for() está ok, mas a saída não está puxando os valores que o looping calcula.

Segue o script no corpo do email e, em anexo, um exemplo de um ano do banco de dados que está sendo utilizado.

Desde já, agradeço a colaboração dos colegas.

Abs

#--------------------Remover objetos salvos na memória do R--------------------#

rm(list = ls())

#--------------------Seleção de diretorios e arquivos de interesse--------------------#

a = read.table("x.txt", sep = "\t", header = T)

#--------------------Calculo das Variaves Astronomicas--------------------#

head(a)

a$DATE = as.Date(levels(a$DATE))[a$DATE]

ano = format(a$DATE, trim = T, '%Y'); a$ano=ano

mes = format(a$DATE, trim = T, '%m') ; a$mes=mes

dia = format(a$DATE, trim = T, '%d') ; a$dia=dia 

dece = ifelse(dia<11, 1,ifelse(dia<21,2,3)); a$decendio=dece

lat = -7.53

corr = pi/180

decl = 23.45*sin(corr*((a$DiaJuliano-80)*360/365)); decl

hn = 1/corr*acos(-tan(corr*lat)*tan(corr*decl))

a$N = 2*hn/15

s = (1+0.033*cos(a$DiaJuliano*360/365))

t = corr*hn*sin(corr*lat)*sin(corr*decl)

u = cos(corr*lat)*cos(corr*decl)*sin(corr*hn)

a$Qo = 37.6*s*(t+u)

head(a)

x = 0.25

b = 0.50

n.est = (a$SRAD/a$Qo-x)*a$N/b; n.est

a$n<- ifelse((n.est)<0,1, n.est); head(a)

a$n.N = a$n/a$N; head(a)

#--------------------Transformar MJ em caloria--------------------#

a$Srad.cal = round((a$SRAD*1000000)/41800,1); head(a) #round é para delimitar o numero de casas apos a virgula

a$sem = paste0(dia, '-', mes)

a$TMED = (((a$TMAX+a$TMIN))/2)

head(a)

#-------------------- Determinacao das datas de simulacao--------------------#

ndc = 140

datas = 1:366 # referentes aos dia juliano

dim(a)

anos = table(a$ano); anos

potencial = matrix(NA, nrow = 366, ncol = 34 )

rownames(potencial) = 1:366

colnames(potencial) = c(1980:2013)

head(potencial); dim(potencial)

head(a)

for (i in (a$DiaJuliano)){

  

  epocas = a[which(  a$sem == '01-01'|a$sem == '11-01' |a$sem =='21-01'|

                     a$sem == '01-02'|a$sem == '11-02' |a$sem =='21-02'|

                     a$sem == '01-03'|a$sem == '11-03' |a$sem =='21-03'|

                     a$sem == '01-04'|a$sem == '11-04' |a$sem =='21-04'|

                     a$sem == '01-05'|a$sem == '11-05' |a$sem =='21-05'|

                     a$sem == '01-06'|a$sem == '11-06' |a$sem =='21-06'|

                     a$sem == '01-07'|a$sem == '11-07' |a$sem =='21-07'|

                     a$sem == '01-08'|a$sem == '11-08' |a$sem =='21-08'|

                     a$sem == '01-09'|a$sem == '11-09' |a$sem =='21-09'|

                     a$sem == '01-10'|a$sem == '11-10' |a$sem =='21-10'|

                     a$sem == '01-11'|a$sem == '11-11' |a$sem =='21-11'|

                     a$sem == '01-12'|a$sem == '11-12' |a$sem =='21-12'& a$ano>=1980),]

  

  semear = length(epocas[,1])

  

  for (j in length(semear)){

    

    ini = as.numeric(rownames(epocas[j]))

    fim = (ini+ndc)-1

    

    pac = a[(ini:fim),]

    

    #--------------------Inserindo as Produtividades Parciais--------------------#

    #----------Inserindo influencia da nebulosidade----------#

    pac$ctc[j] = -4.16+(0.4325*pac$TMED[j])-(0.00725*((a$TMED[j])^2))

    pac$ctn[j] = -1.064+(0.173*a$TMED[j])-(0.0029*((a$TMED[j])^2))

    pac$PPBc[j] = (107.2+(0.36*pac$Qo[j]*1000000)/41800)*pac$ctc[j]*pac$n.N[j]

    pac$PPBn[j] = (31.7+(0.219*pac$Qo[j]*1000000)/41800)*pac$ctn[j]*(1-(pac$n.N[j]))

    pac$PPB[j] = round((pac$PPBc[j])+(pac$PPBn[j]),1); head(pac)

    

    #--------------------Calculo da Produtividade - base seca --------------------#

    pac$c.resp[j] = ifelse(a$TMED[j] < 20, 0.6, 0.5)

    pac$c.colheita[j] = 0.5

    umidade = 13

    c.um = (1/(1-0.01*umidade))

    pac$iaf[j] = 5

    pac$c.iaf[j] = (0.0093+0.185*pac$iaf[j])-(0.0175*pac$iaf[j]^2)

    pac$PP.Liquida_seca = NA

    pac$PP.Liquida_seca[1] = round(pac$PPB[1]*pac$c.resp[1]*pac$c.colheita[1]*pac$c.iaf[1],1); head(pac) # em kg/ha

    

    for (w in 2:ndc){

      

      pac$PP.Liquida_seca[w] = round(pac$PPB[w]*pac$c.resp[w]*pac$c.colheita[w]*pac$c.iaf[w]+

                                       pac$PP.Liquida_seca[w-1],1)

    }

    #--------------------Soma da PP--------------------#

    pac$PP_final = round(pac$PP.Liquida_seca[w]*c.um,1)

    dim(pac)

    PP = pac[ndc,32]

    potencial[i,j] = PP

    }

  print(i)

}

View(pac)

View(potencial)