Program Biomassa
      Implicit None
      
      !Declarando Vari veis do modelo
      
      Real RUE, K,IC,U,IAF,MS,PP,PPTotal, PATotal, a, r, t, aPAR(365)
      Real PAR(365), tPAR(365), rPAR(365), Qg(365), Qo(365),NN
      Real Alfa, DensAgua, CalorEvapo, Patm, TetensA, TetensB, TetensC
      Real ConstPsicr, Delta, FatorRed
      Real UmCC, UmPMP, UmCrit, ProfRad, Armaz, Umidade
      Real Chuva(365), Rn(365), Temp(365), Es
      Real ETpt, ETreal, PPObs(11), PAObs(11)

      
      Integer DJ, NI, DJSow, Ciclo, i
      
      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
            UmCC = 0.28 !m3/m3
      UmPMP = 0.10
      UmCrit = 0.19
      ProfRad = 0.5 !m
      
      !Constantes 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
      
      !Produtividades observadas
       PPObs(1)=10030
       PPObs(2)=8861
       PPObs(3)=9012
       PPObs(4)=9718
       PPObs(5)=8605
       PPObs(6)=8305
       PPObs(7)=9525
       PPObs(8)=10177
       PPObs(9)=7739
       PPObs(10)=10129
       PPObs(11)=11311

       PAObs(1)=8999
       PAObs(2)=7884
       PAObs(3)=8229
       PAObs(4)=8058
       PAObs(5)=6161
       PAObs(6)=6498
       PAObs(7)=7539
       PAObs(8)=9898
       PAObs(9)=7461
       PAObs(10)=8730
       PAObs(11)=8851
       
      !Arquivos com dados meteorol¢gicos

      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"
      

      !Defini‡Æo do ciclo da cultura
      DJSow=60
      Ciclo=120


      !Arquivo com dados de sa¡da
      Open (unit=2, File="Saida_Modelo_Biomassa.txt", status="replace")
      
      
      !Leitura (importa‡Æo) dos dados meteorol¢gicos
      
      do i=1,11

        Umidade = 0.25 !umidade inicial em m3/m3
        Armaz = Umidade * ProfRad * 1000  !mm
        

        Open (unit=1, File=File1(i), Status="old")
        read (1, "(A200)") Lixo !lendo o cabecalho do arquivo Diarios2013.prn e jogando na vari vel Lixo
                !Lendo os dados do arquivos File1
        Do DJ=1,365
         read (1,*) NI, NI, NI, Temp(DJ), NN, NN, NN, NN, NN, NN, NN, 
     C    Chuva(DJ), Qg(DJ),Rn(DJ)
        Enddo

        PPtotal = 0.
        PAtotal = 0.

      
      
        Do DJ=1,365

          IF (DJ > DJSow .and. DJ<(DJSow+Ciclo)) Then
            !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
            !write (*,*) Chuva(DJ), ETreal
            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
             PpTotal = PPTotal + Pp
             PaTotal = PaTotal + Pp*FatorRed

        End if
             
             
      Enddo

        write (*,"(2F12.1,F12.3)") PpTotal, PaTotal, PaTotal/PpTotal
        write (2,"(2F12.1),F12.3") PpTotal, PaTotal, PaTotal/PpTotal


        Close(1)

      enddo
      
      read*

      Close(2)
      End Program Biomassa