program evap
c**********************************************************************
c  This program computes evapotranspiration from methods described in *
c  Thornthwaite and Mather 1957, using tables from their report and   *
c  monthly mean temperatures that were thissen weighted for individual*
c  subbasins at the Great Lakes Environmental Research Laboratory.    *
c  Designed and Programmed by:  Cynthia E. Sellinger, Physical        *
c  Scientist, Great Lakes Environmental Research Laboratory, National *
c  Oceanic and Atmospheric Administration, Department of Commerce.    *
c                                                                     *
c  20 Aug 2013 - Slight modifications by Tim Hunter, IT Specialist,   *
c  NOAA/GLERL.  In response to user-reported issues, I noted and      *
c  made minimal fixes including the following:                        *
c    1) Replaced imbedded TAB characters with spaces.  This was       *
c       required for consistent compiler performance.                 *
c    2) Fixed some formatting issues in the main program. Significant *
c       issues remain with the indentation. This code is not very     *
c       easy to read due to the indentation issues, but I do not      *
c       have time to finish repairing that.                           *
c    3) Very minor changes in file handling to make it a tiny bit     *
c       more user-friendly.                                           *
c    4) Recompiled with gfortran 4.8.0 and placed revised executable  *
c       on our ftp site.  Hoping that this resolves the user issues.  *
c**********************************************************************

      implicit none
      real avg,temp,ae,stt,
     +  heat1,year_heat,total,pet,evapo,lat,correction,
     +  heat,precip,temp_1,precip_1,latitude,p_pe,
     +  x,y
      integer i,j,year,col,mon,count,count2,
     +  year2,ns,year_1,mon_1,moist,row
      character*50 file1,file3
c      
c
      PRINT *, ' '
      PRINT *, ' '
      PRINT *, ' '
      PRINT *, ' '
      PRINT *, ' '
      print *, '         Enter the filename of the monthly '
      print *, '         averaged temperature and precipitation.  '
      print *, '   NOTE: temperatures should be in degrees of'
      print *, '         Celsius and precipitation values should  '
      print *, '         be monthly totals in units of millimeters.'
      PRINT *, ' '
      PRINT *, ' '
      PRINT *, ' '
      PRINT *, ' '
      PRINT *, ' '
      read(5,2)file1
      print *, 'Enter outfile name  '
      read(5,2)file3
c
      print *, 'enter "1" for northern hemisphere and "2" for '
      print *, 'southern hemisphere '
      read(5,*)ns
      print *, 'enter the latitude of your area  '
      read(5,*)latitude
      print *, '  '
      print *, '  '
      print *, '  '
      print *, ' Choose a number 1-11 for the moisture holding'
      print *, ' capacity of the soil '
      print *, '  '
      print *, '  '
      print *, '  1.   25 millimeters'
      print *, '  2.   50 millimeters'
      print *, '  3.   75 millimeters'
      print *, '  4.  100 millimeters'
      print *, '  5.  125 millimeters'
      print *, '  6.  150 millimeters'
      print *, '  7.  200 millimeters'
      print *, '  8.  250 millimeters'
      print *, '  9.  300 millimeters'
      print *, ' 10.  350 millimeters'
      print *, ' 11.  400 millimeters'
      print *, '  '
      print *, '  '
      print *, '  '
      read(5,*)moist
c
c
c
c Three steps must be taken to find evapotranspiration, first must obtain 
c the heat index given monthly temperatures and using table 2 of section 1.
c
c  read in heat index,unadjusted potential evap files, and correction factor
c  files
c
c
c open files
c
      open(1,file=file1,status='old')
      open(2,status='scratch',access='direct',recl=16,
     +  form='unformatted')
      open(3,file=trim(file3),status='replace')
      open(4,status='scratch',access='direct',recl=8,
     +  form='unformatted')
      open(7,status='scratch',access='direct',recl=4,
     +  form='unformatted')
c
c
c
  2   format(a)
c
      count = 0
      do i = 1,50000
       read(1,*,end=784)year,mon,avg,precip
       count = count + 1
       if(avg.lt.0) then
         avg = 0.0
       endif
       write(2,rec=count)year,mon,avg,precip
      enddo
 784  continue
      close(1)
      print *, 'number of records = ', count
