!   
!
!  FUNCTIONS:
!  Crop Yield for several seasons - Entry point of console application.
!

!****************************************************************************
!
!  PROGRAM: CropY_Seasonal
!
!  PURPOSE:  Entry point for the console application.
!
!****************************************************************************

    program Crop_growth

    implicit none

        !Declarando Vari veis do modelo
      Real RUE, K, IC, U, MS, PP, PPTotal, PATotal, a, r, t, aPAR(4015)
      Real PAR(4015), tPAR(4015), rPAR(4015), Qg(4015), Qo(4015),dum, Tinf, Tot1, Tot2, Tsup
      Real Alfa, DensAgua, CalorEvapo, Patm, TetensA, TetensB, TetensC
      Real ConstPsicr, Delta, FatorRed, ftemp
      Real UmCC, UmPMP, UmCrit, ProfRad, Armaz, Umidade, taxaPP
      Real Chuva(4015), Rn(4015), Temp(4015), Es
      Real ETpt, ETreal, PARtotal, IAF(4)
      Integer DJ,  DJSow, Ciclo, i, nd, dia , das, safra
      Character*200 WeatherF(20), dumt
      
      	!Temperaturas Cardinais
	    Tinf = 9.0 ! Temperatura Basal Inferior (Tb)
	    Tot1 = 26.0 !Temperatura Otima 1 (Tot1)
	    Tot2 = 32.0 !Temperatura Otima 2 (Tot2)
	    Tsup = 38.0 ! Temperatura Basal Superior (TB)
      
      !Constantes para calculo da evapotranspiracao pelo metodo de 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 !m3/m3
      UmCrit  = 0.19 !m3/m3
      ProfRad = 0.5  !m
      IAF(1)  = 1.   !m2/m2
      IAF(2)  = 2.5  !m2/m2
      IAF(3)  = 4.5  !m2/m2
      IAF(4)  = 3.5  !m2/m2

      !Ciclo = 120 !dias
      !DJSow = 300 !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

      WeatherF(1)  = "C:\Mod\diarios2001.prn"
      WeatherF(2)  = "C:\Mod\diarios2002.prn"
      WeatherF(3)  = "C:\Mod\diarios2003.prn"
      WeatherF(4)  = "C:\Mod\diarios2005.prn"
      WeatherF(5)  = "C:\Mod\diarios2006.prn"
      WeatherF(6)  = "C:\Mod\diarios2007.prn"
      WeatherF(7)  = "C:\Mod\diarios2009.prn"
      WeatherF(8)  = "C:\Mod\diarios2010.prn"
      WeatherF(9)  = "C:\Mod\diarios2011.prn"
      WeatherF(10) = "C:\Mod\diarios2013.prn"
      WeatherF(11) = "C:\Mod\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.50   !adimensional
      !IAF  = 3.5    !adimensional
      
      Open (unit=2, File="C:\Mod\Saida_Crop.out", status="replace")
      write (2,*) 'semeadura', DJSow
      write (2,*) 'colheita' , DJSow+ciclo
      write (2,*) 'Safra#     Yp         Yw            Yw/Yp '
     

100   Format (A200)   
101   Format (4I5, 6F12.3)  
      
      nd=1
      Do i=1,11
         Armaz = Umidade * ProfRad * 1000  !mm 
          
         Open (unit=1, File=WeatherF(i), Status="old")
         
         Read (1, 100) dumt

         !Lendo os dados do arquivos WeatherF
          Do DJ=1,365
             read (1,*) dum, dum, dum,Temp(nd),dum,dum, dum, dum, dum, dum, dum, Chuva(nd), Qg(nd),Rn(nd)
             nd = nd + 1
           Enddo
           Close (1)
      Enddo

      PPtotal = 0.
      PAtotal = 0.

      i = 0
      dia = 0
      
      safra = 1
      Do DJ= 1,4015
          i = i + 1
          If(i > 365+ciclo) Then
              dia  = dia + 365
              i = 1 + ciclo
              PpTotal = 0.
              PaTotal = 0.
              PARtotal = 0.
              das = 0.
          End if
          
          			
          !Effect of Temperature on Photosynthesis
			IF(temp(dj)<= Tinf .OR. temp(dj)>= Tsup) THEN
				ftemp = 0.
			ELSE IF (temp(dj)>= Tot1 .AND. temp(dj)<= Tot2) THEN
				ftemp = 1.
			ELSE IF (temp(dj)> Tinf .AND. temp(dj)< Tot1) THEN
				ftemp = (temp(dj)-Tinf)/(Tot1-Tinf)
			ELSE
				ftemp = (temp(dj)-Tot2)/(Tsup-Tot2)
            ENDIF
            
          If (DJ > (DJSow+dia).and.DJ<(DJSow+Ciclo+dia)) Then
          
            !Calculo ET por Priestley-Taylor
            das = das + 1
            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/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)
            
             !Estimando a interceptaçao de PAR para as diferentes fases fenologicas
             If (das < 30) Then
                 aPAR(DJ) = PAR(DJ) * (1-EXP(-K*IAF(1)))
             Elseif(das>30 .and. das<60) Then
                 aPAR(DJ) = PAR(DJ) * (1-EXP(-K*IAF(2)))
             Elseif(das>60 .and. das<100) Then
                 aPAR(DJ) = PAR(DJ) * (1-EXP(-K*IAF(3)))
             Else
                 aPAR(DJ) = PAR(DJ) * (1-EXP(-K*IAF(4)))
             Endif

             MS = aPAR(DJ) * RUE
             taxaPP = ms * ic * (1/(1-u))*ftemp* 10
             
             PpTotal  = PPTotal + taxaPP
             PaTotal  = PaTotal + taxaPP*FatorRed

             write (*,*) i,(DJSow+dia),(DJSow+Ciclo+dia),dj,taxaPP,PAR(DJ)
             If (DJ == (DJSow+Ciclo+dia)-1) Then
                write (2,102) safra, PpTotal,PaTotal,PaTotal/PpTotal
                safra = safra + 1
             End If
        End if
      Enddo
102   Format (1I5, 6F12.3) 
      
      write (*,*) 'Fim do processamento. Digite qualquer tecla para encerrar o programa.'
      read*
      Close(2)

end program