Ground Water Recharge by Vetiver Hedgerow

; Vetiver hedgerow Model
; Version 1.0.0 Date 28 Aug 2011  
; created for netlogo 4.1.3 by Vinod Kumar, C-DAC Mumbai, India
; transitioned to netlogo 5.0.5 of 19 Dec 2013
; 
; vim settings
; set ts=2
; set ai


; versions
; avetqtoh0d4.nlogo.14sep2011  model with surface runoff reasonably ok.
;              now try diffusion and then preferential flow  from Nimmo
; avetqtoh0d4.nlogo.sep182215
;     routines to relate water content (theta), suction (psi) and Hydraulic 
;     conductivity (KTheta) for walla walla
; avet.nlogo.20sep1040 colors for patches cleaned up; surface patch blackening
;              removed.  routine for d theta/ d psi added. Diffusion
;              not added
;avet.nlogo.26sep1730 with moisture transport by darcy's modulated by
;              available water content and moisture equalization
;              fingers, patches of low and high moisture occur
;
;avet.nlogo.27sep1231 attempt to set constant qi by computing new theta from
;              old values of other parameters and the new qi. The idea
;              is to adapt the steady state flux equation to
;              slowly varying inputs. But in the running it seems to
;              generate bands and patches of high saturation. screenshot
;              seen in avet.nlogo.27sep1231.jpg
;              the solution is highly unstable. when the theta
; exceeds some value, the outflow computed is so high that the theta
; becomes negative and common sense brings it back to residual moisture.
; thus the scheme working with psi values is unstable. so work with
; theta and D(theta).
;
; avet.nlogo.29sep2400 darcy-buckingham with qi equated to 
;   qo + moisture increase. plot of moisture vs depth. 
;   in 60 mins the moisture diffused around 12 cms. According to
;   Ks (345.6 mm/day =  1.44 cm/hour) this is about 10 times.
;   Probably the Diffusion and matric suction are the cause.

; avet.nlogo.4oct2300 
;   root is actually a cylinder with thin surface film flow.
;   assume RLD =2e5 m^-2, R = .33 mm uniform film thickness and measured
;   max velocity =0.15 e-3 m/s. film thickness hu = 5.516 e-6 m
;   and qcyl = 0.2306 e-6 m^3/s