c
c   look up the heat index
c
      total = 0.0
      count2 = 0
      do i = 1,count
        read(2,rec=i)year,mon,temp,precip
c
        print *, 'computing the heat index for record #',i,'     of ',
     +    count
c        call index(temp,heat1,row,col,x,y)
        call index(temp,heat1,row,col)
c
c
  100  format(1x,4i6,5f10.2)
c
        if(i.gt.1) then
           read(2,rec=i-1)year_1,mon_1,temp_1,precip_1
           if(mon.gt.mon_1.and.year.eq.year_1) then
              total = total + heat1
              if(i.eq.count) then
                 count2 = count2 + 1
                 year_heat = total
                 year2 = year_1
                 write(4,rec=count2)year2,year_heat
                 total = 0.0 + heat1
              endif
           else
              count2 = count2 + 1
              year_heat = total
              year2 = year_1
              write(4,rec=count2)year2,year_heat
              total = 0.0 + heat1
            endif
        elseif(i.eq.1) then
           total = 0.0 + heat1
        endif
      enddo
c
c
c  find unadjusted evapotranspiration by using the yearly summed heat index
c  and the monthly air temperature to apply to table 4
c 
c
c
      do i = 1,count
c
         read(2,rec=i)year,mon,temp,precip
c
c         print *, 'count2 = ',count2
         do j = 1,count2
           read(4,rec=j)year2,year_heat
           if(year2.eq.year) then
             heat = year_heat
c
             if(temp.ge.0.and.temp.lt.26.5) then
               call table4(temp,heat,pet)
c
             elseif(temp.ge.26.5) then
               call table5(temp,pet)
             endif
c
c   adjust evapotranspiration for latitude and month of the year by multiplying
c   by the adjustment factor given in table 6
c


c   adjust evapotranspiration for latitude and month of the year by multiplying
c   by the adjustment factor given in table 6
c
             print *, 'performing evap lat correction for rec # ',i,
     +           '  of',count
             lat = latitude
             col = mon
             if(ns.eq.1) then
                call table6(lat,col,correction)
             else
                call table7(lat,col,correction)
             endif
c
             evapo = pet * correction
c
           endif
         enddo
         write(7,rec=i)evapo
      enddo
      
c
c  compute actual evapotranspiration
c
      call moisture(count,moist)
      
c
      write(3,16)
      do i = 1,count
         read(2,rec=i)year,mon,temp,precip
         read(7,rec=i)evapo
         read(14,rec=i)ae,stt
         read(13,rec=i)p_pe
         write(3,15)year,mon,precip,temp,evapo,ae
      enddo
 15   format(2x,i4,3x,i4,2x,4f8.2)
 16   format(2x,'Year',3x,'Month',4x,'Precip',3x,'Temp',
     +       3x,'  PET',2x,'Actual ET',/)
c
      close(2)
      close(3)
      close(4)
      close(7)
      close(13)
      close(14)
      stop
      end
c
c
c
c
c
c this subroutine reads table2 and returns a heat index
c
c      subroutine index(temp,heat1,row,col,x,y)
      subroutine index(temp,heat1,row,col)
      implicit none
      real heat(40,10),x,heat1,temp
      integer row,col,i,j,dummy,y
      
c
      open(1,file='table2.dat',status='old',recl=135)
      do i = 1,40
        read(1,*)dummy,(heat(i,j),j=1,10)
      enddo
c
c
c     
      y = int(temp)
      x = (temp - float(y)) * 10.
      row = y + 1
      col = int(x) + 1
      heat1 = heat(row,col)
c
c
      close(1)
      return
      end
c
c
      subroutine table4(temp,heat,pet)
      implicit none
      real evp(107,47),temp,heat,pet,index(1,47),temperature(107),
     +temp1,temp2,heat1,heat2
      integer i,j
c
c
c
c select table
c
      open(1,file='table4.dat',status='old')
c
      read(1,10)(index(1,j),j=1,47)
c
      do i = 1,107
      read(1,11)temperature(i),(evp(i,j),j=1,47)
      enddo
c
c   select the right pet for give heat index and temperature
c
      do i = 1,107
        do j = 1,47
