Program Biomassa
      Implicit None
      
      !Declarando Vari veis do modelo
      Real RUE, K,IC,U,IAF,MS,PP,PPTotal, PATotal, a, r, t, aPAR(4015)
      Real PAR(4015), tPAR(4015), rPAR(4015), Qg(4015), Qo(4015),NN
      Real Alfa, DensAgua, CalorEvapo, Patm, TetensA, TetensB, TetensC
      Real ConstPsicr, Delta, FatorRed
      Real UmCC, UmPMP, UmCrit, ProfRad, Armaz, Umidade
      Real Chuva(4015), Rn(4015), Temp(4015), Es
      Real ETpt, ETreal,PARtotal
      Integer DJ, NI, DJSow, Ciclo, i, nd, dia, j
      Character*200 File1(20), Lixo
      
      !Constantes para c lculo do ET-Priestley Taylor
      Alfa = 1.26
      DensAgua = 1000. !kg/m3
      CalorEvapo = 2.45e6 !J/kg
      Patm = 0.94e5  !Pa
      TetensA = 610.8  !Pa
      TetensB = 17.27
      TetensC = 237.3
      ConstPsicr = 6.66e-4 * Patm  !Pa/K
      
      !Constantes para balan‡o h¡drico
      Umidade = 0.25 !umidade inicial em m3/m3
      UmCC = 0.28 !m3/m3
      UmPMP = 0.10
      UmCrit = 0.19
      ProfRad = 0.5 !m

      Ciclo = 120 !dias
      DJSow = 320 !Definindo data de semeadura para 1/3
      !Print *, 'Informe a Data de Semeadura:'
      ! Read *, DJSow
      ! Print *, 'Informe o ciclo da cultura:'
      ! Read *, Ciclo

      !Criando as unidades para leitura de dados de entrada e grava‡Æo
      !dos dados de sa¡da

      File1(1) = "diarios2001.prn"
      File1(2) = "diarios2002.prn"
      File1(3) = "diarios2003.prn"
      File1(4) = "diarios2005.prn"
      File1(5) = "diarios2006.prn"
      File1(6) = "diarios2007.prn"
      File1(7) = "diarios2009.prn"
      File1(8) = "diarios2010.prn"
      File1(9) = "diarios2011.prn"
      File1(10) = "diarios2013.prn"
      File1(11) = "diarios2014.prn"

      !definindo os valores das vari veis da cultura
      RUE = 3.84 !(g[carboidrato/MJ de aPar] para cultura do milho
      IC = 0.33  !adimensional
      U = 0.13   !adimensional
      K = 0.5    !adimensional
      IAF = 3.5  !adimensional
      
      Open (unit=2, File="Saida_Modelo_Biomassa.txt", status="replace")
      
      nd=1
      Do i=1,11
         Open (unit=1, File=File1(i), Status="old")
         Armaz = Umidade * ProfRad * 1000  !mm
         Read (1, 100) Lixo !lendo o cabecalho do arquivo Diarios2013.prn e jogando na vari vel Lixo
100      Format (A200)
         !Lendo os dados do arquivos File1
           Do DJ=1,365
             read (1,*) NI, NI, NI,Temp(nd),NN,NN, NN, NN, NN, NN, NN,
     C       Chuva(nd), Qg(nd),Rn(nd)
             nd = nd + 1
           Enddo
           Close (1)
      Enddo
101   format (4I5, 6F12.3)
      PPtotal = 0.
      PAtotal = 0.

      j = 1
      i = 0
      dia = 0
      write (2,*) 'data de semeadura:',DJSow
      write (2,*) 'ciclo da cultura :',ciclo
      Write (2,*) ' Arquivo/Ano          Ciclo      DJSem    DJCol
     c  PP         PA          PA/PP     PARTotal'

      Do DJ= 1,4015
      !    i = i + 1
      !    If(i > 365+ciclo) Then
      !        dia  = 0
      !        i = 1 + ciclo
      !        PpTotal = 0.
      !        PaTotal = 0.
      !        PARtotal = 0.
      !        DJSow = DJSow + 365
      !    End if
          If (DJ > (DJSow).and.DJ<(DJSow+Ciclo)) Then
          dia = dia + 1
          !C lculo ET-Priestley-Taylor
            If (Rn(DJ)>0) then
              Es=TetensA*exp(TetensB*Temp(DJ)/(Temp(DJ)+TetensC))
              Delta = Es * TetensB*TetensC / ((TetensC+Temp(DJ))**2)
              ETpt = 1000* alfa*(Delta/(Delta+ConstPsicr))*Rn(DJ)*1e6/
     c                                       DensAgua/CalorEvapo       !mm/d
            Else
              ETpt = 0
            Endif
            If (Umidade>=UmCrit) then
               FatorRed = 1 
            Else
               FatorRed = (Umidade-UmPMP)/(UmCrit-UmPMP)
            Endif
            
            ETreal = FatorRed * ETpt
            Armaz = Armaz + Chuva(DJ) - ETreal  !mm

            Umidade = Armaz / (ProfRad*1000)
            if (Umidade>UmCC) then
                Umidade = UmCC
                Armaz = Umidade * ProfRad * 1000 !mm
            endif
             PAR(DJ) = 0.5 * Qg(DJ)
             aPAR(DJ) = PAR(DJ) * (1-EXP(-K*IAF))
             MS = aPAR(DJ) * RUE
             Pp = MS * IC * (1/(1-U)) * 10
             PARtotal = PAR(DJ) + PARtotal
             PpTotal = PPTotal + Pp
             PaTotal = PaTotal + Pp*FatorRed


        ! write (*,*) file1(j)
              !read*
          If (dia == ciclo-1) Then
               !pause
               write (2,102) file1(j),dia+1,DJSow,DJSow+Ciclo,PpTotal,
     c         PaTotal,PaTotal/PpTotal,PARtotal
               write (*,103) file1(j),dia+1,DJSow,DJSow+Ciclo,PpTotal,
     c         PaTotal,PaTotal/PpTotal,PARtotal
               j = j + 1
               dia  = 0
      !        i = 1 + ciclo
               PpTotal = 0.
               PaTotal = 0.
               PARtotal = 0.
               DJSow = DJSow + 365
          End If
        End if
      Enddo
102   format (A15,3I11, 4F12.2)
103   format (A11,3I6, 4F10.2)
      read*
      Close(2)
      End Program Biomassa