;   If the cylindir surface is assumed flat
;   fad = RLD*2pi*R = 414.678 m^-2. assume standard vulu 5.51599 e-10 m^2/s
;   qflat = 0.2294 e-6 m^3/s
;   Use flatttened surface model instead of cylindrical surface
;
;   Assume the water table to be far below and our volume of interest
;   to be sufficiently above it. So we assume that any inflow
;   downward to the last layer will move downward into the nether
;   regions as is whereis. This is similar to the assumption made
;   on the surface where the outflow from the last patch on the right
;   is assumed to vanish into the unknown.
;   The water crossing the last layer at bottom is to be monitored
;   and indicates how effective vetiver hedgerows and their roots
;   affect the water recharge
;   This simplifies the flow processes. We can move from input to output
;   for diffusion and surface flow without worry about any backlash.
;   However for the preferential flow, we still have to consider
;   saturation at the bottom of the root zone and subsequent backlash.
;   as a first step, assume root zone to extend to lowermost patch
;   so that backlash is avoided even for preferential flow.
; 
; avet.nlogo.4oct1700
;   preferential flow working. But no graphs. no coupling between
;   two domains
; avet.nlogo.5oct20111800
;   preferential flow metered and graphed. preferential flow 
;   about 13.8 times more intense than diffuse flow
;   After 70 mins, 63.6 cm^2 of water diffused over 40 cm wide area.
;   Height 1.59 cm in 70 mins = 32 cm/day= 0.022714 cm/min
;   and 44 cm^2 of  water filmed over 2.4 cm of vetiver width.
;   Height of water through film = 18.33 cms in 70 mins 
;   = 0.2619 cm/min film flow/diffuse flow = 13.8365
;   film flow rate = 2.76 mm/min leads to 19.32 cm in 70 mins. 
;   (ok 18.33 cms mesured) 
;   diffuse flow rate = 0.23 mm/min = 1.61 cm in 70 mins 
;   ok. Measured 1.59 cms after 70 mins diffuse 63.6 cm^2, film 44 cm^2
; avet.nlogo.6oct20111815
;   serious doubts about surface flow. this version is with old scheme.
;   the problems listed in modelGWRabstract.odt. now modify in light of
;   moisture and film flow experience
; avet.nlogo.7oct20111035
;   zeroRef drain introduced to simplify diffuse and preferential flow
;   and to take care of surface water level, especially behind the hedgerow
;   surface flow corrected but still not certain. surface water height
;   taken into consideration but no observable effect on diffused
;   water but effect on KTheta observable.
;   vetiver hedgerow density can now be changed and its effect on
;   water height at vetiver and the backwater is visible.
; avet.nlogo.7oct20111700
;   be convinced about the correct handling of qi, height h for surface flow 
; avet.8oct20111000.nlogo
;   avet before adding meters for water balance, removal of rain
;   over area of interest (rain only in hinterland) and 
;   reduction of hinterland to area of interest to study
;   water balance when the rain roughly matches the infiltration
; avet.8oct20111600.nlogo
;   first remove rain from area of interest. rain only in hinterland
; avet9oct0040.nlogo
;   quantities like totalrain, surface runout expressed as cm^2
;   still rain cannot balance with s+so+d+f which should have 
; mismatch and loss due to 1)initial transients and transients due to
; variation in rain. Reduced variation from 10% to 2%
;
; avet9oct1530.nlogo
;   Rationalized graphs for flows and moisture contents.
;   Enabled cleaning the slate for continueing simulation after
;   transients have died down. For this, release go, press cleanSlate.
;   Then the simulation restarts with meters reset, but moisture contents
;   and levels are steady. The world can then be saved and reused for
;   other simulations .
; avet10oct0030.nlogo
;   water height at vetiver for low heights taken care off.
;   area reduction included in Rv calculation. 
; avet.10oct0030.hinterland1m.nlogo
;   roughly match input to diffusive and film flow so that their
;   effects can be compared. Otherwise almost 100 % of the input
;   just ran out at the surface
;
;   avet.10oct0030.hinterland1m.vet100mins140rain1cm.world.csv
;     problems. total water flowed surfaceout 90.2 % Diffusein 54 %
;     filmin 11.7 %. Cannot be. Check whether diffusein and filmin
;     are overstated 10 times
;     1 cm/hr over 1 m will lead to 1.6667 cm^2 per min
;     or 233.33 cm^2 in 140 mins
;     but diffused is 131 cm^2 filmed is 29 cm^2 and surface out 220.82 cm^2
;     totalrain 243.39 cm^2. (about 5% above. why?)
;     inflow 0.06958 (max 0.0701, min 0.06875) ok. at inlet
;     outflow at right (0.068699) 
;     Check diffused. at patch 75,58 h =9.202e-5 m = 0.009201 cm
;     Kseffective Ks = 0.0004 cm/s. 
;   qit=ks
;   diffu=qit*4*.4  per patch [cm/s * 4 s* .4 cm = cm^2] per patch
;   diffu
;   .00064000000000000000 cm^2 in 4 secs per patch
;   diffu*100*140*60/4  cm^2/patch * 100 pactches * 140 * 60 secs/ 4 secs
;   134.400 cm^2 in 140 mins over the full field (compares well with 131 cm^2)
;   check whether diffusion is not subtracted from input before output is set
;   Problem is that the surface flow is Q driven and hence the diffusion
;   should subtract from the inflow not h
;
; avet.11oct0800.nlogo
;   many water balances resolved. interface rationalized. hinterland slider.
;   many csvs   based on this. 
;
; avet.16oct2400.nlogo
;   version with uncontrollable vetiver root length density Lv. most of the
;   csvs are built using this version. Later versions have switch
;   for  varying Lv and procedure called from Go to scan the switch and set
;   variables. Such dynamic change cannot happen for Slope as the dimensions
;   and patches change with it. The update can also handle vetiver hedgerow
;   geometry change N_vs within the update routine instead of within
;   the surface runoff as it is now.
; avet.nlogo
;   for vetiver hedgerow Rv calculation there was a constant D that was
;   set only in initailly (setup). Now do it in lookAtSwitches D that was
;   set only in initailly (setup). Now do it in lookAtSwitches. 
;   May have to rerun many csvs?
; --------------------------------------------------------------
; sizing of the simulation
; assume newtonian mechanics of frictionless balls falling
; on inclined plane of slope theta. from energy considerations
; vn = sqrt (g l tan(theta))  where l is the horizontal distance
; find velocity at 3/4 l where l is the horizontal length.
; and assume that the ball should still be within l with that velocity
; this gives l = 12 delT**2 g tan(theta)
; for tan theta = 1/3  and deltaT = 1 s, l = 39 m
; for tan theta = 1   and deltaT =1, l = 120 m
; so assume 120 m, deltaT = 1.0
;  then tan theta 6/120 and l becomes 5.88 theta around 3 degree
; so assume l = 60 m h = 6 m and max 1/10 slope
; 
globals [ 
  ; for preferential flow as thin film on vetiver root surface
  ; assume location independent 
  ;Lv        ; root length density m^-2. depends on location
  rootR     ; root radius m
  Vu        ; max avg film flow velocity m/s
  VuMax     ; max film flow velocity m/s
  HuC       ; film thickness on root cylinder m
  rhc       ; radius of root and fill m
  HuF       ; film thickness on flat surface m = Lu
  VuHuC     ; of cyl thin film flow m^2/s(water
  VuHuF     ; of flat thin film flow m^2/s(water)
  fad       ; facial area density m^-1

  thetaRoots        ; volume of roots in 1 m^3 [m^3 m^-3]
  thetaFilmMax      ; max volume of film water +root per unit vol (m^3 m^-3) 
  thetaSatWithRoots ; thetas - thetaFilmMax [m^3 m^-3]
  qFilmMax          ; m s-1

  ; Field Cross section dimensions in m
  ;hinterland ; m
  cswidth csheight ; width and height of cross section in m
  pDimen ; patch measures pDimen m X pDimen m
  pDcm ; pDimen inn cm
  horizonL  ; land horizon at left 
  ;slpAngle ; slope down to right in degrees
  mSlope ; tan (slpAngle)
  v1 ; velocity from patch to patch+1 m/s
  ; horizonY = horizonL - mSlope*X
  noVetiverRootsHere ; no vetiver roots here at the bottom of cross section
  zeroRef           ; top of drain. Pressure here is air pressure.
  drainPatches ;
  ;rainfallIntensity ; cm/hour

  ; colors of various patches
  vetiverRootPColor   
  diffusiveSoilPColor   
  vetiverShootPColor    
  surfacePColor     
  infoPColor
  aakaashPColor
  topPColor
  drainColor

  rainFallsHerePS ; agentSet of patches.  when 'rain drops keep falling
                  ;  on my head, therainFallsHere is head.
  patchSetWithDiffusion ; agentset of patches where diffusion takes place
                        ; includes patches under vetiver
  patchSetWithFilmFlow ; agentset of patches where Film flow takes place
  surfPatchInfoRow ; the 0 row contains info about the rainFallsHerePS
  ;tickTime ; s
  ; simulateDuration; mins
  simulatedTime; mins
  ; surfaceRun ; switch on surfaceRun
  ;     mannings formula 
  ; avg velocity v = (kn / n) R ^ (2/3) S ^ (1/2)
  kn ; 1.0 in SI units
  n  ; manning coeff
  g  ; acceleration gravity m/s**2
  D  ; stands for Kn/n sqrt(mSlope) (1 - 2 * vetiverpm * rootR)

  ; vars for backwater for vetiver hedgerow
  ; vetiver hedge is modeled as a vertical plane on line at vetiverXPos
  ; the hedgeWidth comes before the line
  hedgeWidth ; m
  vetiverXPos ; m vetiver hedge here 
  vetLpxcor ;
  vetRpxcor ;
  deadHeight ; m th econceptual dead water height at the vetiver hedge
             ; = Hv - H. (Hv height with modified R due to vetiver
             ;            H height at vetiver with old R.
             ; With vetiver
             ; R at vetiver hedge = 1/2i m .
             ; i  vetiver tillers/m
             ; without vetiver
             ; R  =H the height of water without vetiver
  backwaterLength ; m = deadHeight/sin(slpAngle) . 
  backwaterStartPos ; m 
  ; within the backwaterLength  at x from the start f backwater
  ; height = height from manning equation + dead water height
  ;        = meadHeight*x/backwaterLength
  ; 
  ; vetiverpm ; slider vetiver tillers per m length of hedgerow
  hv ; height of water at hedgerow m
  Rv ; manning radius for hedge 1/(2 vetiverpm)    metres
  vetiverShootRadius ; m

  ; vars for soil parameters
  thetas    ; saturated water content m^3 m^-3 = porosity
  thetaResidual ; residual moisture
  psid      ; owen dry suction; head units; cm of water
  Ks        ; saturated hydraulic conductivity cm/s
  c         ; parameter of Rossi-Nimmo 3 parameter sum model
  alpha     ; parameter of Rossi-Nimmo 3 parameter sum model
  psii      ; cm parameter of Rossi-Nimmo 3 parameter sum model
  psio      ; cm  parameter of Rossi-Nimmo 3 parameter sum model
  lamda     ; parameter of Rossi-Nimmo 3 parameter sum model
  powerLaw  ; component of theta from power law  in th emid region
  logLaw    ; component of theta from log law in the low theta region
  parabolaLaw   ; component of theta from parabola law in the saturated
  
  ; end vars for soil parameters

  ; vars for plots
  totalRain ; cm^2 monitored at field input
  totalSurfaceWater ; cm ^2
  soilMoisture      ; cm^2
  filmMoisture      ; cm^2
  surfaceRunout ; cm^2 
  diffusedWater ; cm^2
  filmedWater ; cm^2
  filmRunout ; cm^s
  waterLevelBeforeVetiver ; cm
  waterLevelAfterVetiver ; cm
  waterLevelAtVetiver ; cm
  leftPatchForPlot ; 
  rightPatchForPlot ; 
  plotAt; tick at which to plot
  bigBang ; the time at which the world starts
  ; simulatedTime; mins
]

patches-own
[
  ; for preferential flow as thin film on vetiver root surface
  qipref      ; inflow head units cm/s
  qopref      ; outflow head units cm/s
  thetapref   ; volumetric water content
  aaf         ; active area fraction

  surfacepYcor ; the pycor of the surface stored in row 0
  qi          ; input flow to patch in cm of water at begin of tick
  qo          ; output flow from the patch at the end of the tick. input
              ; to next patch at the next instant
  ; diffuse flow parameters for below surface patches
  thetaDifus  ; volumetric water content m^3/m^3
  theta1Dash  ;
  thetaRatio  ; theta/thetas
  psiPatch    ; head units cm of water. = -suction. matric potential
  transWater  ; water in head units cm to be moved down
  z           ; depth cm below surface level
  KTheta      ; unsat hydraulic conductivity in head units cm/s
  KThetaFactor; to take care of water level above ground
  dThetaBydPsi ; in head units 1/cm
  DTheta ; diffusivity in head units cm^2/s
  ; end diffuse flow parameters for below surface patches

  ; surface flow parameters
  waterInPatch ; water cm  in each patch. filled and water gets higher. 
  hm ;  m normal manning height of water depends only on slope and roughness
  h ; surface water height in m. 
  hprev ; surface water in m at prev tick
  v ; velocity of water col m/s
  runoff ; in terms of cm per tickTime
  ;surfacepYcor ; the pycor of the surface stored in row 0
  originalColor ; setup initializes
  inFlow ; cm/s  outFlow of up patch in previous instant. For first flow
          ; from hinterland
  outFlow ; cm/s outFlow = inFlow - change in storage
  ; end surface flow parameters


]

; routines for soil states
; to-report KThetaRatioMethod2 [ thetah ]
; to-report thetaRatioToPsi [thetaRatio]
; to-report psiToThetaRatio [ psi ]
; to-report dThetaRatioBydPsi [ psi ]
; to-report DThetaReporter [ thetah ]

to setup
  ;; (for this model to work with NetLogo's new plotting features,
  ;; __clear-all-and-reset-ticks should be replaced with clear-all at
  ;; the beginning of your setup procedure and reset-ticks at the end
  ;; of the procedure.)
  __clear-all-and-reset-ticks
  
  init-soil-params;
  ; colors of various patches
  set vetiverRootPColor   (yellow + 1) 
  set diffusiveSoilPColor   (brown + 1)   
  set surfacePColor     (brown + 2)
  set vetiverShootPColor    (green - 2)
  set infoPColor        (magenta)
  set aakaashPColor     (gray)
  set topPColor       (gray - 1)
  set drainColor      sky

  ; ground cross section paramaters
  set hinterland 1 ; m
  set cswidth .4 ;m
  ;set csheight 12 ; m
  ;set pDimen 0.5 ; m
  set vetiverXPos (cswidth * 3 / 4) ; m
  set backwaterStartPos vetiverXPos
  set noVetiverRootsHere .04 ; m
  set zeroRef noVetiverRootsHere  ;m   Pressure here is air pressure.
  ; set slpAngle 30 ; slider slope down to right in degrees
  set mSlope  tan slpAngle
  set hedgeWidth 0.02 ; m
  ;set max-pxcor cswidth / pDimen
  ;set max-pycor csheight / pDimen
  set pDimen cswidth / (max-pxcor + 1 )
  set pDcm  pDimen * 100 ; cm
  set csheight pDimen * (max-pycor + 1)
  set horizonL csheight - 0.04 ; m
  ;set rainfallIntensity 20 ; slider cm/hour
  ; set tickTime 1 ; sec slider
  ;set simulateDuration 20 ; mins
  set surfaceRun true ; switch on surfaceRun
  ; for checking water balance the quantities will be measured as cm*cm
  ; totalRain is the only source of water and it will be shown as 100 %
  ; all other water contents will be expressed and shown as % of totalRain
  set totalRain 10 ; cm^2 should be ideally 0 but then at time 0 div by 0
  set totalSurfaceWater 0 ; cm ^2
  set soilMoisture 0 ; cm ^2
  set filmMoisture 0 ; cm ^2
  set surfaceRunout 0 
  set filmRunout 0 
  set v1 pDimen / tickTime ; m/s
  set diffusedWater 0.0;
  set filmedWater 0.0;
  
  ; for preferential flow as thin film on vetiver root surface
  ;set Lv 0 ; 2.0e5                    ; root length density m^-2
  set rootR 0.33e-3             ; root radius m
  set Vu 1.0e-4                 ; max avg film flow velocity m/s
  set VuMax 1.5e-4              ; max film veleocity m/s
  set HuC  100e-6 ; 5.51599e-6  ; film thickness on root cylinder m
  set rhc rootR + HuC           ;
  set HuF  5.530957e-6          ; film thickness on flat surface m = Lu
  set VuHuC 150e-10   ; 5.51599e-10  ; of cyl thin film flow m^2/s(water) 
  set VuHuF 5.530957e-10        ; of flat thin film flow m^2/s(water) 
  set fad (Lv * 2 * pi * rootR) ; facial area density m^-1

  set thetaRoots        (Lv * pi * rootR * rootR)
  set thetaFilmMax      (Lv * pi * rhc * rhc)  ;(fad * HuF)
  set thetaSatWithRoots (thetas - thetaFilmMax)
  set qFilmMax          (thetaFilmMax - thetaRoots) * VuMax  ;Lv * pi * (rhc * rhc - rootR * rootR) * VuMax;


  ; avg velocity v = (kn / n) R ^ (2/3) S ^ (1/2)
  ; from The Engineering ToolBox www.EngineeringToolBox.com
  ; http://www.engineeringtoolbox.com/mannings-formula-gravity-flow-d_800.html
  ;  Manning's Formula for Gravity Flow
  ;  Manning's equation for calculating gravity flow in open channels 
  ; http://www.engineeringtoolbox.com/mannings-roughness-d_799.html
  ;   Manning's Roughness Coefficient
  ;   Manning's roughness coefficient values for some common materials 
  set kn 1 ; 1.0 in SI units
  set n  0.035 ; manning coeff for flood plains (pasteur, farmland)
               ; or Earth channel (stony, cobbles)
  set g 9.81 ; acceleration gravity m/s**2
  set vetiverShootRadius (1.2 / 1000) ;1.2 mm radius for each shoot
  ;set Rv  (0.5 / vetiverpm - vetiverShootRadius)   ; manning radius for hedge
  set D  (Kn / n * (sqrt mSlope) * (1 - 2 * vetiverpm * rootR))
  
  set-default-shape turtles "square"
  
  set vetRpxcor ( floor (vetiverXPos / pDimen))
  set vetLpxcor ( floor ((vetiverXPos - hedgeWidth) / pDimen))
  set hedgeWidth ((vetRpxcor - vetLpxcor + 1) * pDimen )
  set leftPatchForPlot vetLpxcor - 40  
  set rightPatchForPlot vetRpxcor + 10 
  ask patches [
    let horY (horizonL - mSlope * pxcor * pDimen )
    ; Y = horizonL - mSlope*X
    let groundSurfaceYcor floor ( horY / pDimen )

    if (pycor <  groundSurfaceYcor) and 
        (pycor > (noVetiverRootsHere / pDimen) - 1) [
      ifelse ( pxcor >= vetLpxcor) and ( pxcor <= vetRpxcor) and
          (pycor > (noVetiverRootsHere / pDimen) - 1) [
        set pcolor vetiverRootPColor 
        set aaf 1.0  ; 0.8                   ; active area fraction
      ][
        set pcolor diffusiveSoilPColor
      ]   
    ]
    if (pycor = (noVetiverRootsHere / pDimen) - 1)
        or (pycor = (noVetiverRootsHere / pDimen) - 2) [
      set pcolor drainColor
    ]
    if pycor >  groundSurfaceYcor [
      ifelse (( pxcor >= vetLpxcor) and ( pxcor <= vetRpxcor)) [
        set pcolor vetiverShootPColor
      ][
        set pcolor aakaashPColor
      ]
    ]

    if (pycor = groundSurfaceYcor) [
      set pcolor surfacePColor
      ask patch pxcor 0 [ set surfacepYcor groundSurfaceYcor]
    ]
    set originalColor pcolor 
  ]
  set drainPatches patches with [pcolor = drainColor]
  set rainFallsHerePS patches with [pcolor = surfacePColor]
  set surfPatchInfoRow patches with [pycor = 0]
  ask surfPatchInfoRow [ 
    set pcolor infoPColor
    set originalColor pcolor]
  set patchSetWithDiffusion (patches with 
      [(pcolor = diffusiveSoilPColor) or (pcolor = vetiverRootPColor)])
  set patchSetWithFilmFlow (patches with [pcolor = vetiverRootPColor])
  ask patches with [pycor = max-pycor]
    [set pcolor topPColor]
  ask patchSetWithDiffusion [
    set z (pycor + 0.5) * pDcm
    set thetaDifus thetaResidual
    set qi 0.0
    set qo 0.0
    setOtherDiffusiveParams self
  ]
  my-setup-plots
  my-update-plots
  set plotAt ((ticks - bigBang) * tickTime) + 60 ; plot at 60 s intervals

  show word "thetaDifus 0.35 k(theta) should be  59.7588 mm/day. Is  " 
      (345.6 * (KThetaRatioMethod2  0.35)) 
  show word "psi 104.6  cm; thetaDifus should be  0.380571 Is  " 
      (thetas * (psiToThetaRatio  104.6)) 
  show word "psi 836 cm; thetaDifus should be  0.202975 Is  " 
      (thetas * (psiToThetaRatio  836)) 
  show word "psi 836 cm; d thetaDifus / d psi should be -1.88545e-05 . Is  " 
      (thetas * (dThetaRatioBydPsi  836)) 
  show word "thetaDifus .3806  psi should be  104.439 cm. Is  " 
       (thetaRatioToPsi  (0.3806 / thetas)) 
  show word "thetaDifus .3806  D(theta) should be  1.02741 cm^2/s. Is  " 
       (DThetaReporter  0.3806) 
  show word "thetaDifus .22  psi should be  609.407 cm. Is  " 
       (thetaRatioToPsi  (0.22 / thetas)) 

  let vetpmsaved vetiverpm
  let rootrsaved rootR
  let dsaved D
  set vetiverpm 100
  set rootR 3.3e-4
  set D 7.893229367244781
  show word word "water at 100/m vetiver hedge radius 0.33 mm Q 0.00027616 m^3/s"
    " should be 0.0025207m. Is " heightFromQForVetiverHedge (0.000276160)
  set D dsaved
  set rootR rootrsaved
  set vetiverpm vetpmsaved
end 

to init-soil-params
  ; for walla walla
  set thetas 0.39
  set thetaResidual 0.04 ; 0.07
  set psid 1e7  ; cm of water
  ;set Ks (345.6 / 1000 / (24 * 60 * 60)) ; from mm d ^-1 to m s^-1
  set Ks (345.6 / 10 / (24 * 60 * 60)) ; from mm d ^-1 to cm s^-1
  set psii 104.6    ; fitted parameters
  set psio 35.6
  set lamda 0.61
  set alpha 0.0399487   ; computed parameter
  set c  0.00280054
end 

to setOtherDiffusiveParams [aPatch]
  ; thetaDifus is known. set others
  ask aPatch [
    ; mainly parameters for eqn [3] nimmo2010
    if (thetaDifus < thetaResidual) [ set thetaDifus thetaResidual]
    if (thetaDifus > thetas) [set thetaDifus thetas]
    set thetaRatio (thetaDifus / thetas)
    set psiPatch ( - (thetaRatioToPsi thetaRatio))
    set KTheta (KThetaRatioMethod2 thetaDifus) * Ks
    set DTheta (DThetaReporter thetaDifus)
  ]
end 

; Q m^3/s
; reports height m

to-report heightFromQForVetiverHedge [ Q ]
  let h1 0.0
  let h2 1.0
  if (Q = 0.0) [ report 0.0]
  set Rv (0.5 / vetiverpm )
  let howclose (abs (h1 / h2))
  while [(howclose < 0.999) or (howclose > 1.001)] [
    set h2 (Q / (D * exp (0.6667 * ln(Rv))))
    set h1 (0.5 * (h2 + h1))
    set Rv ((h1 * (1 - (2 * rootR * vetiverpm)))  ; 
          / (1 + (2 * vetiverpm * (h1 - rootR))))
    set howclose (abs (h1 / h2))
  ]
  report h1
end 

to-report KThetaRatioMethod2 [ thetah ]
  ; Measured and modeled unsaturated hydraulic conductivity
  ; C. Chen and WA Payne
  ; method2 eqn 12
  ; returns relative conductivity KTheta/Ks
  let sx 0
  let sx175 0
  let sx1b175 0
  let ftemp 0
  let k 0
  
  set sx  ((thetah - thetaResidual) / (thetas - thetaResidual) )
  if (sx != 0.0) [
    set sx175  exp(1.75 / (1.75 - 1.0) * ln(sx))]
  if ((1.0 - sx175) != 0.0) [
    set sx1b175  exp((1.0 - 1.0 / 1.75) * ln(1.0 - sx175))]
  set ftemp  (1.0 - sx1b175)
  if (sx != 0.0) [
    set k  exp(0.555 * ln(sx)) * ftemp * ftemp]
  report k
end 

; psi in cm of water (head units)

to-report thetaRatioToPsi [thetaRatParam]
    set powerLaw 0.0
    let psipower 0.0
    let thetaiRatio 0.0
    let psi 0.0
    let psitrial 0.0
    let thetaRatParamtrial 0.0
    set thetaiRatio (psiToThetaRatio psii)
    ifelse (thetaRatParam >= thetaiRatio)[
        set psi (psio * sqrt((1.0 - thetaRatParam) / c))
    ][
        set psitrial psii
        set thetaRatParamtrial thetaiRatio
        while [thetaRatParam < thetaRatParamtrial] [
            set psitrial (psitrial * 2)
            set thetaRatParamtrial (psiToThetaRatio psitrial);
        ]
        set psi psitrial
        while [(abs (thetaRatParam - thetaRatParamtrial)) > 0.00001] [
            set powerLaw (powerLaw * thetaRatParam / thetaRatParamtrial)
            set psipower (psio / (pow (powerLaw + (pow (psio / psid) lamda))
                                  (1.0 / lamda)))
            set psitrial psipower
            set psi psipower
            set thetaRatParamtrial (psiToThetaRatio psitrial);
        ]
    ]
    report psi
end 

to-report pow [x y]
    let ret 0.0
    if (x != 0.0) [
        set ret (exp(y * ln (x)))
    ]
    report ret
end 

to-report psiToThetaRatio [ psi ]
    let ftemp 0.0
    let thetar 0.0
    set powerLaw 0.0
    set logLaw 0.0
    set parabolaLaw 0.0

    ifelse (psi <= psii) [
        set ftemp (psi / psio)
        set parabolaLaw ( 1.0 -  (c * ftemp * ftemp) )
        set thetar   parabolaLaw
    ][
        set powerLaw ( (pow (psio / psi) lamda) - (pow (psio / psid) lamda))
        set logLaw ( alpha * ln (psid / psi ))
        set thetar (powerLaw + logLaw)
    ]
    report thetar
end 

to-report dThetaRatioBydPsi [ psi ]
    let ftemp 0.0

    ifelse (psi <= psii) [
        set ftemp (-2.0 * c * psi / (psio * psio))
    ][
        set ftemp (- lamda * (pow psio lamda) / (pow psi (3.0 - lamda))
               - (alpha / psi))
    ]
    report ftemp
end 

to-report DThetaReporter [thetah]
  let dtemp 0.0
  let psitemp 0.0
  let dthetabydpsitemp 0.0
  let kthetatemp 0.0

  if (thetah >= thetas) [set thetah (thetas - 0.0001)]
  set kthetatemp((KThetaRatioMethod2 thetah) * Ks)
  set psitemp (thetaRatioToPsi (thetah / thetas))
  set dthetabydpsitemp (thetas * (dThetaRatioBydPsi psitemp))

  if (dthetabydpsitemp != 0.0) [
    set dtemp (kthetatemp / (- dthetabydpsitemp))
  ]
  report dtemp
end 

to go
  lookAtSwitches
  setInputFromOutputForSurfaceFlow
  setInputFromOutputForDiffuseFlow
  setInputFromOutputForFilmFlow
  surfaceRunoff
  diffuseMoisture
  surfaceFilmFlowOnRoots
  tick
  showSurfaceWaterLevel
  showMoistureLevel
  set simulatedTime simulatedTime + (tickTime / 60)
  my-update-plots
  if (simulatedTime >  simulateDuration  )
  [ stop ]
end 

to lookAtSwitches
  ; mSlope sitch can be handled at setup alone
  ; look at vetiverpm switch. 
  set D  (Kn / n * (sqrt mSlope) * (1 - 2 * vetiverpm * rootR))

  ; look at Lv (Root length density) switch
  set fad (Lv * 2 * pi * rootR) ; facial area density m^-1
  set thetaRoots        (Lv * pi * rootR * rootR)
  ;set thetaFilmMax      (fad * HuF)
  ;set thetaSatWithRoots (thetas - thetaRoots)
  ;set qFilmMax          (fad * VuHuF)     ; by this we underestimate
  
  set thetaFilmMax      (Lv * pi * rhc * rhc)  ;(fad * HuF)
  set thetaSatWithRoots (thetas - thetaFilmMax)
  set qFilmMax          (thetaFilmMax - thetaRoots) * VuMax   ; m s-1 
end 

; reset all simulation variables related to plots while retaining
; moisture and water levels obtained after the transient period.
; Idea is to run the simulation for sufficient time so that the
; model is in a steady flow state. Then export the world. Later to
; run another simulation, import the saved world, cleanSlate and go

to cleanSlate
  set totalRain 0
  set totalSurfaceWater 0
  set soilMoisture 0
  set filmMoisture 0
  set surfaceRunout 0
  set diffusedWater 0
  set filmedWater 0
  set filmRunout 0
  set waterLevelBeforeVetiver 0
  set waterLevelAtVetiver 0
  set bigBang ticks
  set simulatedTime 0
  set plotAt 0 ;((ticks - bigBang) * tickTime) + 60 ; plot at 60 s intervals

  clearPlots
end 

to setInputFromOutputForDiffuseFlow
  ask patchSetWithDiffusion [
    let qitemp 0.0
    ask patch pxcor (pycor + 1) [
      ifelse member? self rainFallsHerePS [
        set  KThetaFactor 
          ((h * 100) + (pycor * pDcm) - (zeroRef * 100)) 
          / ((pycor  * pDcm) - (zeroRef * 100)) 
        let Kseffective (Ks * KThetaFactor)
        ifelse (Kseffective < inFlow) [
          set qitemp Kseffective
        ][
          set qitemp inFlow
        ]
        set inFlow (inFlow - qitemp)
        set  diffusedWater (diffusedWater + (qitemp * tickTime * pDcm))
      ][
        set qitemp qo
      ]
    ]
    set qi qitemp
  ]
end 

to diffuseMoisture
  ; diffuse moisture according to unsteady unsaturated flow with
  ; hydraulic diffusivity (soil-water diffusivity)
  ; eqns (6) and (7) from [Nimmo] Unsaturated Zone Flow Processes

  ; diffuse accprding to steady state diffusion model
  ; eqn [3] in nimmo2010
  ; relationship between thetaDifus (of a patch) and outflow in head units
  ; assume patch with pDcm cm dimension and thetaDifus moisture content;
  ; assume q cm (head units) of water enters the patch
  ; thetaDash = thetaDifus + q/pDcm
  ; similarly if moixture content changes from thetaDifus to thetaDash
  ; q (water entering in cm of head units) = (thetaDash - theta)* pDcm

  ; assuming steady state flow
  ; qi' = qo' = -K(theta) * ((psi2-psi1')/delZ +1)
  ; assume qi' and psi2 at this t and all others at t-1, solve for psi1'
  ; new theta1'. then qo' from qi and the water to change theta1
  ; 28 sep 2011
  ; 

; psi is somewhat uncontrollable. better to go with thetaDifus which has
; small range and limits.
; example t1 = .074,kt1 = 8.48817e-15 dt1 = 5.65117e-08
; t2 = .07 
; patch 1 now gets qi = .0016. 
;assume qo = qi and qo is result of diffusion + gravity (prev values) and
; t2. get t1 new and recompute mositure change + outflow
; qi = qo = -Dt1 * (t2 - t1new)/dz + Kt1
; t1dash=t2 - dz*(kt1-qi)/dt1
; t1dash computed as 11325
; -dt1*(t2-t1dash)/dz + kt1 = .0015999999999999999
; so this scheme does not work

; assume change in moisture content too
; qi = -Dt1(t2-t1dash)/dz + Kt1 + (t1dash - t1)*dz
; t1dash = (qi + dt1*t2/dz -kt1 + t1*dz)/(dt1/dz + dz)
; gives t1dash as .07799999717439477757
; appears ok.
; so from prev values of D(theta1), theta2, theata1 and K(theta1) 
; calculate theta1Dash for all patches and qo
; screenshot  in archive/avet.nlogo.29sep1200.jpg
; 
  set soilMoisture 0.0
  ask patchSetWithDiffusion [
    ; diffusion
    let surfycor    0
    let waterOnSurf 0.0 ; cm
    ask patch pxcor 0 [set surfycor surfacepYcor]
    ask patch pxcor surfycor [set waterOnSurf (h * 100)]
    set KThetaFactor ((waterOnSurf + (surfycor + 1) * pDcm 
      - (zeroRef * 100)) / ((surfycor + 1) * pDcm - (zeroRef * 100)) )
    let effKTheta KTheta * KThetaFactor

    ;let theta1Dash 0.0
    let theta2 thetaResidual
    ifelse (pycor > (noVetiverRootsHere / pDimen) ) [
      ask patch pxcor (pycor - 1) [
        set theta2  thetaDifus
      ]
    ][
      set theta2 0.0
    ]
    set theta1Dash ((qi + (DTheta * theta2 / pDcm) - effKTheta 
      + (thetaDifus * pDcm / tickTime)) / ((DTheta / pDcm) + pDcm / tickTime))
    let thetaSatForThis thetas
		if member? self patchSetWithFilmFlow [
			set thetaSatForThis thetaSatWithRoots
		]
		if (theta1Dash > thetaSatForThis) [ set theta1Dash thetaSatForThis ]
  ]  
  ask patchSetWithDiffusion [
    ; it is better to postpone setting of thetaDifus, soilMoisture to
    ; another iteration over all the patches. This is similar to setting
    ; q_i from q_o of up patches at the previous instant for all
    ; patches before proceeding further
    set qo (qi - (theta1Dash - thetaDifus) * pDcm / tickTime)
    set thetaDifus theta1Dash
    set soilMoisture (soilMoisture + (thetaDifus * pDcm * pDcm))
    setOtherDiffusiveParams self
  ]
end 

to setInputFromOutputForFilmFlow
  ; set qipref from qopref of patch above or from the surface water

  let qTempPref 0.0

  ask patchSetWithFilmFlow [
    let qitemp 0.0
    set qTempPref (qFilmMax * aaf) 
    set qTempPref (qTempPref * 100 ); / (pDimen * 1)) ; head units cm/s
    ask patch pxcor (pycor + 1) [
      ifelse member? self rainFallsHerePS [
        ifelse (qTempPref < inFlow) [
          set qitemp qTempPref
        ][
          set qitemp inFlow
        ]
        set inFlow (inFlow - qitemp)
        set  filmedWater (filmedWater + (qitemp * tickTime * pDcm)  )
      ][
        set qitemp qopref
      ]
    ]
    set qipref qitemp
  ]
  ; done set qipref from qopref of patch above or from the surface water
end 

; to flow water as thin surface film on the vertical root surfaces
; ref Theory for Source-Responsive and Free-Surface Film Modelling 
; of Unsaturated Flow. John R. Nimmo; Vadose Zone J. 2010, Vol 9 p295-306

to  surfaceFilmFlowOnRoots

  ; film flow across the patches
  set filmMoisture 0.0
  ask patchSetWithFilmFlow [
    ;ifelse (thetapref >= thetaFilmMax) [
    ;  set qopref qipref
    ;  ; later check whether we have to dip into the film
    ;][
      let qiprefAsRatio (qipref * tickTime / pDcm )
      let qopreftempAsRatio  0.0
      set thetapref (thetapref + qiprefAsRatio)
      if (thetapref > thetaFilmMax) [
        set qopreftempAsRatio (thetapref - thetaFilmMax)
        set thetapref thetaFilmMax
      ]
      set qopref (qopreftempAsRatio * pDcm / tickTime)
    ;]
    set filmMoisture (filmMoisture + (thetapref * pDcm * pDcm))
    if (pycor = (noVetiverRootsHere / pDimen)) [
      set filmRunout (filmRunout + (qopref * tickTime * pDcm))
    ]
  ]
  ; end film flow across the patches
end 

to-report surfacePatchFinder [aPatch stepsDown]
  if not member? aPatch rainFallsHerePS [report nobody]
  let downPatch nobody
  ask aPatch [
    if ((pxcor + stepsDown) <= max-pxcor) and
        ((pxcor + stepsDown) >= 0) [
      let yOfAnotherSurfPatch 0
      ask patch (pxcor + stepsDown) 0 [
        set yOfAnotherSurfPatch surfacepYcor
      ]
      set downPatch patch (pxcor + stepsDown) yOfAnotherSurfPatch
      if not member? downPatch rainFallsHerePS [
        set downPatch nobody]
    ]
  ]
  report downPatch
end 

to setInputFromOutputForSurfaceFlow
  ask rainFallsHerePS [
    let ourInflow 0
    set hprev h
    let upPatch (surfacePatchFinder self  ( - 1))
    ifelse upPatch != nobody [
      ask upPatch [ set ourInflow outFlow ]
    ][
      ; 0 th patch. Calculate rain in hinterland and treat as inflow
      let ourIn ( rainfallIntensity / 3600) * hinterland / pDimen 
        * (0.99 + random-float 0.02)  ; cm/sec
      set ourInflow ourIn
      set totalRain (totalRain + (ourInflow * tickTime * pDcm))
    ]
    set inFlow ourInflow
  ]
end 

to surfaceRunoff
  let QiTemp 0
  let Q 0

  set totalSurfaceWater 0
  let vpatch rainFallsHerePS with [pxcor = vetRpxcor] 

  ask vpatch [
    set QiTemp inFlow  ; water inflow head units cm/s
    set Q (QiTemp * pDimen * 1 / 100) ; m^3/s
    set hv 0
    set Rv  (0.5 / vetiverpm - vetiverShootRadius)
    set hm 0
    if  (mSlope > 0) and (Q > 0) [
      set hm exp (0.6 * ln (Q * (n / kn) / sqrt (mSlope) ))]
    ifelse (vetiverpm = 1) [
      set hv hm
    ][
      set hv (heightFromQForVetiverHedge Q)
    ]
    set h hv
    set deadHeight (hv - hm + (pycor * pDimen))
    set waterLevelAtVetiver h * 100
    set totalSurfaceWater totalSurfaceWater + (waterLevelAtVetiver * pDcm)
    set backwaterLength (hv - hm) / (tan slpAngle) 
    set backwaterStartPos (vetiverXPos - backwaterLength)
  ]
  ask rainFallsHerePS [
    set QiTemp inFlow  ; water inflow head units cm/s
    set Q (QiTemp * pDimen * 1 / 100) ; m^3/s
    set hm 0
    if  (mSlope > 0) and (Q > 0) [
      set hm exp (0.6 * ln (Q * (n / kn) / sqrt (mSlope) ))
    ]
    set v 0.0
    if (mSlope > 0) and (hm > 0) [
      set v (kn / n) * exp (0.6667 * ln hm ) * sqrt (mSlope)]
    ifelse surfaceRun [
      if pxcor != vetRpxcor [
        ; not at vetiver hedge
        ifelse pxcor < vetRpxcor [
          let yGround (pycor) * pDimen
          ifelse deadHeight > yGround [
            set h (deadHeight + hm - yGround)
          ][
             set h hm
          ]
        ][
           set h hm
        ]
      ]
      set outFlow (inFlow + (hprev - h) / tickTime)  ; steady state condition
    ][
      set outFlow 0
    ]
    set totalSurfaceWater totalSurfaceWater + (h * 100 * pDcm)
    if pxcor =  (vetLpxcor -  5) [ 
      set waterLevelBeforeVetiver (h * 100) ; cm
    ]
    if pxcor =  (vetRpxcor + 5) [
      set waterLevelAfterVetiver (h * 100) ; cm
    ]
    if pxcor =  max-pxcor [
      set surfaceRunout surfaceRunout + outFlow * tickTime * pDcm
    ]
  ]
end 

to showSurfaceWaterLevel
  let thePatchAbove nobody
  ask rainFallsHerePS [
    let wabove (h * 100) ;
    if wabove < 0 [set wabove 0]
    let exitloop false
    let waterFilledPatches ceiling (wabove / pDcm)
    let patchAbove 0
    let waterForRest wabove
    if debug and (wabove > 200) [ ; enough to drown
      show word "water " wabove
      show word "fills " waterFilledPatches
    ]
    repeat waterFilledPatches 
    [ set thePatchAbove (patch-at 0 patchAbove)
      if thePatchAbove != nobody
      [ ask thePatchAbove [
          ifelse (waterForRest > pDcm)
          [ set waterInPatch pDcm
            set waterForRest waterForRest - pDcm
          ][
            set waterInPatch waterForRest
            set waterForRest 0
          ]
          set pcolor (blue  + 5 - (waterInPatch * 5 / pDcm))
        ]
        set patchAbove (patchAbove + 1)
      ]
    ]
    while [not exitloop] [
      set thePatchAbove (patch-at 0 patchAbove)
      ifelse thePatchAbove = nobody [
        set exitloop true
      ][
        ask patch-at 0 patchAbove [ 
          ifelse (pcolor = originalColor) [
            set exitloop true
          ][
            set waterInPatch 0
            set pcolor originalColor 
            set patchAbove (patchAbove + 1)
          ]
        ]
      ]
    ]
  ]    
end 

to showMoistureLevel
  let thePatchAbove nobody
  ask patchSetWithDiffusion [
		ifelse member? self patchSetWithFilmFlow [
    	set pcolor (originalColor  + 2 - ((thetaDifus + thetapref) * 5  
			/ (thetas + thetaFilmMax - thetaRoots)))
		][
    	set pcolor (originalColor  + 2 - ((thetaDifus ) * 5  
			/ thetas ))
		]
  ]    
end 
;;; plotting procedures

to my-setup-plots
  set-current-plot "Total water flowed"
  set-plot-x-range 0 simulateDuration / 2
  set-plot-y-range 0 100; water contents will be expressed as % of totalRain
  ;
  set-current-plot "water level"
  set-plot-x-range 0 10 ; simulateDuration / 2
  set-plot-y-range 0 2 ; 4 ;(10 * ceiling (rainfallIntensity * 10 / 60 / 10))
end 

to clearPlots
  ;; update averaged rain
  set-current-plot "water at vetiver"
  clear-plot
  set-current-plot "Total water flowed"
  clear-plot
  set-current-plot "water level"
  clear-plot
  set-current-plot "diffusive and film flow"
  clear-plot
  set-current-plot "thetaDifus vs depth"
  clear-plot
  set-current-plot "water diffused filmed"
  clear-plot
  set-current-plot "Moisture stores"
  clear-plot
end 

to my-update-plots
  ;; update averaged rain
  set-current-plot "Total water flowed"

  set-current-plot-pen "rainIn"
  plotxy (simulatedTime) 100  ; (totalRain * 100 / totalRain)

  set-current-plot-pen "surfaceOut"
  plotxy (simulatedTime) (surfaceRunout * 100 / totalRain)

  ;set-current-plot-pen "filmOut"
  ;plotxy (simulatedTime) (filmRunout * 100 / totalRain)
  
  set-current-plot-pen "filmIn"
  plotxy (simulatedTime) (filmedWater * 100 / totalRain)

  set-current-plot-pen "diffuseIn"
  plotxy (simulatedTime) (diffusedWater * 100 / totalRain )

  ;; plot water height before and after hedgerow
  set-current-plot "water level"
  set-current-plot-pen "hedge-"
  plotxy (simulatedTime) waterLevelBeforeVetiver * 10
  set-current-plot-pen "hedge+"
  plotxy (simulatedTime) waterLevelAfterVetiver * 10
  set-current-plot-pen "hedge"
  plotxy (simulatedTime) waterLevelAtVetiver   * 10

  ;; plot near vetiver hedgerow
  if ((ticks - bigBang) * tickTime) >= plotAt [
    update-waterAtVetiver
    update-moistureDiffusive
    update-DiffusiveAndFilmFlowFluxes
    update-DiffusiveAndFilmFlowWater
    plotMoistureStores
    set plotAt ((ticks - bigBang) * tickTime)  + 60 
  ]
end 

to update-waterAtVetiver 
  ;; plot near vetiver hedgerow
  set-current-plot "water at vetiver"
  clear-plot
  let ystart 0
  ; leftPatchForPlot ; 
  ; rightPatchForPlot ; 
  ask patch rightPatchForPlot 0 [set ystart (surfacepYcor  * pDimen * 100)]
  let ipx leftPatchForPlot
  repeat (rightPatchForPlot - leftPatchForPlot + 1) [
    let thePatchycor  0
    ask patch ipx 0 [ set thePatchycor surfacepYcor]
    ask patch ipx thePatchycor [
      if (pxcor >= leftPatchForPlot ) and (pxcor <= rightPatchForPlot) [
        let grlevel ((pycor + 1) * pDimen   * 100 - ystart)
        ifelse pxcor = vetRpxcor [
          set-current-plot-pen "vetiver" ; 
          plotxy (pxcor + 1  - leftPatchForPlot) grlevel
          set-current-plot-pen "vetiver"
          plotxy (pxcor + 1  - leftPatchForPlot) (h * 100 + grlevel)
          set-current-plot-pen "waterstep"
          plotxy (pxcor + 1  - leftPatchForPlot) (h * 100 + grlevel)
        ][
          set-current-plot-pen "waterstep"
          plotxy (pxcor + 1  - leftPatchForPlot) (h * 100 + grlevel)
          set-current-plot-pen "ground"
          plotxy (pxcor + 1  - leftPatchForPlot)   grlevel
          set-current-plot-pen "deadLevel"
          let dhplot deadHeight * 100 - ystart
          if (dhplot  > grlevel) and (pxcor <= vetRpxcor + 1) [
            plotxy (pxcor + 1  - leftPatchForPlot)    dhplot
          ]
        ]
      ]
    ] 
    set ipx ipx + 1
  ]
end 

to  update-moistureDiffusive
  set-current-plot "thetaDifus vs depth"
  clear-plot
  let xpatch 12
  let ymax 0
  ask patch xpatch 0 [ set ymax ( surfacepYcor - 1)]
  let ythis  (noVetiverRootsHere / pDimen)
  let topz (ymax * pDcm + 0.5 * pDcm)
  let ydepth 0.0
  let thetaatyx 0.0
  while [ythis <= ymax] [
    ask patch xpatch ythis [
      set ydepth z
      set thetaatyx thetaDifus
    ]
    set-current-plot-pen "atxcor12"
    plotxy thetaatyx ydepth ;(topz -  ydepth)
    set ythis (ythis + 1)
  ]
end 

to  update-DiffusiveAndFilmFlowFluxes
  set-current-plot "diffusive and film flow"
  let xpatch (vetRpxcor - 2)
  let ypatch 30
  ask surfPatchInfoRow with [pxcor = xpatch] [
    set ypatch (surfacepYcor - 5)
  ]
  let qd 0.0
  let qf 0.0
  ask patch xpatch ypatch [
    set qf qopref
    set qd qo
  ]
  set-current-plot-pen "diffuse"
  plotxy (simulatedTime) (qd * 10 * 60)
  ;plot-pen-reset
  set-current-plot-pen "film"
  plotxy (simulatedTime) (qf * 10 * 60)
end 

to  update-DiffusiveAndFilmFlowWater
  set-current-plot "water diffused filmed"

  set-current-plot-pen "diffuse"
  plotxy (simulatedTime) (diffusedWater )

  set-current-plot-pen "film"
  plotxy (simulatedTime) (filmedWater )
end 

to plotMoistureStores
  set-current-plot "Moisture stores"

  set-current-plot-pen "surface"
  plotxy (simulatedTime) (totalSurfaceWater)

  set-current-plot-pen "soil"
  plotxy (simulatedTime) (soilMoisture)

  set-current-plot-pen "root"
  plotxy (simulatedTime) (filmMoisture)
end