c
      if(i.eq.1.and.j.eq.1) then
         if(temp.ge.temperature(i).and.temp.lt.temperature(i+1).
     +        and.heat.ge.index(1,j).and.heat.lt.index(1,j+1)) then
c  compute the nearest value
                     temp1 = temp - temperature(i)
           temp2 = temperature(i+1) - temp
           heat1 = heat - index(1,j)
           heat2 = index(1,j+1) - heat
             if(temp1.le.temp2.and.heat1.le.heat2) then
               pet = evp(i,j)
             elseif(temp1.le.temp2.and.heat1.gt.heat2) then
               pet = evp(i,j+1)
             elseif(temp1.gt.temp2.and.heat1.gt.heat2) then
               pet = evp(i+1,j+1)
             elseif(temp1.gt.temp2.and.heat1.le.heat2) then
               pet = evp(i+1,j)
             endif
                endif
c
      elseif(j.eq.1.and.i.gt.1) then
             if(temp.lt.temperature(i).and.temp.ge.temperature(i-1).
     +       and.heat.ge.index(1,j).and.heat.lt.index(1,j+1)) then
c  compute the nearest value
                    temp1 = temp - temperature(i-1)
          temp2 = temperature(i) - temp
          heat1 = heat -index(1,j)
          heat2 = index(1,j+1) - heat
            if(temp1.le.temp2.and.heat1.le.heat2) then
              pet = evp(i-1,j)
            elseif(temp1.le.temp2.and.heat1.gt.heat2) then
              pet = evp(i-1,j+1)
            elseif(temp1.gt.temp2.and.heat1.gt.heat2) then
              pet = evp(i,j+1)
            elseif(temp1.gt.temp2.and.heat1.le.heat2) then
              pet = evp(i,j)
            endif
        else
        endif        
c
      elseif(j.gt.1.and.i.eq.1) then
        if(temp.ge.temperature(i).and.temp.lt.temperature(i+1).
     +       and.heat.lt.index(1,j).and.heat.ge.index(1,j-1)) then
c  compute the nearest value
                  temp1 = temp - temperature(i)
        temp2 = temperature(i+1) - temp
        heat1 = heat - index(1,j-1)
        heat2 = index(1,j) - heat
         if(temp1.le.temp2.and.heat1.le.heat2) then
          pet = evp(i,j-1)
         elseif(temp1.le.temp2.and.heat1.gt.heat2) then
          pet = evp(i,j)
         elseif(temp1.gt.temp2.and.heat1.gt.heat2) then
          pet = evp(i+1,j)
         elseif(temp1.gt.temp2.and.heat1.le.heat2) then
          pet = evp(i+1,j-1)
         endif
        endif
c
      elseif(i.gt.1.and.j.gt.1) then
        if(temp.lt.temperature(i).and.temp.ge.temperature(i-1).
     +       and.heat.lt.index(1,j).and.heat.ge.index(1,j-1)) then
c
c  compute the nearest value
                 temp1 = temperature(i) - temp
            temp2 = temp - temperature(i-1)
       heat1 = index(1,j) - heat
       heat2 = heat - index(1,j-1)
c
                       if(temp1.le.temp2.and.heat1.le.heat2)  then
               pet = evp(i,j)
             elseif(temp1.le.temp2.and.heat1.gt.heat2) then
               pet = evp(i,j-1)
             elseif(temp1.gt.temp2.and.heat1.gt.heat2) then
               pet = evp(i-1,j-1)
             elseif(temp1.gt.temp2.and.heat1.le.heat2) then
               pet = evp(i-1,j)
             endif
        else
        endif
c
       endif
      enddo
      enddo
  10  format(9x,47f8.2)
  11  format(1x,f6.2,2x,47f8.2)
      close(1)
      return                    
      end
c
c
c
      subroutine table5(temp,pet)
      implicit none
      real temp,evp(13,10),x,pet
      integer row,col,i,j,dummy,y,temp2(13)
      
c
      data temp2/26,27,28,29,30,31,32,33,34,35,46,37,38/
c
      open(1,file='table5.dat',status='old',recl=128)
      read(1,2)
      do i = 1,12
      read(1,*)dummy,(evp(i,j),j=1,10)
      enddo
c
c
   2  format(a)
