C**********************************************************************C
C                                                                      C
C     *** SUBROUTINE WARMPET ***                                       C
C     INPUT  MONTHLY PRECIPITATION AND TEMPERATURE, AND CALCULATE      C
C     MONTHLY PET AND POTENTIAL AET                                    C
C                                                                      C
C**********************************************************************C
      SUBROUTINE WARMPET(I, J, M, MNDAY)
      
      COMMON/BASIC/NGRID,NYEAR,NLC,BYEAR,IYSTART,IYEND,POP_FLAG,REGION,
     &CLIMATE_FLAG
       
      COMMON/CELLINFO/LADUSE(4000,10),HUCNO(4000),
     >  LATUDE(4000),LONGI(4000),TAREA(4000)
     
      COMMON/LANDCHANGE/FPERD, FOR_CONV_TO, LAI_CHG
              
      COMMON/OUTPUT/ PET(200,12,10),APET(12),PAET(200,12,10),APAET(12),
     &AET(12), AVUZTWC(12), AVUZFWC(12), AVLZTWC(12), AVLZFPC(12),
     &AVLZFSC(12),RO(2200,200,12)            
      
      COMMON/CLIMATE/ RAIN(4000,200,12), TEMP(4000,200, 12)
      
      COMMON/LAI/LAI(4000,12,10)   
                
      INTEGER NLC, HUCNO(4000), M, BYEAR, WRR
      
      INTEGER MONTHD(12),MONTHL(12), MJD(12), MMD
      
      REAL DEGLAT, LATUDE(4000), PET(200, 12, 10),TAREA(4000)
 
      REAL APET(12), PE, LADUSE(4000,10), TEMP(4000,200,12)
      
      REAL TPET, RAIN(4000,200,12), FOREST
      
      REAL FPERD
      
      REAL LAI(4000,12,10)
      
      REAL PAET(200,12,10), TPAET, APAET(12)
      
      CHARACTER*5 CLIMATE_FLAG
      
      CHARACTER*2 REGION
     
             
C --- Number of days for each month during regular year
      DATA MONTHD/31,28,31,30,31,30,31,31,30,31,30,31/
C 
C --- Number of days for each month during leap year
      DATA MONTHL/31,29,31,30,31,30,31,31,30,31,30,31/
C-----JULIAN DAY FOR MID-DAY OF EACH MONTH
      DATA MJD/15,46,76,107,137,168,198,229,259,290,321,351/      

C --- Calculate Monthly potential evapotranspiration
                 
        
            DTEMP=TEMP(I,J,M)
            
            TPET=0. 
            
            TPAET=0.         
                         
        DO 40 K=1, NLC
        
              
C -- MMD = JULIAN DATE FOR MONTH M
                               
            MMD=MJD(M)
                     
            DEGLAT = LATUDE(I)
            
C -K= LANDUSE TYPE, DTEMP = AIR TEMP, HPEC = CORRECTION PARAMETER, PE =CALCUALTED PET (MM)
           
            CALL HAMON(K,DTEMP,M,MMD,DEGLAT,PE)
            
            
C ----PET (YEAR, MONTH, LANDUSE) FOR THE CURRENT CELL)
           
             PET(J,M,K) = PE  
                  
             PET(J,M,K) = PET(J,M,K) * MNDAY 

C ----CALCULATE PAET 
C ----PAET (YEAR, MONTH, LANDUSE) FOR THE CURRENT CELL)
C ----PAET IS THE POTENTIAL AET, ASSUMING SOIL MOISTURE NOT LIMITING

             IF (REGION.EQ.'US') THEN

                WRR=INT(HUCNO(I)/1000000)
             
C --- COMPUTE THE TOTAL FOREST LAND COVER PERCENTAGE (PRIOR TO CHANGE BY USER)
C --- THIS IS NECESSARY TO PREVENT USING BOTH EQUATIONS FOR THE SAME HUC
C --- FOR ORIGINAL FOREST COVER AND REDUCED FOREST COVER IF REDUCING THE FOREST
C --- COVER CAUSES FOREST TO GO FROM >0.2 TO <0.2

C -- SET PAET FOR OPEN WATER LAND COVER CLASS TO PET

                IF(K.EQ.8) THEN
                
                PAET(J,M,K)=PET(J,M,K)
                
                ELSE

