c$Author: fredfox $
c$Date: 2003-02-27 19:47:39 $
c$Revision: 1.3.2.1 $
c$Source: /weru/cvs/weps/weps.src/hydro/dvolw.for,v $
      subroutine dvolw(neqn,tsec,volw,wfluxn)

c     + + + PURPOSE + + +
c     Returns the values of net flux when given a simulation time and
c     state of the system. It is used by LSODA to integrate forward in time

c     + + + KEYWORDS + + +
c     soil water redistribution, evaporation, deep percolation, runoff

c     + + + ARGUMENT DECLARATIONS + + +
      integer neqn(*)
      real tsec, volw(*), wfluxn(*)

c     + + + ARGUMENT DEFINITIONS + + +
c     neqn(1) - number of Ordinary Differential Equations to be integrated
c     neqn(2) - flag to activate internal printing
c     tsec    - time in seconds since lsoda was initialized
c     volw(*) - array of state variables
c               (1-5) - rain, runoff, evap, infil, pond (m)
c               (6-imaxlay) - volume of water in each soil layer (m)
c               (imaxlay+1) - drainage (m)
c     wfluxn(*) - net flux rates corresponding to volw(*) (m/s)

C*** FUNCTION DECLARATIONS ***
      real*4 soilrelhum,
     &       unsatcond_bc,
     &       vaporden, diffusive,
     &       airtempsin, satvappres

c*** INCLUDED COMMON BLOCKS ***
      include 'p1werm.inc'
      include 'hydro/dvolwparam.inc'
      include 'hydro/vapprop.inc'

C*** LOCAL DECLARATIONS ***
      integer     lrx, imaxlay
      real theta(layrsn), wfluxr(layrsn+1),
     &     conda(layrsn), potm(layrsn), evap,
     &     intenspeak, diffua(layrsn), fluxv(layrsn), fluxw(layrsn),
     &     airtemp, airvappres, airsatvappres, airhumid, airvapden,
     &     rainduration, tday, max_infil_rate, frac_pond_area,
     &     max_evap_rate, soil_evap_rate, pond_evap_rate,
     &     pond_infil, rain_infil, rain_pond
c     saved between calls
c     cond(layrsn), swm(layrsn),
c     soilrh(layrsn), soilvapden(layrsn), soildiffu(layrsn)

      real*4 pi
      parameter( pi = 3.1415927 )

C*** LOCAL DEFINITIONS ***
c     locally, the soil layers go from 1 to layrsn. In the passed arrays,
c     the soil layers go from 6 to layrsn+5

      if( neqn(2).gt.0 ) then
          write(*,*) 'time ', tsec
          do lrx = 1,layrsn
              write(*,*) 'lay volw', lrx, volw(lrx+5)
          end do
      end if

c     set time of the beginning of the day
      tday = tsec - beginday

c     soil state
      imaxlay = layrsn + 5
      do lrx = 1,layrsn
          if( volw(lrx+5).ne.lastvolw(lrx+5) ) then
              theta(lrx) = volw(lrx+5)/tlay(lrx)
c             Brooks and Corey soil properties
              call matricpot_bc(theta(lrx), thetar(lrx), thetas(lrx),
     &                      airentry(lrx), lambda(lrx), thetaw(lrx),
     &                      theta80rh(lrx), soiltemp(lrx),
     &                      potm(lrx), soilrh(lrx) )
              swm(lrx) = potm(lrx) - depth(lrx)
              soilvapden(lrx) = vaporden( soiltemp(lrx), soilrh(lrx) )
              soildiffu(lrx) = diffusive(theta(lrx), thetas(lrx),
     &                               soiltemp(lrx), atmpres )
              cond(lrx) = unsatcond_bc(theta(lrx),thetar(lrx),
     &                    thetas(lrx), ksat(lrx),lambda(lrx))
          end if
      end do