c   
      y = int(temp)
      do i = 1,12
      if(y.eq.temp2(i)) then
      row = i
      else
      endif
      enddo
c
      x = (temp - float(y)) * 10.
      col = nint(x + 1)
      pet = evp(row,col)
c
c
      close(1)
      return
      end
c      
c
      subroutine table6(lat,col,correction)
      implicit none
      real corr(51,12),correction,lat
      integer row,col,dummy,i,j
c
c
      open(1,file='table6.dat',status='old',recl=125)
c
      do i = 1,51
      read(1,*)dummy,(corr(i,j),j=1,12)
      enddo
c
      row = nint(lat) + 1
      correction = corr(row,col)
      close(1)
      return
      end
c      
      subroutine table7(lat,col,correction)
      implicit none
      real corr(51,12),correction,lat
      integer row,col,dummy,i,j
c
c
      open(1,file='table7.dat',status='old',recl=125)
c
      do i = 1,51
      read(1,*)dummy,(corr(i,j),j=1,12)
      enddo
c
      row = nint(lat) + 1
      correction = corr(row,col)
      close(1)
      return
      end
c      

c
c
      subroutine moisture(c2,moist)
c
c  this program reads moisture of thornwaite's soil moisture retention
c   and provides values for actual evapotranspiration
c
      implicit none
      real pet,precip,p_pe,act_pot_wl(3000),
     + delta_st(3000),ae,sum,sum2,sum_pos,
     + pe1(1500),w_ret(200,10),w_ret2,remainder,sum_neg,pe1_2,
     + w_ret3(200),water_loss,st2,sum7,
     + ae1,st(3000),stt,dummy
       integer i,j,k,col,row,start,finish,begin,check,c2,begin2,
     + c3,c4,b_c(200),year,mon,e_c(200),c5,year_1,mon_1,
     + c6,years,c8,c7,moist,c10,count,c9,c99
c
      open(13,status='scratch',access='direct',recl=4,
     +     form='unformatted')
      open(14,status='scratch',access='direct',recl=8,
     +     form='unformatted')
c
c compute p_pe
c
      do i = 1,c2
      read(2,rec=i)year,mon,dummy,precip
      read(7,rec=i)pet
      p_pe = precip - pet
      write(13,rec=i)p_pe
      enddo
c

c
      if(moist.eq.1) then
      open(1,file='table23.dat',status='old')
      c10 = 10
      elseif(moist.eq.2) then
      open(1,file='table24.dat',status='old')
      c10 = 23
      elseif(moist.eq.3) then
      open(1,file='table25.dat',status='old')
       c10 = 37
      elseif(moist.eq.4) then
      open(1,file='table26.dat',status='old')
      c10 = 51
      elseif(moist.eq.5) then
      open(1,file='table27.dat',status='old')
      c10 = 35
      elseif(moist.eq.6) then
      open(1,file='table28.dat',status='old')
      c10 = 75
      elseif(moist.eq.7) then
      open(1,file='table29.dat',status='old',recl=140)
      c10 = 106
      elseif(moist.eq.8) then
      open(1,file='table30.dat',status='old')
      c10 = 126
      elseif(moist.eq.9) then
      open(1,file='table31.dat',status='old')
      c10 = 120
      elseif(moist.eq.10) then
      open(1,file='table32.dat',status='old')
      c10 = 115
      elseif(moist.eq.11) then
      open(1,file='table33.dat',status='old')
      c10 = 156
      endif
c
c
      do i = 1,c10
      read(1,*,end=22)pe1(i),(w_ret(i,j),j=1,10)
      enddo
   22 continue
c
c  break data into yearly sets
c
      c3 = 0
      c4 = 0
      c5 = 0
      c6 = 0
      c7 = 0
      c8 = 0
c
c      print *, 'c2  = ',c2
      do i = 1,c2
      read(2,rec=i)year,mon,dummy,dummy
      if(i.eq.1) then
      read(2,rec=i+1)year_1,mon_1,dummy,dummy
      if(year.lt.year_1) then
      c3 = c3 + 1      
      c4 = c4 + 1
      c5 = c5 + 1
      b_c(c3) = 1
      e_c(c3) = 1
      elseif(year.eq.year_1) then
      c4 = c4 + 1
      c5 = c5 + 1
      endif
      elseif(i.gt.1) then
      read(2,rec=i-1)year_1,mon_1,dummy,dummy