C --- COMPUTE THE TOTAL FOREST LAND COVER PERCENTAGE (PRIOR TO CHANGE BY USER)
C --- THIS IS NECESSARY TO PREVENT USING BOTH EQUATIONS FOR THE SAME HUC
C --- FOR ORIGINAL FOREST COVER AND REDUCED FOREST COVER IF REDUCING THE FOREST
C --- COVER CAUSES FOREST TO GO FROM >0.2 TO <0.2

                FOREST=(LADUSE(I,2)+LADUSE(I,3)+LADUSE(I,4))/(1+FPERD)
             
                IF (WRR.EQ.1.OR.(WRR.EQ.2.AND.LATUDE(I).GT.40).OR.
     &          WRR.EQ.4.OR.WRR.EQ.5) THEN
             
                   IF (FOREST.GT.0.2) THEN
     
                      PAET(J,M,K) = 0.00169*PET(J,M,K)*RAIN(I,J,M) + 
     &                0.4*PET(J,M,K)+7.83*LAI(I,M,K)
     
                   ELSE
             
                     PAET(J,M,K) = (PET(J,M,K))*LAI(I,M,K)*0.0222+0.174*
     &               RAIN(I,J,M)+0.502*(PET(J,M,K))+5.31*LAI(I,M,K)  
     
                   ENDIF
                
                ELSE
             
                   PAET(J,M,K) = (PET(J,M,K))*LAI(I,M,K)*0.0222+0.174*
     &             RAIN(I,J,M)+0.502*(PET(J,M,K))+5.31*LAI(I,M,K)                           

                ENDIF
                
                ENDIF
             
             ELSEIF (REGION.EQ.'MX'.OR.REGION.EQ.'RW') THEN
             
                PAET(J,M,K) = (PET(J,M,K))*LAI(I,M,K)*0.0222+0.174*
     &          RAIN(I,J,M)+0.502*(PET(J,M,K))+5.31*LAI(I,M,K)   
     
             ENDIF   

             
         


C ----CALCULATE TOTAL PET AND PAET FOR THE HUC FOR A GIVEN YEAR AND MONTH
        
             
             TPET = TPET + PET(J,M,K) * LADUSE(I,K)/TAREA(I)
             
             TPAET = TPAET + PAET(J,M,K) * LADUSE(I,K)/TAREA(I)
             

   40      CONTINUE                 

C ------APET =AVERAGE PET FOR CURRENT CELL, FOR ALL YEAR, MONTH

              APET(M) = TPET    
              
C ------APAET =AVERAGE PAET FOR CURRENT CELL, FOR ALL YEAR, MONTH

              APAET(M) = TPAET    

C ------TEST OUTPUT

C      WRITE(77, 5040) HUCNO(I), J, M, RAIN(I,J,M), APET(M), APAET(M)
      
5040  FORMAT (3I10, 3F20.1)
                              
C --- Return
      RETURN
      END
C**********************************************************************C
C                                                                      C
C     *** SUBROUTINE HAMON ***                                         C
C     Calculate potential evapotranspiration by using Hamon's          C
C     equation                                                         C
C                                                                      C
C**********************************************************************C
      SUBROUTINE HAMON(K, TEMP,MONTH,J,DEGLAT,PE)

      
      REAL PE, TEMP, DEGLAT
      
      REAL SOLDEC, SSANG, DAY
      
      REAL ESAT, RHOSAT, PI
      
      INTEGER J
      
      PI = 3.14159265

C --- ESAT = saturated vapor pressure at TEMP
C --- RHOSAT = saturated vapor density at TEMP
      
      ESAT = 6.108*EXP(17.2693882*TEMP/(TEMP+237.3))
      RHOSAT = 216.7*ESAT/(TEMP+273.3)
      
C --- CALCULATE MEAN MONTHLY DAY LENGTH (DAY) BASED ON LATITUDE AND MID-MONTH JULIAN DATE
C --- SHUTTLEWORTH, W.J. 1993. Evaporation, in Handbook of Hydrology, 
C --- MAIDMENT, D.R. ED., MCGRAW HILL, NY, PP. 4, 1-53

C ---    CALCULATE THE ANGLE OF SOLAR DECLINATION IN RADIANS     
 
         SOLDEC = 0.4093*SIN(2*PI*J/365.-1.405)
      
C ---    CALCULATE THE SUNSET HOUR ANGLE IN RADIANS

         SSANG = ACOS(-1*TAN(DEGLAT*.0174529)*TAN(SOLDEC))
      
C ---    CALCULATE THE ADJUSTMENT IN LENGTH OF DAY FOR 12 HR PERIOD
      
         DAY = 2*SSANG/PI

C --- Calculate Daily PE
     
      PE = 0.1651*DAY*RHOSAT
      


      RETURN

      END