C     calc the conductivity and rate of flow in for each layer
      do lrx = 2,layrsn
          conda(lrx) = (cond(lrx-1)*tlay(lrx-1)+cond(lrx)*tlay(lrx))
     &               / (2*dist(lrx))
          diffua(lrx) = (soildiffu(lrx-1)*tlay(lrx-1)
     &                + soildiffu(lrx)*tlay(lrx)) / (2*dist(lrx))
          fluxw(lrx) = (swm(lrx-1)-swm(lrx)) * conda(lrx)
     &                / dist(lrx)
          fluxv(lrx) = ((soilvapden(lrx-1)-soilvapden(lrx))
     &              * diffua(lrx) / denwat) / dist(lrx)
          wfluxr(lrx) = fluxw(lrx) + fluxv(lrx)
      end do
      wfluxr(layrsn+1) = cond(layrsn)

c     boundary fluxes
c     rainfall intensity = wfluxn(1)
c     runoff rate = wfluxn(2)
c     evaporation rate = wfluxn(3)
c     infiltration rate = wfluxn(4)
c     change in ponded depth = wfluxn(5)
c     drainage rate = wfluxn(imaxlay+1)

c     rainfall intensity
      if(     (tday.ge.rainstart)
     &   .and.(tday.le.rainend)
     &   .and.(raindepth.gt.0.0) ) then
          rainduration = rainend - rainstart
          intenspeak = 2.0*raindepth/rainduration
          if( tday.lt.rainmid ) then
              wfluxn(1) = intenspeak
     &                  * (tday-rainstart)/(rainmid-rainstart)
          else
              wfluxn(1) = intenspeak
     &                         * (tday-rainend)/(rainmid-rainend)
          end if
      else
          wfluxn(1) = 0.0
      end if

      pondmax = max(0.001, pondmax) !avoid div by zero
c     generate runoff if above retention limit (kind of like WEPP)
      if(volw(5).ge.pondmax) then
          wfluxn(2) = ((volw(5)-pondmax)**1.5)
     &              * (8.0*gravconst*soilslope/dw_friction)**0.5
     &              /  slopelength
      else
          wfluxn(2) = 0.0   ! no runoff
      end if

c     determine fraction of soil area covered by ponding
      if( volw(5).ge.0.0 ) then
          frac_pond_area = min(1.0,(volw(5)/pondmax)**0.5)
      else
          frac_pond_area = 0.0
      end if

c     calculate maximum evaporation rate
c     this method using a lenday that is 12 hours (43200 sec) allows
c     running time continuously through many days
c      wfluxn(3) = max(0.0,evapamp * sin(pi*(tday-sunrise)/lenday))
c     this method assumes that tday and sunrise and sunset are in the
c     same day and lenday can be actual daylight hours
      if( tday.gt.sunrise.and.tday.lt.sunset ) then
          max_evap_rate = evapamp * sin(pi*(tday-sunrise)/lenday)
      else
          max_evap_rate = 0.0
      endif

c     find air relative humidity from diurnal air temperature
c     and dewpoint temperature
      if( tday .lt. 21600 ) then
          airtemp = airtempsin( tday, airtmaxprev, airtmin )
      else if( tday .lt. 64800 ) then
          airtemp = airtempsin( tday, airtmax, airtmin )
      else
          airtemp = airtempsin( tday, airtmax, airtminnext )
      end if

      airvappres = satvappres(tdew)
      airsatvappres = satvappres(airtemp)
      airhumid = airvappres/airsatvappres
      airvapden = vaporden( airtemp, airhumid )

c     calculate reduction in evaporation rate due to dry soil
      if( (soilrh(1).gt.airhumid).and.(airhumid.le.0.99) ) then
          soil_evap_rate = max_evap_rate
     &                   * (soilrh(1)-airhumid)/(1.0-airhumid)
      else
          soil_evap_rate = 0.0
      end if

c*** try this
c      soil_evap_rate = min( max_evap_rate,
c     &                  (soilvapden(1)-airvapden)
c     &                 * 10. * soildiffu(1) / denwat) / dist(1))
c
c      if( neqn(2).gt.0 ) then
c          write(*,*)'dvolw:tcddf',
c     &    tday, airtemp, soilvapden(1), airvapden, wfluxn(3)
c      endif
c      if( neqn(2).gt.0 ) then
c          if(max_evap_rate.gt.0.0)
c     &    write(*,*)'dvolw:evap',soilrh(1),airhumid, max_evap_rate
c      endif
c*** end try this