c       print *, 'year,year1,mon,mon1  ',year,year_1,mon,mon_1
          if(i.eq.c2) then
      c3 = c3 + 1
      b_c(c3) = c5 - c4 + 1
      e_c(c3) = c2
          elseif(year.eq.year_1.and.mon.gt.mon_1) then
          c4 = c4 + 1
      c5 = c5 + 1
      else
      c3 = c3 + 1
      b_c(c3) = c5 - c4 + 1
      e_c(c3) = c5
      c4 = 1
      c5 = c5 + 1
          endif
      endif          
      enddo
      
c
c
c
      do years = 1,c3
c
      c6 = e_c(years)
      c7 = b_c(years)
      print *, 'computing evapotranspiration for year#',years,'   of',
     + c3
c
c
c      print *, 'c7  , c6   ',c7,c6
      sum = 0.0
      do i = c7,c6
      read(13,rec=i)p_pe
      sum = sum + p_pe
      enddo
c
c
c  compute the accumulative water loss potential if the total p-pe is positive
c
      if(sum.ge.0.00) then
            c9 = 0
       c99 = 0
           sum2 = 0.0
      do i = c7,c6
      read(13,rec=i)p_pe
      if(p_pe.gt.0.0) then
      c9 = c9 + 1
      else
      goto 9999
      endif
      enddo
 9999   continue
c
           do i = c7,c6
      read(13,rec=i)p_pe
           c99 = c99 + 1 
             if(p_pe.ge.0.0.and.c99.le.c9) then
              act_pot_wl(i) = 0.0
              elseif(p_pe.lt.0.0) then
              sum2 = sum2 + p_pe
              act_pot_wl(i) = sum2
         elseif(p_pe.ge.0.0.and.c99.gt.c9) then
         act_pot_wl(i) = p_pe
              endif
           enddo
c
C  FIND THE STORAGE VALUES FOR THE POS SUMS
c
            do i = c7,c6
            sum7 = abs(act_pot_wl(i))
       row = int(sum7/10.)
       remainder = sum7 - (float(row) * 10.)
       col = int(remainder)
       col = col + 1
       row = row + 1
C
            st(i) = w_ret(row,col)
            enddo
c
c
c  find change in storage
c
          do i = c7,c6
      if(i.eq.1) then
        delta_st(i) = 0.0
      else
        delta_st(i) = st(i) - st(i-1) 
      endif
          enddo
      endif
c
c
c  find the water loss if p-pe is negatice
c
c
      if(sum.lt.0.00) then
c
            sum_neg = 0.0
       sum_pos = 0.0
       do i = c7,c6
       read(13,rec=i)p_pe
       if(p_pe.lt.0.00) then
         sum_neg = sum_neg + p_pe
       else
         sum_pos = sum_pos + p_pe
       endif
       enddo
c
c  find the value of soil moisture retention equal to the total negative
c  value of p - pe from table 
c
            sum_neg = abs(sum_neg)
c
c  find the proper column and row in the mositure table
c
       row = int(sum_neg/10.)
       remainder = sum_neg - (float(row) * 10.)
       col = int(remainder)
       col = col + 1
       row = row + 1
c
            w_ret2 = w_ret(row,col)
       w_ret3(1) = w_ret2
c
c  add the positive sum to this value to obtain a value at the end of
c  the moist period
c
            do k = 2,100
            w_ret2 = w_ret2 + sum_pos
c       print *, 'sum pos  = ',sum_pos
c
c find the value potential evaporation value on the table's border
c
            pe1_2 = 0.0
            do i = 1,c10
       do j = 1,10
             if(j.eq.1) then
          if(w_ret2.eq.w_ret(i,j)) then
          pe1_2 = pe1(i)
               goto 16
               else
          endif
        else
         if(w_ret2.eq.w_ret(i,j)) then
         pe1_2 = pe1(i) + float(j-1)
         goto 16
         else
         endif
           endif
      enddo
      enddo
  16       continue
c
c  now add the total negative p-pe to this value
c
c           print *, 'pet1_2  ',pe1_2
           pe1_2 = pe1_2 + sum_neg
