	PROGRAM DEFL_PC
C	Programa contendo modelo de umidade do solo.

	REAL EP(25,12),XP(520),YP(520),XPP(25),YPP(25)
	REAL FCHUVA(500),CHUVA(500),CC(520),PP(520),I0(520),I0L
        REAL PERCOL(500),INFF1,INFI,INFF,UMIDADE(500)
	REAL ICMES(500),PESCOA(500),PEVAPO(500)
	CHARACTER*25 NOMEPOSTO(25),POSTO(520),DUMMY,SOLOS(520)

	BETA=1.

C	Abertura dos arquivos I/O
	OPEN(UNIT=10,FILE='evapo.ent')
	OPEN(UNIT=20,FILE='chuva.ent')
	OPEN(UNIT=30,FILE='solos.ent')
	OPEN(UNIT=60,FILE='defluvio.sai')
	OPEN(UNIT=70,FILE='prev.ent')
	OPEN(UNIT=80,FILE='umidade.sai')

C	Define algumas variaveis
	GLAT=111.1
	PI=3.141592
	WRITE(80,*)'UMIDADE DO SOLO (% CC)'

C	Leitura dos dados de evapo-transpiracao potencial.
	READ(10,*)NPOSTOS
	DO I=1,NPOSTOS
	   READ(10,100)NOMEPOSTO(I),IYG,IYM,IXG,IXM,(EP(I,J),J=1,12)
	   RLAT=FLOAT(IYG)+FLOAT(IYM)/60.
 	   RLONG=FLOAT(IXG)+FLOAT(IXM)/60.
 	   RRLAT=(RLAT*PI)/180.
	   GLONG=2.0*6378.*COS(RRLAT)*PI/360.
	   XPP(I)=-1.*GLONG*RLONG
	   YPP(I)=-1.*GLAT*RLAT
	ENDDO

C	Leitura dos dados de solo das estacoes.
	READ(30,*)NSOLOS
	DO II=1,NSOLOS
	   READ(30,410)SOLOS(II),CC(II),PP(II),I0(II)
	ENDDO	

C	Come‡o dos calculos para cada posto.
	READ(20,*)NMES,IDIAF

C	Come‡o do LOOP envolvendo todos os postos.
	DO I=1,300
	   READ(20,200)DUMMY,IYG,IYM,IXG,IXM
	   IF(DUMMY.EQ.'FIMDEARQUIVO             ') GOTO 50
	   POSTO(I)=DUMMY
	   WRITE(6,*)POSTO(I)
	   RLAT=FLOAT(IYG)+FLOAT(IYM)/60.
 	   RLONG=FLOAT(IXG)+FLOAT(IXM)/60.
 	   RRLAT=(RLAT*PI)/180.
	   GLONG=2.0*6378.*COS(RRLAT)*PI/360.
	   XP(I)=-1.*GLONG*RLONG
	   YP(I)=-1.*GLAT*RLAT
	   RLAT=RLAT*100.
	   RLONG=RLONG*100.

C	Leitura dos dados de chuva.
	   DO J=1,NMES
	      L=(J-1)*31
	      IF(J.NE.NMES) THEN
	         READ(20,300)(FCHUVA(K),K=L+1,L+31)
	       ELSE
	         READ(20,300)(FCHUVA(K),K=L+1,L+IDIAF)
	      ENDIF
	   ENDDO
	   NFDIAS=(NMES-1)*31+IDIAF
	   NDIAS=0
	   INICIAL=0
	   NNDIAS=0
	   NNMES=10
	   DO N=1,NFDIAS
	      NNDIAS=NNDIAS+1
	      IF(NNDIAS.GT.31) THEN
	         NNDIAS=1
	         NNMES=NNMES+1
	         IF(NNMES.GT.12) NNMES=1
	      ENDIF
c******************************************************************
c	      Esta parte e' opcional, para caso de falta de dados
c******************************************************************
	      if(fchuva(n).gt.990.) fchuva(n)=0.
c******************************************************************
	      IF(FCHUVA(N).NE.888.8) THEN
	         NDIAS=NDIAS+1
	         CHUVA(NDIAS)=FCHUVA(N)
	         ICMES(NDIAS)=NNMES
	      ENDIF
	   ENDDO
	   INICIAL=93

C	Definicao dos parametros de solo do posto em questao.
	   CCL=999.
	   DO J=1,NSOLOS
	      IF(POSTO(I).EQ.SOLOS(J)) THEN
	         CCL=CC(J)
	  	 PPL=PP(J)
		 I0L=I0(J)
	      ENDIF
	   ENDDO

