!  Rad_Vegetacao.f90 
!
!  FUNCTIONS:
!  Rad_Vegetacao - Entry point of console application.
!  Written by Clerverson Freitas - 4/4/2019


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

    program Rad_Vegetacao

    implicit none

    integer ano(730), doy(730), das, dj, NI, djsowing, ciclo, i, nd, j, nda
    real rue, k, ic, u, iaf, ms, taxapp, pptotal, ftemp, Tinf, Tot1, Tot2, Tsup, dum
    real PAR(730), apar(730), qg(730), chuva(730), rn(730), temp(730)
    character*100 title
    
    open(unit=1, file="Diarios2013_2014.prn", status="old")
    open(unit=2, file="Saida.out", status="replace")
    
    !perguntando ao usuario a data e o ciclo
    print *,'Informe a data de semeadura:'
    read *,djsowing
    
    print *,'Informe o ciclo da cultura:'
    read *,ciclo
    
    !definindo os valores dos parametros
    ic  = 0.33
    rue = 3.84
    u   = 0.13
    k   = 0.5
    iaf = 3.5
	
	!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)
  
100 format (a100)
  
    pptotal = 0.
    read (1,100) title


    !lendo o arquivo de dados climaticos
    do dj = 1, 730
        read (1,*) ano(dj), doy(dj), dum, temp(dj), dum, dum, dum, dum, dum, dum, dum, chuva(dj), qg(dj), rn(dj)
		!write(*,*) ano(dj), temp(dj)
	enddo
	
	das = 1
    do dj = 1, 730
        if (dj >= djsowing .and. dj<(djsowing+ciclo)) then
            Par(dj) = 0.5 * Qg(dj)
            aPAR(dj) = par(dj) * (1-exp(-k*iaf))
            ms = apar(dj) * rue
           
			!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
			
			!WRITE(*,*) dj, temp(dj), ftemp
				
			taxaPP = ms * ic * (1/(1-u))*ftemp* 10
            
			pptotal = pptotal + taxapp
            
			!if (das == (ciclo - 1)) then
                !write (2,102) das, djsowing, djsowing+ciclo, pptotal
                !write (*,102) das, djsowing, djsowing+ciclo, pptotal
            !    das = 1
            !    pptotal = 0.
            !endif
            write (*,102) das, djsowing, djsowing+ciclo, aPAR(dj),pptotal
			write (2,102) das, djsowing, djsowing+ciclo, aPAR(dj),pptotal
            das = das+1
        endif
    enddo

102 format (3i11, 4f12.2)
    close(1)
    close(2)
    !write (*,102) das, djsowing, djsowing+ciclo, pptotal
    !pause
	
 end program Rad_Vegetacao