c           print *, 'pet1_2  ',pe1_2
c
c  now find the value of the soil moisture retention at the end of the
c  dry period in the body of the table
c
       row = int(pe1_2/10.)
       remainder = pe1_2 - (float(row) * 10.)
       col = int(remainder)
       col = col + 1
       row = row + 1
c       print *, 'row col = ',row,col
c
            w_ret2 = w_ret(row,col)
       w_ret3(k) = w_ret2
c
c  check to see if the moisture retention is the same
            if(w_ret3(k).eq.w_ret3(k-1)) then
       row  =int( (w_ret2 + sum_pos)/10.)
       remainder = (w_ret2 + sum_pos) - (float(row)*10.)
       col = int(remainder)
       row = row + 1
       col = col + 1
c            print *, 'row,col, w_ret2  ',row,col,w_ret2
       water_loss = -w_ret(row,col)
       go to 188
       else
       endif
            enddo
  188  continue
c
c
c  find the actual potential water loss
           sum2 = water_loss
           count = 0
           do i = c7,c6
           read(13,rec=i)p_pe
              if(p_pe.ge.0.0) then
        act_pot_wl(i) = 0.0
              else
         count = count + 1
         if(count.eq.1) then
         act_pot_wl(i-1) = sum2
         sum2 = sum2 + p_pe
              act_pot_wl(i) = sum2
         else
              sum2 = sum2 + p_pe
              act_pot_wl(i) = sum2
              endif
         endif
           enddo
c
c  Find storage by finding the act_pot_wl in the body of the table and
c  obtaining it from the border then obtaining the positive ones by 
c  subtracting or adding to the water loss potential
c
      start = 0
      begin = 0
      check = 0
      finish = 0
      do i = c7,c6
      read(13,rec=i)p_pe
      if(p_pe.gt.0.0) then
      else
      check = check + 1
      if(check.eq.1) then
      start = i
      else
      endif
      endif
      enddo
      finish = start + check
      begin = start - 1
c
c
c   look up storage values
c
      sum7 = 0.0
      do i = c7,c6
      if(act_pot_wl(i).lt.0.0) then
c           
            sum7 = abs(act_pot_wl(i))
       row = int(sum7/10.)
       remainder = sum7 - (float(row) * 10.)
       col = int(remainder)
       col = col + 1
       row = row + 1
            st(i) = w_ret(row,col)
c
       else
       st(I) = 0.0
       endif
       enddo
c
      st2 = st(begin)
      begin2 = begin - 1
      do i = begin2,c7,-1
      read(13,rec=i+1)p_pe
      st2 = st2 - p_pe
      st(i) = st2
      enddo
c     
      st2 = st(finish-1)
      do i = finish,c6
      read(13,rec=i)p_pe
      st2 = st2 + p_pe
      st(i) = st2
      enddo

C
C
c
c  compute change in storage which is zero if the storage is above 0
c  and is the difference from one month to the next.
c
      do i = c7,c6
      read(13,rec=i)p_pe
      if(p_pe.lt.0.0.and.i.gt.begin) then
      delta_st(i)= st(i) - st(i-1)
      else
      delta_st(i) = p_pe      
      endif
      enddo
c
      ENDIF
c
c
c
c  compute the actual evapotranspiration when precipitation is greater
c than the potential evapotranspiration, the soil remains full of water
c and the actual evapotranspiration will equal the potential.  When the
c precipitation drops below the potential evapotranspiration the soil begins 
c to dry out and actual evapotranspiration becomes less that the potentially
c  possible.  In those months the actual evapotranspiration equals the 
c precipitation plus the amount of water drawn from the soil moisture 
c storage (delta_st)
c
      do i = c7,c6
      read(2,rec=i)year,mon,dummy,precip
      read(13,rec=i)p_pe
      read(7,rec=i)pet
      if(p_pe.ge.0.0) then
      ae1 = pet
      else 
      ae1 = precip -  delta_st(i)
      endif
      c8 = c8 + 1
      ae = ae1 
      stt = delta_st(i)
      write(14,rec=c8)ae,stt
c      print *, 'c8    ',c8
      enddo
c
      ENDDO
c
c
      return
      end