C	Evita calculos desnecess rios.
	IF(CCL.EQ.999.) GOTO 53
	
C	Determina‡Ůo da estacao de evapotranspiracao potencial a ser 
c       utilizadas no calculo da evp de cada posto.
	   DMIN=10000000.
	   DO 13 L=1,NPOSTOS
	      DIST=SQRT((XP(I)-XPP(L))**2.+(YP(I)-YPP(L))**2.)
	      IF(DIST.LT.DMIN) THEN
		   DMIN=DIST
		   IDMIN=L
	      ENDIF
  13      CONTINUE

C      Comeco do calculo da umidade do solo e demais parcelas do balanco
	  EVAPTOTAL=0.
	  PERTOTAL=0.
	  ESCOATOTAL=0.
	  CHUVATOTAL=0.
	  HSI=0.
	  ERRO=0.2

C	Come‡o do LOOP para todos os dias do posto em questao.
	  DO K=1,NDIAS
	
C	Calculo da evapotranspira‡Ůo potencial para o dia.
	     EVAPO=(0.7/30.)*EP(IDMIN,ICMES(K))

C 	Calculo dos valores iniciais.
	     HSF=HSI+CHUVA(K)
	     INFI=FUNCINFILT(HSI,PPL,I0L,CCL,BETA)
	     INFF=FUNCINFILT(HSF,PPL,I0L,CCL,BETA) 

C 	Comeco dos varios casos de infiltracao.
	     ICASO=0

C 	*********CASO 1***********************
    	     IF(INFI.GE.CHUVA(K).AND.INFF.GE.CHUVA(K)) THEN
	        ESCOA=0.
	        PERI=FUNCPER(HSI,PPL,CCL,BETA)
	        PERF=FUNCPER(HSF,PPL,CCL,BETA)
	        PERT=(PERI+PERF)/2.
		IF(PERT.GE.HSF) PERT=HSF
	        HSF=HSF-PERT
	        HSM=(HSI+HSF)/2.
	        EVAPORA=FUNCEVAP(HSM,HSF,EVAPO,CCL,BETA)
	        HSF=HSF-EVAPORA
	  	ICASO=1
	     ENDIF

C 	*********CASO 2***********************
  	      IF(INFI.GT.CHUVA(K).AND.INFF.LT.CHUVA(K)
     *					.AND.HSF.LT.CCL) THEN
  1	         X=(CHUVA(K)-INFF)/(INFI-INFF)
		 VESCOA=X*(CHUVA(K)-INFF)/2.
		 IF(VESCOA.GT.CHUVA(K)) VESCOA=CHUVA(K)
		 HSF1=HSF-VESCOA
	         INFF1=FUNCINFILT(HSF1,PPL,I0L,CCL,BETA)
	         DESVIO=ABS(INFF1-INFF)
		 IF(DESVIO.GE.ERRO) THEN
		     INFF=INFF1
		     GOTO 1
		 ENDIF
	         ESCOA=VESCOA
	         IF(ESCOA.GT.CHUVA(K)) ESCOA=CHUVA(K)
		 HSF=HSF-ESCOA
	         PERI=FUNCPER(HSI,PPL,CCL,BETA)
	         PERF=FUNCPER(HSF,PPL,CCL,BETA)
	         PERT=(PERI+PERF)/2.
		 IF(PERT.GE.HSF) PERT=HSF
		 HSF=HSF-PERT
		 HSM=(HSI+HSF)/2.
	         EVAPORA=FUNCEVAP(HSM,HSF,EVAPO,CCL,BETA)
	         HSF=HSF-EVAPORA
	         ICASO=2
	      ENDIF