c     split evaporation between soil and pond
c     fluxv(1) is negative since this is a loss to the soil layer
c     consequently it mus be subtracted to make it an addition to evap.
      fluxv(1) = -soil_evap_rate * (1.0-frac_pond_area)
      pond_evap_rate = max_evap_rate * frac_pond_area
      wfluxn(3) = -fluxv(1) + pond_evap_rate

c     calculate max infiltration rate !note, the potential difference
c     should be (0.0-swm(1)), but this only works if we account for
c     hysteresis. As it is written, infiltration will not exceed
c     saturation. When hysteresis is incorporated, sorbing soil will
c     have an airentry potential of 0.0 again making them match.
      max_infil_rate = (airentry(1) - swm(1))
     &               * 0.5 * (ksat(1)+cond(1))/dist(1)

c     add infiltration from both ponded area and rainfall area
      pond_infil = max_infil_rate * frac_pond_area
      rain_infil = min(wfluxn(1),max_infil_rate) * (1.0-frac_pond_area)
      wfluxn(4) = pond_infil + rain_infil
      fluxw(1) = wfluxn(4)

c     rate rain is added to the pond
      rain_pond = wfluxn(1) - rain_infil

c     rate change in pond depth (rain - infil - evap - runoff)
      wfluxn(5) = rain_pond - pond_infil - pond_evap_rate - wfluxn(2)

c     resultant flux at surface layer
      wfluxr(1) = fluxw(1) + fluxv(1)

c     drainage
      wfluxn(imaxlay+1) = wfluxr(layrsn+1)

c     calc the net flow into (+) or out of (-) each layer
      do lrx = 1, layrsn
          wfluxn(lrx+5) = wfluxr(lrx) - wfluxr(lrx+1)
          if( neqn(2).gt.0 ) then
              write(*,*) 'lay ', lrx, wfluxr(lrx), wfluxn(lrx+5)
          end if
          lastvolw(lrx+5) = volw(lrx+5)
      end do

c      if( neqn(2).gt.0 ) then
c          if( (wfluxr(1).lt.0.0).or.(wfluxr(2).gt.0.0) )
c     &    write(*,*) 'flx 1:',wfluxr(1), wfluxr(2)
c      endif

      return
      end

c     dummy jacobian matrix
      subroutine jac (neq, t, y, ml, mu, pd, nrowpd)
      integer neq, ml, mu, nrowpd
      real*4 t, y, pd
      dimension y(*), pd(nrowpd,*)
c     the full matrix case (jt = 1), ml and mu are ignored,
c     and the jacobian is to be loaded into pd in columnwise manner,
c     with df(i)/dy(j) loaded into pd(i,j).
c     in this case, dwfluxn(i)/dvolw(j) in pd(i,j)
      return
      end

c table lookup function archived here for reference
c      real*4 function afgen(numpts, table, xval)
c     Finds two entries in table( ,1) which bracket the value of the independent
c     variable, xval, and then linearly interpolates between the two corresponding
c     entries in Table( ,2) to find the corresponding value of the dependent
c     variable.
c
c      integer   numpts
c      real*4    table(numpts,2), xval
c
c      integer maxindex, minindex, midindex
c
c      maxindex = numpts
c      minindex = 0  !Place limit on one end of output.
c	  Set limits on values the function may take.
c	  if( xval.lt.table(1,1)) then
c          afgen = table(1,2)
c      else if( xval.gt.table(numpts,1)) then
c          afgen = table(numpts,2)
c      else
c          do while((maxindex-minindex).gt.1)
c              midindex = (maxindex+minindex)/2
c              if((table(numpts,1).gt.table(1,1))
c     &            .eqv.(xval.gt.table(midindex,1))) then
c                  minindex = midindex
c              else
c                  maxindex = midindex
c              end if
c          end do
c          afgen = (xval-table(minindex,1))
c     &          * (table(minindex+1,2)-table(minindex,2))
c     &          / (table(minindex+1,1)-table(minindex,1))
c     &          + table(minindex,2)
c          write(*,*) 'afgen int:', afgen
c      end if
c
c      return
c      end
c