C 	*********CASO 3***********************
  	      IF(INFI.GT.CHUVA(K).AND.INFF.LT.CHUVA(K)
     *					.AND.HSF.GE.CCL) THEN
		 TERMO1=CCL-HSI
		 TERMO2=(0.5+(INFI-CHUVA(K))/(2.*(INFI-PPL)))*
     * 		        (CHUVA(K)-PPL)
		 TERMO3=PPL+TERMO2
	         X=TERMO1/TERMO3
		 IF(X.GT.1.) THEN
		    INFF=PPL
  3	            XX=(CHUVA(K)-INFF)/(INFI-INFF)
		    VESCOA=XX*(CHUVA(K)-INFF)/2.
		    IF(VESCOA.GT.CHUVA(K)) VESCOA=CHUVA(K)
		    HSF1=HSF-VESCOA
	            INFF1=FUNCINFILT(HSF1,PPL,I0L,CCL,BETA)
	            DESVIO=ABS(INFF1-INFF)
		    IF(DESVIO.GE.ERRO) THEN
		        INFF=INFF1
		        GOTO 3
		    ENDIF
	            ESCOA=VESCOA
	            IF(ESCOA.GT.CHUVA(K)) ESCOA=CHUVA(K)
		    HSF=HSF-ESCOA
	            PERI=FUNCPER(HSI,PPL,CCL,BETA)
	            PERF=FUNCPER(HSF,PPL,CCL,BETA)
	            PERT=(PERI+PERF)/2.
		    IF(PERT.GE.HSF) PERT=HSF
		    HSF=HSF-PERT
		    HSM=(HSI+HSF)/2.
	            EVAPORA=FUNCEVAP(HSM,HSF,EVAPO,CCL,BETA)
	            HSF=HSF-EVAPORA
	            ICASO=32
	          ELSE
		    ESCOA=CHUVA(K)-TERMO1-(1.-X)*PPL
	            IF(ESCOA.GT.CHUVA(K)) ESCOA=CHUVA(K)
		    HSF=HSF-ESCOA
	            PERI=FUNCPER(HSI,PPL,CCL,BETA)
		    PERT=(PERI+PPL)*X/2.+PPL*(1.-X)
		    IF(PERT.GE.HSF) PERT=HSF
		    HSF=HSF-PERT
		    HSM=(HSI+HSF)/2.
	            EVAPORA=FUNCEVAP(HSM,HSF,EVAPO,CCL,BETA)
	            HSF=HSF-EVAPORA	     
		    ICASO=3
		 ENDIF
 	      ENDIF

C 	*********CASO 4***********************
  	      IF(INFI.LT.CHUVA(K).AND.INFF.LT.CHUVA(K)
     *                                  .AND.(HSI+INFI).LT.CCL) THEN
  2		 HSF1=HSI+(INFI+INFF)/2.
	         INFF1=FUNCINFILT(HSF1,PPL,I0L,CCL,BETA)
	         DESVIO=ABS(INFF1-INFF)
		 IF(DESVIO.GE.ERRO) THEN
		     INFF=INFF1
		     TERMO1=0.
		     GOTO 2
		 ENDIF
		 HSF=HSI+(INFI+INFF)/2.
		 ESCOA=CHUVA(K)-(INFI+INFF)/2.
		 PERI=FUNCPER(HSI,PPL,CCL,BETA)
	         PERF=FUNCPER(HSF,PPL,CCL,BETA)
	         PERT=(PERI+PERF)/2.
		 IF(PERT.GE.HSF) PERT=HSF
		 HSF=HSF-PERT
		 HSM=(HSI+HSF)/2.
	         EVAPORA=FUNCEVAP(HSM,HSF,EVAPO,CCL,BETA)
	         HSF=HSF-EVAPORA
	   	 ICASO=4
	      ENDIF

C 	*********CASO 5***********************
  	      IF(INFI.LT.CHUVA(K).AND.INFF.LT.CHUVA(K)
     *                                  .AND.HSF.GE.CCL) THEN
		 TERMO1=CCL-HSI
		 IF(INFI.EQ.0.) THEN
		    X=0.
		    GOTO 72
		 ENDIF
		 TERMO2=(INFI-PPL)/2.+PPL
		 X=TERMO1/TERMO2
 72		 IF(X.GT.1.) THEN
		     INFF=PPL
  4		     HSF1=HSI+(INFI+INFF)/2.
	             INFF1=FUNCINFILT(HSF1,PPL,I0L,CCL,BETA)
	             DESVIO=ABS(INFF1-INFF)
		     IF(DESVIO.GE.ERRO) THEN
		       INFF=INFF1
		       GOTO 4
		     ENDIF
		     HSF=HSI+(INFI+INFF)/2.
		     ESCOA=CHUVA(K)-(INFI+INFF)/2.
		     PERI=FUNCPER(HSI,PPL,CCL,BETA)
	             PERF=FUNCPER(HSF,PPL,CCL,BETA)
	             PERT=(PERI+PERF)/2.
		     IF(PERT.GE.HSF) PERT=HSF
		     HSF=HSF-PERT
		     HSM=(HSI+HSF)/2.
	             EVAPORA=FUNCEVAP(HSM,HSF,EVAPO,CCL,BETA)
	             HSF=HSF-EVAPORA
	   	     ICASO=54
		   ELSE
		     ESCOA=CHUVA(K)-(1.-X)*PPL-TERMO1
	             IF(ESCOA.GT.CHUVA(K)) ESCOA=CHUVA(K)
		     HSF=HSF-ESCOA
		     PERI=FUNCPER(HSI,PPL,CCL,BETA)
		     PERT=(PERI+PPL)*X/2.+PPL*(1.-X)
   		     IF(PERT.GE.HSF) PERT=HSF
		     HSF=HSF-PERT
		     HSM=(HSI+HSF)/2.
	             EVAPORA=FUNCEVAP(HSM,HSF,EVAPO,CCL,BETA)
	             HSF=HSF-EVAPORA
		     ICASO=5
		ENDIF
	      ENDIF
	      INFF=FUNCINFILT(HSF,PPL,I0L,CCL,BETA)

C	Imprime resultados parciais para cada dia.
	   UMIDADE(K)=HSI
	   PERCOL(K)=PERT
	   PEVAPO(K)=EVAPORA
	   PESCOA(K)=ESCOA
	   IF(K.GE.INICIAL) THEN
	      WRITE(60,800)K-INICIAL+1,HSI,CHUVA(K),ESCOA,PERT,EVAPORA,
     *        HSF,ICASO,INFI,INFF,EVAPO
C	      Atualiza valores totais.
	      ESCOATOTAL=ESCOATOTAL+PESCOA(K)
	      PERTOTAL=PERTOTAL+PERCOL(K)
	      EVAPTOTAL=EVAPTOTAL+PEVAPO(K)
	      CHUVATOTAL=CHUVATOTAL+CHUVA(K)
	    ENDIF

C	Seta valor de umidade inicial para novo dia.
	    HSI=HSF
  	  ENDDO
	  UMIDADE(NDIAS+1)=HSF
	  SIUMIDADE=(HSF/CCL)*100.

C	 Imprime resultados gerais para o posto e come‡a um novo posto.	
	  WRITE(80,150)POSTO(I),NINT(SIUMIDADE)
	  WRITE(60,140)POSTO(I),CHUVATOTAL,ESCOATOTAL,PERTOTAL,EVAPTOTAL
	  WRITE(70,120)POSTO(I),RLAT,RLONG,HSF,CCL,EVAPO,PPL
  53	ENDDO
  50	CONTINUE


 100	FORMAT(A25,2X,2I2,2X,2I2,12F8.3) 
 110	FORMAT(2I4,I6)
 120    FORMAT(A25,2F7.1,4F10.4)
 130	FORMAT(A25,2F7.1)
 140	FORMAT(A25,4F7.2)
 150    FORMAT(A25,2X,I4)
 200	FORMAT(A25,2X,2I2,2X,2I2)
 300	FORMAT(31F5.1)
 400 	FORMAT(A25,F6.1)
 410	FORMAT(A25,3f6.1)
 800	FORMAT(I3,15X,6F7.2,2X,I2,3F7.2)
 900	FORMAT('TOTAL',6X,5f9.3) 
	STOP
	END
	   

C	Definicao da infiltracao em funcao da umidade (HS~X).
	FUNCTION FUNCINFILT(X,PPL,I0L,CCL,BETA)
	REAL I0L
	IF(X.LE.CCL) THEN
	    FUNCINFILT=(I0L-PPL)*((CCL-X)/CCL)**BETA+PPL
	  ELSE
	    FUNCINFILT=PPL
	ENDIF
	RETURN
	END

C	Definicao da percolacao em funcao da umidade (HS~X).
	FUNCTION FUNCPER(X,PPL,CCL,BETA)
	IF(X.LE.CCL) THEN
	    FUNCPER=PPL*(X/CCL)**BETA
	  ELSE
	    FUNCPER=PPL
	ENDIF
	RETURN
	END
	
C	Definicao da evapotranspiracao em funcao da umidade (HS~X).
	FUNCTION FUNCEVAP(X,Y,EVAPO,CCL,BETA)
	IF(X.GT.(0.7*CCL)) THEN
	    FUNCAUX=EVAPO
	  ELSE
	    FUNCAUX=EVAPO*(X/(CCL*0.7))**BETA
	ENDIF
	IF(FUNCAUX.GE.Y) THEN
	    FUNCEVAP=Y
	  ELSE
	    FUNCEVAP=FUNCAUX
	ENDIF
	RETURN
	END

