PEAT VOXEL Model

;; PEAT VOXEL MODEL
;;
;; A voxel-based model of peatland hydrology
;;
;; Aryo Adhi Condro @ 2017
;;
;; Supervisors: Idung Risdiyanto, S. Si, M.Sc. & Prof. Dr. Hidayat Pawitan
;;
;; GFM - IPB, Bogor, Indonesia
;;
;; Updated: 10-May-2017


breed[waters water]

extensions
[
  array
  csv
  gis
  palette
  pathdir
]

waters-own
[
  water-height
]

globals
[
  input-array-asc
  input-array-csv
  input-path-asc
  input-path-csv

  output-path-asc
  output-path-csv
  output-path-png

  lookup-parameters-peat
  lookup-parameters-climate
  lookup-parameters-land
  lookup-parameters-scenario

  map-attributes-n-cols
  map-attributes-n-rows
  map-attributes-corner-xll
  map-attributes-corner-yll
  map-attributes-cell-size-degree
  map-attributes-cell-size-meter
  map-attributes-cell-area-ha

  timeseries-precipitation ;;mm
  timeseries-potential-evapotranspiration ;;mm

  ;;peat
  parameter-max-peat-depth ;; meter
  parameter-z-size ;; meter
  parameter-peat-layer-number ;;layers

  parameter-hydraulic-conductivity-mean-acrotelm ;; meter per day
  parameter-hydraulic-conductivity-mean-catotelm ;; meter per day
  parameter-hydraulic-conductivity-cv-acrotelm ;; %
  parameter-hydraulic-conductivity-cv-catotelm ;; %

  parameter-initial-water-height-mean-acrotelm ;;mm
  parameter-initial-water-height-cv-acrotelm ;;%
  parameter-initial-water-height-mean-catotelm ;;mm
  parameter-initial-water-height-cv-catotelm ;;%

  parameter-hydraulic-conductivity-min-acrotelm ;; meter per day
  parameter-hydraulic-conductivity-min-catotelm ;; meter per day
  parameter-hydraulic-conductivity-max-acrotelm ;; meter per day
  parameter-hydraulic-conductivity-max-catotelm ;; meter per day

  parameter-max-depth-acrotelm-mean ;;m
  parameter-max-depth-acrotelm-cv ;;%
  parameter-max-depth-acrotelm-min ;;m
  parameter-max-depth-acrotelm-max ;;m
  parameter-saturated-capacity-mean-acrotelm ;;mm
  parameter-saturated-capacity-cv-acrotelm ;;%
  parameter-saturated-capacity-min-acrotelm ;;mm
  parameter-saturated-capacity-max-acrotelm ;;mm
  parameter-saturated-capacity-mean-catotelm ;;mm
  parameter-saturated-capacity-cv-catotelm ;;%
  parameter-saturated-capacity-min-catotelm ;;mm
  parameter-saturated-capacity-max-catotelm ;;mm


  parameter-field-capacity-mean-acrotelm ;;mm
  parameter-field-capacity-cv-acrotelm ;;%
  parameter-field-capacity-min-acrotelm ;;mm
  parameter-field-capacity-max-acrotelm ;;mm
  parameter-field-capacity-mean-catotelm ;;mm
  parameter-field-capacity-cv-catotelm ;;%
  parameter-field-capacity-min-catotelm ;;mm
  parameter-field-capacity-max-catotelm ;;mm

  parameter-overland-flow-velocity ;;m/day
  parameter-overland-flow-routing-distance-per-day-meter;;meter
  parameter-overland-flow-routing-distance-per-day-cell;;xy of voxel

  ;;climate
  parameter-precipitation-gross
  parameter-potential-evapotranspiration

  ;;scenario
  parameter-scenario-number

  ;;voxels
  voxels-zone
  voxels-surface
  voxels-peat
  voxels-acrotelm
  voxels-catotelm

  ;;user message
  message-status

  ;;run
  run-day
  run-day-number
]

patches-own
[
  ;;voxels
  patch-base?
  patch-layer
  patch-surface?
  patch-zone?

  ;;peat
  patch-peat?
  patch-acrotelm?
  patch-depth
  patch-peat-depth
  patch-peat-hydraulic-conductivity
  patch-water-height
  patch-peat-water-diffusion-darcy
  patch-peat-water-diffusion-outflow-target-voxels
  patch-peat-water-diffusion-outflow-target-voxels-fraction
  patch-peat-water-diffusion-inflow
  patch-peat-water-diffusion-inflow-list
  patch-peat-water-diffusion-inflow-delta-list
  patch-peat-water-diffusion-outflow-delta-list
  patch-peat-water-diffusion-outflow
  patch-peat-water-diffusion-inflow-source-voxels-list
  patch-peat-saturated-capacity
  patch-peat-field-capacity
  patch-water-height-peat


  ;;elevation
  patch-elevation

  ;;land
  patch-land
  patch-land-pet-correction-factor ;;%
  patch-land-interception-potential ;;mm
  patch-land-drought-resistance ;;%
  patch-land-interception-effect-on-transpiration ;;%

  ;;climate
  patch-precipitation-gross


  ;;reporting maps
  patch-water-height-ascii

  ;;observation
  patch-observation-id
]

to-report gis-read-map[m]
 let r gis:load-dataset m
 gis:set-world-envelope gis:envelope-of r
 report r
end 

to initialise-input-parameters
  csv-read-input

  set parameter-z-size item 0 item 1 lookup-parameters-peat ;; meter
  set parameter-initial-water-height-mean-acrotelm item 1 item 1 lookup-parameters-peat;;mm
  set parameter-initial-water-height-cv-acrotelm item 2 item 1 lookup-parameters-peat;;%
  set parameter-initial-water-height-mean-catotelm item 3 item 1 lookup-parameters-peat;;mm
  set parameter-initial-water-height-cv-catotelm item 4 item 1 lookup-parameters-peat;;%
  set parameter-max-depth-acrotelm-mean item 5 item 1 lookup-parameters-peat;;m
  set parameter-max-depth-acrotelm-cv item 6 item 1 lookup-parameters-peat;;%
  set parameter-max-depth-acrotelm-min item 7 item 1 lookup-parameters-peat;;m
  set parameter-max-depth-acrotelm-max item 8 item 1 lookup-parameters-peat;;m
  set parameter-overland-flow-velocity item 9 item 1 lookup-parameters-peat;;m/day

  set parameter-overland-flow-routing-distance-per-day-meter parameter-overland-flow-velocity;; meter/day
  set parameter-overland-flow-routing-distance-per-day-cell round(parameter-overland-flow-routing-distance-per-day-meter / map-attributes-cell-size-meter)
  set parameter-scenario-number item 0 item 1 lookup-parameters-scenario


  set timeseries-precipitation item 1 lookup-parameters-climate;;mm
  set timeseries-potential-evapotranspiration item 2 lookup-parameters-climate;;mm

  set run-day 1
  set run-day-number length (item 0 lookup-parameters-climate)
end 

to initialise-declare-input
 set input-path-asc "data/input/asc/"
 set input-path-csv "data/input/csv/"

 set input-array-asc array:from-list
 [
   "map-dem.asc"
   "map-land.asc"
   "map-observation.asc"
   "map-peat-depth.asc"
 ]

 set input-array-csv array:from-list
 [
   "parameters-peat.csv"
   "parameters-climate.csv"
   "parameters-land.csv"
   "parameters-scenario.csv"
 ]
end 

to initialise-declare-output
  set output-path-csv word"/data/output/" word "run-" word parameter-scenario-number "/csv/"
  pathdir:create output-path-csv

  set output-path-png word"/data/output/" word "run-" word parameter-scenario-number "/png/"
  pathdir:create output-path-png

  if(report-ascii-map?)
  [
    set output-path-asc word"/data/output/" word "run-" word parameter-scenario-number "/asc/"
    pathdir:create output-path-asc
  ]

  output-write-main-header-csv1 word output-path-csv "output-gwl-average.csv"
  output-write-main-header-csv2 word output-path-csv "output-gwl-observation1.csv"
  output-write-main-header-csv3 word output-path-csv "output-gwl-observation2.csv"
  output-write-main-header-csv4 word output-path-csv "output-gwl-observation3.csv"
  output-write-main-header-csv5 word output-path-csv "output-gwl-observation4.csv"
  output-write-main-header-csv6 word output-path-csv "output-gwl-observation5.csv"

  output-write-main-header-csv7 word output-path-csv "output-wc-average.csv"
  output-write-main-header-csv8 word output-path-csv "output-wc-observation1.csv"
  output-write-main-header-csv9 word output-path-csv "output-wc-observation2.csv"
  output-write-main-header-csv10 word output-path-csv "output-wc-observation3.csv"
  output-write-main-header-csv11 word output-path-csv "output-wc-observation4.csv"
  output-write-main-header-csv12 word output-path-csv "output-wc-observation5.csv"

  output-write-main-header-csv13 word output-path-csv "output-wc-landcover-forest.csv"
  output-write-main-header-csv14 word output-path-csv "output-wc-landcover-shrub.csv"
  output-write-main-header-csv15 word output-path-csv "output-wc-landcover-burntland.csv"
  output-write-main-header-csv16 word output-path-csv "output-wc-landcover-bareland.csv"

  output-write-main-header-csv17 word output-path-csv "output-gwl-landcover-forest.csv"
  output-write-main-header-csv18 word output-path-csv "output-gwl-landcover-shrub.csv"
  output-write-main-header-csv19 word output-path-csv "output-gwl-landcover-burntland.csv"
  output-write-main-header-csv20 word output-path-csv "output-gwl-landcover-bareland.csv"
end 

to gis-read-map-attributes
 let m word input-path-asc array:item input-array-asc 0

 file-open m
  set map-attributes-n-cols read-from-string remove "ncols " file-read-line;;cells
  set map-attributes-n-rows read-from-string remove "nrows " file-read-line;;cells
  set map-attributes-corner-xll read-from-string remove "xllcorner " file-read-line;;degree
  set map-attributes-corner-yll read-from-string remove "yllcorner " file-read-line;;degree
  set map-attributes-cell-size-degree read-from-string remove "cellsize " file-read-line;;degree
 file-close

 set map-attributes-cell-size-meter map-attributes-cell-size-degree * 111000 ;;meter
 set map-attributes-cell-area-ha (map-attributes-cell-size-meter ^ 2) / 10000 ;;ha
end 

to gis-resize-world
 let max-z round max[patch-elevation]of patches
 gis-read-map-attributes
 resize-world 0 (map-attributes-n-cols - 1) 0 (map-attributes-n-rows - 1) 0 (parameter-peat-layer-number)
 set-patch-size 18

 orbit-down 45
 orbit-left 210
end 

to gis-read-input-maps
  ask patches
  [
    set patch-elevation -9999
    set patch-land -9999
    set patch-observation-id -9999
    set patch-peat-depth -9999
  ]

  gis:apply-raster (gis-read-map word input-path-asc array:item input-array-asc 0) patch-elevation
  gis:apply-raster (gis-read-map word input-path-asc array:item input-array-asc 1) patch-land
  gis:apply-raster (gis-read-map word input-path-asc array:item input-array-asc 2) patch-observation-id
  gis:apply-raster (gis-read-map word input-path-asc array:item input-array-asc 3) patch-peat-depth

  set parameter-max-peat-depth round((max[patch-peat-depth]of patches) / 100);; meter
  set parameter-peat-layer-number round(parameter-max-peat-depth / parameter-z-size)
end 

to initialise-land
  ask voxels-surface
  [
    let i position patch-land item 1 lookup-parameters-land
    set patch-land-pet-correction-factor item i item 2 lookup-parameters-land ;;%
    set patch-land-interception-potential item i item 3 lookup-parameters-land ;;mm
    set patch-land-drought-resistance item i item 4 lookup-parameters-land ;;%
    set patch-land-interception-effect-on-transpiration item i item 5 lookup-parameters-land ;;%

    set parameter-hydraulic-conductivity-mean-acrotelm item i item 6 lookup-parameters-land ;;m/d
    set parameter-hydraulic-conductivity-cv-acrotelm item i item 7 lookup-parameters-land ;;%
    set parameter-hydraulic-conductivity-min-acrotelm item i item 8 lookup-parameters-land ;;m/d
    set parameter-hydraulic-conductivity-max-acrotelm item i item 9 lookup-parameters-land ;;m/d

    set parameter-hydraulic-conductivity-mean-catotelm item i item 10 lookup-parameters-land ;;m/d
    set parameter-hydraulic-conductivity-cv-catotelm item i item 11 lookup-parameters-land ;;%
    set parameter-hydraulic-conductivity-min-catotelm item i item 12 lookup-parameters-land ;;m/d
    set parameter-hydraulic-conductivity-max-catotelm item i item 13 lookup-parameters-land ;;m/d

    set parameter-field-capacity-mean-acrotelm item i item 14 lookup-parameters-land;;mm
    set parameter-field-capacity-cv-acrotelm item i item 15 lookup-parameters-land;;%
    set parameter-field-capacity-min-acrotelm item i item 16 lookup-parameters-land;;mm
    set parameter-field-capacity-max-acrotelm item i item 17 lookup-parameters-land;;mm

    set parameter-field-capacity-mean-catotelm item i item 18 lookup-parameters-land;;mm
    set parameter-field-capacity-cv-catotelm item i item 19 lookup-parameters-land;;%
    set parameter-field-capacity-min-catotelm item i item 20 lookup-parameters-land;;mm
    set parameter-field-capacity-max-catotelm item i item 21 lookup-parameters-land;;mm

    set parameter-saturated-capacity-mean-acrotelm item i item 22 lookup-parameters-land;;mm
    set parameter-saturated-capacity-cv-acrotelm item i item 23 lookup-parameters-land;;%
    set parameter-saturated-capacity-min-acrotelm item i item 24 lookup-parameters-land;;mm
    set parameter-saturated-capacity-max-acrotelm item i item 25 lookup-parameters-land;;mm

    set parameter-saturated-capacity-mean-catotelm item i item 26 lookup-parameters-land;;mm
    set parameter-saturated-capacity-cv-catotelm item i item 27 lookup-parameters-land;;%
    set parameter-saturated-capacity-min-catotelm item i item 28 lookup-parameters-land;;mm
    set parameter-saturated-capacity-max-catotelm item i item 29 lookup-parameters-land;;mm
  ]
end 

to initialise-peat-properties
  ask patches
  [
    set patch-peat-field-capacity -9999
    set patch-peat-hydraulic-conductivity -9999
    set patch-peat-saturated-capacity -9999
  ]

  ask voxels-acrotelm
  [
    set patch-peat-hydraulic-conductivity ;;meter per day
       (min list parameter-hydraulic-conductivity-max-acrotelm
       (max list parameter-hydraulic-conductivity-min-acrotelm
       (random-normal parameter-hydraulic-conductivity-mean-acrotelm
       ((parameter-hydraulic-conductivity-cv-acrotelm / 100) * parameter-hydraulic-conductivity-mean-acrotelm))))

    set patch-peat-field-capacity ;;mm
        (min list parameter-field-capacity-max-acrotelm
        (max list parameter-field-capacity-min-acrotelm
        (random-normal parameter-field-capacity-mean-acrotelm
        ((parameter-field-capacity-cv-acrotelm / 100) * parameter-field-capacity-mean-acrotelm))))

    set patch-peat-saturated-capacity ;;mm
        (min list parameter-saturated-capacity-max-acrotelm
        (max list parameter-saturated-capacity-min-acrotelm
        (random-normal parameter-saturated-capacity-mean-acrotelm
        ((parameter-saturated-capacity-cv-acrotelm / 100) * parameter-saturated-capacity-mean-acrotelm))))
  ]

  ask voxels-catotelm
  [
    set patch-peat-hydraulic-conductivity ;;meter per day
       (min list parameter-hydraulic-conductivity-max-catotelm
       (max list parameter-hydraulic-conductivity-min-catotelm
       (random-normal parameter-hydraulic-conductivity-mean-catotelm
       ((parameter-hydraulic-conductivity-cv-catotelm / 100) * parameter-hydraulic-conductivity-mean-catotelm))))

    set patch-peat-field-capacity ;;mm
        (min list parameter-field-capacity-max-catotelm
        (max list parameter-field-capacity-min-catotelm
        (random-normal parameter-field-capacity-mean-catotelm
        ((parameter-field-capacity-cv-catotelm / 100) * parameter-field-capacity-mean-catotelm))))

     set patch-peat-saturated-capacity ;;mm
        (min list parameter-saturated-capacity-max-catotelm
        (max list parameter-saturated-capacity-min-catotelm
        (random-normal parameter-saturated-capacity-mean-catotelm
        ((parameter-saturated-capacity-cv-catotelm / 100) * parameter-saturated-capacity-mean-catotelm))))
  ]

  ask voxels-acrotelm
  [
      set patch-water-height ;;mm
        (min list patch-peat-saturated-capacity
        (max list 0
        (random-normal parameter-initial-water-height-mean-acrotelm
        ((parameter-initial-water-height-cv-acrotelm / 100) * parameter-initial-water-height-mean-acrotelm))))
  ]

  ask voxels-catotelm
  [
    set patch-water-height ;;mm
        (min list patch-peat-saturated-capacity
        (max list 0
        (random-normal parameter-initial-water-height-mean-catotelm
        ((parameter-initial-water-height-cv-catotelm / 100) * parameter-initial-water-height-mean-catotelm))))
  ]
end 

to initialise-patches
  ask patches
  [
    set patch-base? false
    set patch-surface? false
    set patch-zone? false
    set patch-water-height -9999
    set patch-peat? false
    set patch-depth -9999
    set patch-acrotelm? false
  ]

  ask patches
  [
    set patch-base? pzcor = 0 and (patch-elevation <= 0 or patch-elevation >= 0) and (patch-land > 0) and (patch-peat-depth <= 0 or patch-peat-depth >= 0)
  ]

  let z parameter-peat-layer-number + 1
  let z-list n-values z[?]

  ask patches with[patch-base?]
  [
    let x pxcor
    let y pycor

    let dem patch-elevation
    let land patch-land
    let obs patch-observation-id
    let peat-depth patch-peat-depth

    foreach z-list
    [
      ask patch x y ?
      [
        set patch-zone? true
        set patch-elevation dem
        set patch-land land
        set patch-observation-id obs
        set patch-peat-depth peat-depth
        set patch-layer parameter-peat-layer-number - ?
        set patch-depth round(patch-layer * parameter-z-size * 10) / 10
      ]
    ]

    ask patch x y (last z-list)
    [
      set patch-surface? true
    ]
  ]

  ask patches
  [
    set patch-peat? patch-zone? and patch-layer >= 1 and patch-layer <= parameter-peat-layer-number
        and patch-depth <= (patch-peat-depth / 100)
  ]

  set voxels-zone patches with[patch-zone?]
  set voxels-surface patches with[patch-surface?]
  set voxels-peat patches with[patch-peat?]

  ask voxels-peat
  [
    let d
        max list parameter-max-depth-acrotelm-min
        (min list parameter-max-depth-acrotelm-max
        (random-normal parameter-max-depth-acrotelm-mean
        (parameter-max-depth-acrotelm-mean * (parameter-max-depth-acrotelm-cv / 100))))

    let rd (round (d * 10)) / 10

    set patch-acrotelm? (patch-depth > 0 and patch-depth <= rd)
  ]

  set voxels-acrotelm voxels-peat with[patch-acrotelm?]
  set voxels-catotelm voxels-peat with[not patch-acrotelm?]

  ask voxels-surface
  [
    set patch-water-height 0
  ]
end 

to gis-display-voxels
  ask voxels-zone
  [
    set pcolor palette:scale-gradient palette:scheme-colors "Divergent" "Spectral" 11 patch-water-height 0 100 ;;more blue more water
  ]
end 

to run-setup
  set message-status ""
  ca
  set message-status "initialising the model... please, wait..."
  initialise-declare-input
  gis-read-map-attributes
  initialise-input-parameters
  gis-read-input-maps
  gis-resize-world
  initialise-patches
  initialise-declare-output
  initialise-land
  initialise-peat-properties
  update-scenario-number
  reset-ticks
  set message-status "model initialisation is done! now, you can simulate!"
end 

to dynamic-flow-surface-water
  let target min-one-of neighbors with[patch-surface?] [patch-elevation + (sum [water-height]of waters-here)]

  if[patch-elevation + (sum [water-height]of waters-here)] of target < (patch-elevation + (sum [water-height]of waters-here))
  [
    move-to target
  ]
end 

to dynamic-simulate-surface-water
   set parameter-precipitation-gross item(run-day - 1) timeseries-precipitation
   set parameter-potential-evapotranspiration item(run-day - 1) timeseries-potential-evapotranspiration


   ask voxels-surface
   [
     set patch-precipitation-gross parameter-precipitation-gross
   ]

   ask voxels-surface
   [
     let x pxcor
     let y pycor
     let z pzcor
     let unsaturated-water-zone patches with[pxcor = x and pycor = y and patch-acrotelm?]

     let h mean[patch-water-height]of unsaturated-water-zone
     let saturated-capacity mean[patch-peat-saturated-capacity]of unsaturated-water-zone
     let field-capacity mean[patch-peat-field-capacity]of unsaturated-water-zone


     let interception-actual patch-land-interception-potential * ( 1 - exp (-1 * (patch-precipitation-gross / patch-land-interception-potential)))
     let evapotranspiration-corected (patch-land-pet-correction-factor / 100) * parameter-potential-evapotranspiration

     let relative-water-availability ifelse-value (field-capacity > 0)
     [
        min list 1 (h / ((patch-land-drought-resistance / 100) * field-capacity))
     ]
     [
       1
     ]

     let evapotranspiration-actual relative-water-availability * (max list 0 (evapotranspiration-corected - ((patch-land-interception-effect-on-transpiration / 100) * interception-actual)))

     ask unsaturated-water-zone
     [
       set patch-water-height max list 0 (patch-water-height - evapotranspiration-actual)
     ]


     ;; elaborate infiltration (MaxInf*t24)

     let infiltration-capacity max list 0 (saturated-capacity - h)
     let total-water-gross max list 0 (patch-precipitation-gross + patch-water-height - interception-actual)
     let infiltration min list total-water-gross infiltration-capacity

     ask unsaturated-water-zone
     [
       set patch-water-height min list patch-peat-saturated-capacity (patch-water-height + infiltration)
     ]

     let total-water-net max list 0 (total-water-gross - infiltration)

     if(total-water-net > 0)
     [
       sprout-waters 1
       [
         ht
         set water-height total-water-net
       ]
     ]
   ]

   ask voxels-surface
   [
     set patch-water-height 0
   ]


   if(any? waters)
   [
     repeat parameter-overland-flow-routing-distance-per-day-cell ;; iteration = routing distance of each cell per day
     [
       ask waters[dynamic-flow-surface-water]
        ask patches with[not patch-surface?]
       [
         if(any? waters-here)
         [
           ask waters-here [die]
         ]
       ]
     ]
   ]

   ask voxels-surface
   [
     if(any? waters-here)
     [
       set patch-water-height sum[water-height]of waters-here
     ]
   ]

  ask waters[die]
end 

to dynamic-diffuse-peat-water
  ;;initialisation
  ask voxels-peat
  [
    set patch-peat-water-diffusion-darcy 0
    set patch-peat-water-diffusion-inflow 0
    set patch-peat-water-diffusion-outflow 0
    set patch-peat-water-diffusion-inflow-list[]
    set patch-peat-water-diffusion-inflow-source-voxels-list []
    set patch-peat-water-diffusion-inflow-delta-list[]
    set patch-peat-water-diffusion-outflow-delta-list[]
  ]

  ask voxels-peat
  [
    let h patch-water-height
    let n neighbors6 with[patch-peat?]
    let z pzcor
    let source self

    ask n
    [
      let d ifelse-value(pzcor = z) ;;mm
         [
           map-attributes-cell-size-meter ;; meter
         ]
         [
           parameter-z-size ;; meter
         ]

      set patch-peat-water-diffusion-darcy
        -1 * patch-peat-hydraulic-conductivity * ((h - patch-water-height) / d)
    ]

    set patch-peat-water-diffusion-outflow-target-voxels n with[patch-peat-water-diffusion-darcy < 0]

    if(any? patch-peat-water-diffusion-outflow-target-voxels)
    [
      let total-outflow sum[-1 * patch-peat-water-diffusion-darcy]of patch-peat-water-diffusion-outflow-target-voxels
      let outflow min list patch-water-height total-outflow

      set patch-peat-water-diffusion-outflow ;;actual outflow from source voxels
          outflow

      ask patch-peat-water-diffusion-outflow-target-voxels
      [
        set patch-peat-water-diffusion-outflow-target-voxels-fraction
            (-1 * patch-peat-water-diffusion-darcy) / total-outflow

        set patch-peat-water-diffusion-inflow ;;actual inflow to target voxels
            outflow * patch-peat-water-diffusion-outflow-target-voxels-fraction

        set patch-peat-water-diffusion-inflow-list lput patch-peat-water-diffusion-inflow patch-peat-water-diffusion-inflow-list

        set patch-peat-water-diffusion-inflow-source-voxels-list lput source patch-peat-water-diffusion-inflow-source-voxels-list
      ]
    ]
  ]

  ;; limit inflow to the maximum capacity of peat retention
  ask voxels-peat
  [
    let unlimited-inflow-list patch-peat-water-diffusion-inflow-list
    let sum-unlimited-inflow sum unlimited-inflow-list
    let sum-maximum-inflow max list 0 (patch-peat-saturated-capacity - patch-water-height)

    if(sum-unlimited-inflow > sum-maximum-inflow)
    [
      let i-fraction (map[? / sum-unlimited-inflow] unlimited-inflow-list)
      let i (map[sum-maximum-inflow * ?] i-fraction)
      set patch-peat-water-diffusion-inflow-delta-list (map[max list 0 (?1 - ?2)] unlimited-inflow-list i)
      set patch-peat-water-diffusion-inflow-list i
    ]
  ]

  ;; delta change outflow by source voxels -- as limited by maximum water retention
  ask voxels-peat
  [
    if(length patch-peat-water-diffusion-inflow-delta-list > 0)
    [
      let ldi patch-peat-water-diffusion-inflow-delta-list
      let lsi patch-peat-water-diffusion-inflow-source-voxels-list

      foreach ldi
      [
        let di ?
        let i position di ldi
        let source item i lsi
        ask source
        [
            set patch-peat-water-diffusion-outflow-delta-list lput di patch-peat-water-diffusion-outflow-delta-list
        ]
      ]

    ]
  ]

  ;; synchronise the outflow from source voxels according to maximum inflow of target voxels
  ask voxels-peat
  [
     set patch-peat-water-diffusion-outflow max list 0 (patch-peat-water-diffusion-outflow - (sum patch-peat-water-diffusion-outflow-delta-list))
  ]

  ;; update water height
  ask voxels-peat
  [
      set patch-water-height
        max list 0
        (patch-water-height + sum patch-peat-water-diffusion-inflow-list - patch-peat-water-diffusion-outflow)
  ]

  ;; convert water level from positive term into negative term
  ask patches
  [
    set patch-water-height-peat -9999
  ]

  ask voxels-peat
  [
    set patch-water-height-peat patch-water-height - patch-peat-saturated-capacity
  ]
end 

to run-simulate
  repeat run-day-number
  [
    dynamic-simulate-surface-water
    dynamic-diffuse-peat-water


    gis-display-voxels
    gis-report-daily

    tick

    if(run-day < run-day-number)
    [
      set run-day run-day + 1
    ]
  ]
end 

to run-warm-up
  repeat warming-up-duration
  [
    dynamic-simulate-surface-water
    dynamic-diffuse-peat-water

    tick

    if(run-day < run-day-number)
    [
      set run-day run-day + 1
    ]
  ]
end 

to run-go
  set message-status "warming up..."
  run-warm-up
  clear-all-plots
  set run-day 1
  reset-ticks

  set message-status "simulating..."
  run-simulate
  set message-status "simulation is done!"
end 

to-report csv-read[file is-lookup? is-multi-string? n-string]
  file-open file
    let first-line file-read-line
    let comma? true
    let col 0

    while [comma?]
    [
      set first-line remove (substring first-line 0 ((position "," first-line) + 1)) first-line
      set comma? member? "," first-line
      set col col + 1
    ]

    let data-list n-values (col + 1) [[]]


    while [not file-at-end?]
    [
       let index 0
       let init file-read-line
       let v substring init 0 (position "," init)

       set data-list ifelse-value (is-lookup?)
       [
         replace-item index data-list (lput v item index data-list)
       ]
       [
         replace-item index data-list (lput (read-from-string v) item index data-list)
       ]

       let r ifelse-value (is-lookup?)
       [
         remove-item (position "," (remove v init)) (remove v init)
       ]
       [
         remove-item (position "," (substring init (position "," init) (length init))) (substring init (position "," init) (length init))
       ]

       ifelse col > 1
       [
         repeat (col - 1)
         [
           set index index + 1
           set v substring r 0 (position "," r)
           set r ifelse-value (is-multi-string?)
           [
             ifelse-value (index < n-string)
             [
               remove-item (position "," (remove v r)) (remove v r)
             ]
             [
               remove-item (position "," (substring r (position "," r) (length r))) (substring r (position "," r) (length r))
             ]
           ]
           [
             remove-item (position "," (substring r (position "," r) (length r))) (substring r (position "," r) (length r))
           ]

           set data-list ifelse-value (is-multi-string?)
           [
             ifelse-value (index < n-string)
             [
               replace-item index data-list (lput v item index data-list)
             ]
             [
               replace-item index data-list (lput (read-from-string v) item index data-list)
             ]
           ]
           [
             replace-item index data-list (lput (read-from-string v) item index data-list)
           ]
         ]
         set data-list replace-item (index + 1) data-list (lput (read-from-string r) item (index + 1) data-list)
       ]
       [
         set index index + 1
         set data-list replace-item index data-list (lput (read-from-string r) item index data-list)
      ]
    ]

  file-close
  report data-list
end 

to csv-read-input
  set lookup-parameters-peat csv-read (word input-path-csv array:item input-array-csv 0) true false 0
  set lookup-parameters-climate csv-read (word input-path-csv array:item input-array-csv 1) false false 0
  set lookup-parameters-land csv-read (word input-path-csv array:item input-array-csv 2) true false 0
  set lookup-parameters-scenario csv-read (word input-path-csv array:item input-array-csv 3) true false 0
end 

to gis-report-daily
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]

  export-view (word output-path-png word "peat-water-day-" word day ".png")

  output-report-csv1 word output-path-csv "output-gwl-average.csv"
  output-report-csv2 word output-path-csv "output-gwl-observation1.csv"
  output-report-csv3 word output-path-csv "output-gwl-observation2.csv"
  output-report-csv4 word output-path-csv "output-gwl-observation3.csv"
  output-report-csv5 word output-path-csv "output-gwl-observation4.csv"
  output-report-csv6 word output-path-csv "output-gwl-observation5.csv"

  output-report-csv7 word output-path-csv "output-wc-average.csv"
  output-report-csv8 word output-path-csv "output-wc-observation1.csv"
  output-report-csv9 word output-path-csv "output-wc-observation2.csv"
  output-report-csv10 word output-path-csv "output-wc-observation3.csv"
  output-report-csv11 word output-path-csv "output-wc-observation4.csv"
  output-report-csv12 word output-path-csv "output-wc-observation5.csv"

  output-report-csv13 word output-path-csv "output-wc-landcover-forest.csv"
  output-report-csv14 word output-path-csv "output-wc-landcover-shrub.csv"
  output-report-csv15 word output-path-csv "output-wc-landcover-burntland.csv"
  output-report-csv16 word output-path-csv "output-wc-landcover-bareland.csv"

  output-report-csv17 word output-path-csv "output-gwl-landcover-forest.csv"
  output-report-csv18 word output-path-csv "output-gwl-landcover-shrub.csv"
  output-report-csv19 word output-path-csv "output-gwl-landcover-burntland.csv"
  output-report-csv20 word output-path-csv "output-gwl-landcover-bareland.csv"

  if(report-ascii-map?)
  [
    let l sort remove-duplicates[patch-depth]of voxels-peat ;;m

    foreach l
    [
      let d-m ?
      let d-cm d-m * 100
      let p voxels-peat with[patch-depth = d-m]

      ask patches with[pzcor = 0]
      [
        set patch-water-height-ascii -9999
      ]

      let z max[pzcor]of voxels-peat with[patch-depth = d-m]

      ask patches with[pzcor = z]
      [
        let x pxcor
        let y pycor
        let w patch-water-height
        ask patches with[pzcor = 0]
        [
          if(pxcor = x and pycor = y)
          [
            set patch-water-height-ascii w
          ]
        ]
      ]

      let m gis:patch-dataset patch-water-height-ascii
      let f word output-path-asc "tmp.asc"
      gis:store-dataset m f

      gis-replace-header-asc ("tmp.asc") (word "peat-water-day-" word day word "-depth-" d-cm "-cm.asc")
    ]
  ]
end 

to gis-replace-header-asc[f1 f2]
  let d1 word output-path-asc f1

  let num 0
  let line ""
  let line-data []

  file-open d1
     while [not file-at-end?]
     [
       set line file-read-line
       set line-data lput line line-data
       set num num + 1
     ]
  file-close

  file-delete d1

  let d2 word output-path-asc f2

  if (file-exists? d2)
  [
    file-delete d2
  ]

  file-open d2
    file-print word "ncols         " map-attributes-n-cols
    file-print word "nrows         " map-attributes-n-rows
    file-print word "xllcorner     " map-attributes-corner-xll
    file-print word "yllcorner     " map-attributes-corner-yll
    file-print word "cellsize      " map-attributes-cell-size-degree
    file-print "NODATA_value  -9999"
    foreach n-values (num - 6)[?]
    [
      file-print item (? + 6) line-data
    ]
  file-close
end 

to update-scenario-number
  set parameter-scenario-number parameter-scenario-number + 1
  let f (word input-path-csv array:item input-array-csv 3)
  file-delete f

  file-open f
      file-type "Scenario parameters"
      file-type ","
      file-print "Value"
      file-type "Scenario number"
      file-type ","
      file-print parameter-scenario-number
  file-close
end 

to output-write-main-header-csv1[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat) ;;cm
  let h["Day"]
  set h lput "Average surface water (mm/pixel)" h
  set h lput "Average acrotelm peat water (mm/pixel)" h
  set h lput "Average catotelm peat water (mm/pixel)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm/pixel)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Pixel average"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv1[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat ;;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm ;;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm ;;m

  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface
     file-type comma

     let a (map[mean[patch-water-height-peat]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height-peat]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height-peat]of voxels-peat with[patch-depth = ?]] l)

     let s[]
     foreach l
     [
       let i position ? l
       set s lput (sum sublist m 0 (i + 1)) s
     ]

     foreach bl s
     [
       file-type ?
       file-type comma
     ]

     file-print last s
  file-close
end 

to output-write-main-header-csv2[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 1]) ;;cm
  let h["Day"]
  set h lput "Surface water observation 1 (mm)" h
  set h lput "Acrotelm peat water observation 1 (mm)" h
  set h lput "Catotelm peat water observation 1 (mm)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Observation point #1"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv2[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 1];;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm with[patch-observation-id = 1];;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm with[patch-observation-id = 1];;m


  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface with[patch-observation-id = 1]
     file-type comma

     let a (map[mean[patch-water-height-peat]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height-peat]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height-peat]of voxels-peat with[patch-depth = ? and patch-observation-id = 1]] l)

     let s[]
     foreach l
     [
       let i position ? l
       set s lput (sum sublist m 0 (i + 1)) s
     ]

     foreach bl s
     [
       file-type ?
       file-type comma
     ]

     file-print last s
  file-close
end 

to output-write-main-header-csv3[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 2]) ;;cm
  let h["Day"]
  set h lput "Surface water observation 2 (mm)" h
  set h lput "Acrotelm peat water observation 2 (mm)" h
  set h lput "Catotelm peat water observation 2 (mm)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Observation point #2"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv3[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 2];;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm with[patch-observation-id = 2];;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm with[patch-observation-id = 2];;m


  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface with[patch-observation-id = 2]
     file-type comma

     let a (map[mean[patch-water-height-peat]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height-peat]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height-peat]of voxels-peat with[patch-depth = ? and patch-observation-id = 2]] l)

     let s[]
     foreach l
     [
       let i position ? l
       set s lput (sum sublist m 0 (i + 1)) s
     ]

     foreach bl s
     [
       file-type ?
       file-type comma
     ]

     file-print last s

  file-close
end 

to output-write-main-header-csv4[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 3]) ;;cm
  let h["Day"]
  set h lput "Surface water observation 3 (mm)" h
  set h lput "Acrotelm peat water observation 3 (mm)" h
  set h lput "Catotelm peat water observation 3 (mm)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Observation point #3"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv4[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 3];;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm with[patch-observation-id = 3];;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm with[patch-observation-id = 3];;m


  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface with[patch-observation-id = 3]
     file-type comma

     let a (map[mean[patch-water-height-peat]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height-peat]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height-peat]of voxels-peat with[patch-depth = ? and patch-observation-id = 3]] l)

     let s[]
     foreach l
     [
       let i position ? l
       set s lput (sum sublist m 0 (i + 1)) s
     ]

     foreach bl s
     [
       file-type ?
       file-type comma
     ]

     file-print last s

  file-close
end 

to output-write-main-header-csv5[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 4]) ;;cm
  let h["Day"]
  set h lput "Surface water observation 4 (mm)" h
  set h lput "Acrotelm peat water observation 4 (mm)" h
  set h lput "Catotelm peat water observation 4 (mm)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Observation point #4"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv5[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 4];;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm with[patch-observation-id = 4];;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm with[patch-observation-id = 4];;m


  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface with[patch-observation-id = 4]
     file-type comma

     let a (map[mean[patch-water-height-peat]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height-peat]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height-peat]of voxels-peat with[patch-depth = ? and patch-observation-id = 4]] l)

     let s[]
     foreach l
     [
       let i position ? l
       set s lput (sum sublist m 0 (i + 1)) s
     ]

     foreach bl s
     [
       file-type ?
       file-type comma
     ]

     file-print last s

  file-close
end 

to output-write-main-header-csv6[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 5]) ;;cm
  let h["Day"]
  set h lput "Surface water observation 5 (mm)" h
  set h lput "Acrotelm peat water observation 5 (mm)" h
  set h lput "Catotelm peat water observation 5 (mm)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Observation point #5"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv6[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 5];;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm with[patch-observation-id = 5];;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm with[patch-observation-id = 5];;m

  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface with[patch-observation-id = 5]
     file-type comma

     let a (map[mean[patch-water-height-peat]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height-peat]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height-peat]of voxels-peat with[patch-depth = ? and patch-observation-id = 5]] l)

     let s[]
     foreach l
     [
       let i position ? l
       set s lput (sum sublist m 0 (i + 1)) s
     ]

     foreach bl s
     [
       file-type ?
       file-type comma
     ]

     file-print last s

  file-close
end 

;;add reporting

to output-write-main-header-csv7[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat) ;;cm
  let h["Day"]
  set h lput "Average surface water (mm/pixel)" h
  set h lput "Average acrotelm peat water (mm/pixel)" h
  set h lput "Average catotelm peat water (mm/pixel)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm/pixel)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Pixel average"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv7[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat ;;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm ;;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm ;;m

  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface
     file-type comma

     let a (map[mean[patch-water-height]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height]of voxels-peat with[patch-depth = ?]] l)

     foreach bl m
     [
       file-type ?
       file-type comma
     ]

     file-print last m
  file-close
end 

to output-write-main-header-csv8[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 1]) ;;cm
  let h["Day"]
  set h lput "Surface water observation 1 (mm)" h
  set h lput "Acrotelm peat water observation 1 (mm)" h
  set h lput "Catotelm peat water observation 1 (mm)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Observation point #1"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv8[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 1];;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm with[patch-observation-id = 1];;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm with[patch-observation-id = 1];;m

  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface with[patch-observation-id = 1]
     file-type comma

     let a (map[mean[patch-water-height]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height]of voxels-peat with[patch-depth = ? and patch-observation-id = 1]] l)

     foreach bl m
     [
       file-type ?
       file-type comma
     ]

     file-print last m
  file-close
end 

to output-write-main-header-csv9[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 2]) ;;cm
  let h["Day"]
  set h lput "Surface water observation 2 (mm)" h
  set h lput "Acrotelm peat water observation 2 (mm)" h
  set h lput "Catotelm peat water observation 2 (mm)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Observation point #2"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv9[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 2];;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm with[patch-observation-id = 2];;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm with[patch-observation-id = 2];;m

  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface with[patch-observation-id = 2]
     file-type comma

     let a (map[mean[patch-water-height]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height]of voxels-peat with[patch-depth = ? and patch-observation-id = 2]] l)

     foreach bl m
     [
       file-type ?
       file-type comma
     ]

     file-print last m

  file-close
end 

to output-write-main-header-csv10[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 3]) ;;cm
  let h["Day"]
  set h lput "Surface water observation 3 (mm)" h
  set h lput "Acrotelm peat water observation 3 (mm)" h
  set h lput "Catotelm peat water observation 3 (mm)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Observation point #3"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv10[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 3];;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm with[patch-observation-id = 3];;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm with[patch-observation-id = 3];;m

  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface with[patch-observation-id = 3]
     file-type comma

     let a (map[mean[patch-water-height]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height]of voxels-peat with[patch-depth = ? and patch-observation-id = 3]] l)

     foreach bl m
     [
       file-type ?
       file-type comma
     ]

     file-print last m

  file-close
end 

to output-write-main-header-csv11[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 4]) ;;cm
  let h["Day"]
  set h lput "Surface water observation 4 (mm)" h
  set h lput "Acrotelm peat water observation 4 (mm)" h
  set h lput "Catotelm peat water observation 4 (mm)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Observation point #4"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv11[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 4];;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm with[patch-observation-id = 4];;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm with[patch-observation-id = 4];;m

  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface with[patch-observation-id = 4]
     file-type comma

     let a (map[mean[patch-water-height]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height]of voxels-peat with[patch-depth = ? and patch-observation-id = 4]] l)

     foreach bl m
     [
       file-type ?
       file-type comma
     ]

     file-print last m

  file-close
end 

to output-write-main-header-csv12[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 5]) ;;cm
  let h["Day"]
  set h lput "Surface water observation 5 (mm)" h
  set h lput "Acrotelm peat water observation 5 (mm)" h
  set h lput "Catotelm peat water observation 5 (mm)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Observation point #5"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv12[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat with[patch-observation-id = 5];;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm with[patch-observation-id = 5];;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm with[patch-observation-id = 5];;m

  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface with[patch-observation-id = 5]
     file-type comma

     let a (map[mean[patch-water-height]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height]of voxels-peat with[patch-depth = ? and patch-observation-id = 5]] l)

     foreach bl m
     [
       file-type ?
       file-type comma
     ]

     file-print last m

  file-close
end 

;; add report by land use

to output-write-main-header-csv13[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat) ;;cm
  let h["Day"]
  set h lput "Average surface water (mm/pixel)" h
  set h lput "Average acrotelm peat water (mm/pixel)" h
  set h lput "Average catotelm peat water (mm/pixel)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm/pixel)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Pixel average - land use/cover type: forest"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv13[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat with[patch-land = 1];;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm with[patch-land = 1];;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm with[patch-land = 1];;m


  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface with[patch-land = 1]
     file-type comma

     let a (map[mean[patch-water-height]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height]of voxels-peat with[patch-depth = ? and patch-land = 1]] l)

     foreach bl m
     [
       file-type ?
       file-type comma
     ]

     file-print last m
  file-close
end 

to output-write-main-header-csv14[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat) ;;cm
  let h["Day"]
  set h lput "Average surface water (mm/pixel)" h
  set h lput "Average acrotelm peat water (mm/pixel)" h
  set h lput "Average catotelm peat water (mm/pixel)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm/pixel)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Pixel average - land use/cover type: shrub"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv14[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat with[patch-land = 2];;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm with[patch-land = 2];;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm with[patch-land = 2];;m


  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface with[patch-land = 2]
     file-type comma

     let a (map[mean[patch-water-height]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height]of voxels-peat with[patch-depth = ? and patch-land = 2]] l)

     foreach bl m
     [
       file-type ?
       file-type comma
     ]

     file-print last m
  file-close
end 

to output-write-main-header-csv15[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat) ;;cm
  let h["Day"]
  set h lput "Average surface water (mm/pixel)" h
  set h lput "Average acrotelm peat water (mm/pixel)" h
  set h lput "Average catotelm peat water (mm/pixel)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm/pixel)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Pixel average - land use/cover type: burnt land"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv15[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat with[patch-land = 3];;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm with[patch-land = 3];;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm with[patch-land = 3];;m


  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface with[patch-land = 3]
     file-type comma

     let a (map[mean[patch-water-height]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height]of voxels-peat with[patch-depth = ? and patch-land = 3]] l)

     foreach bl m
     [
       file-type ?
       file-type comma
     ]

     file-print last m
  file-close
end 

to output-write-main-header-csv16[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat) ;;cm
  let h["Day"]
  set h lput "Average surface water (mm/pixel)" h
  set h lput "Average acrotelm peat water (mm/pixel)" h
  set h lput "Average catotelm peat water (mm/pixel)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm/pixel)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Pixel average - land use/cover type: bareland"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv16[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat with[patch-land = 4];;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm with[patch-land = 4];;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm with[patch-land = 4];;m


  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface with[patch-land = 4]
     file-type comma

     let a (map[mean[patch-water-height]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height]of voxels-peat with[patch-depth = ? and patch-land = 4]] l)

     foreach bl m
     [
       file-type ?
       file-type comma
     ]

     file-print last m
  file-close
end 

;; add GWL by land use

to output-write-main-header-csv17[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat) ;;cm
  let h["Day"]
  set h lput "Average surface water (mm/pixel)" h
  set h lput "Average acrotelm peat water (mm/pixel)" h
  set h lput "Average catotelm peat water (mm/pixel)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm/pixel)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Pixel average - land use/cover type: forest"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv17[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat with[patch-land = 1];;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm with[patch-land = 1];;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm with[patch-land = 1];;m


  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface with[patch-land = 1]
     file-type comma

     let a (map[mean[patch-water-height-peat]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height-peat]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height-peat]of voxels-peat with[patch-depth = ? and patch-land = 1]] l)

     foreach bl m
     [
       file-type ?
       file-type comma
     ]

     file-print last m
  file-close
end 

to output-write-main-header-csv18[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat) ;;cm
  let h["Day"]
  set h lput "Average surface water (mm/pixel)" h
  set h lput "Average acrotelm peat water (mm/pixel)" h
  set h lput "Average catotelm peat water (mm/pixel)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm/pixel)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Pixel average - land use/cover type: shrub"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv18[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat with[patch-land = 2];;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm with[patch-land = 2];;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm with[patch-land = 2];;m


  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface with[patch-land = 2]
     file-type comma

     let a (map[mean[patch-water-height-peat]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height-peat]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height-peat]of voxels-peat with[patch-depth = ? and patch-land = 2]] l)

     foreach bl m
     [
       file-type ?
       file-type comma
     ]

     file-print last m
  file-close
end 

to output-write-main-header-csv19[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat) ;;cm
  let h["Day"]
  set h lput "Average surface water (mm/pixel)" h
  set h lput "Average acrotelm peat water (mm/pixel)" h
  set h lput "Average catotelm peat water (mm/pixel)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm/pixel)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Pixel average - land use/cover type: burnt land"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv19[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat with[patch-land = 3];;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm with[patch-land = 3];;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm with[patch-land = 3];;m


  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface with[patch-land = 3]
     file-type comma

     let a (map[mean[patch-water-height-peat]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height-peat]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height-peat]of voxels-peat with[patch-depth = ? and patch-land = 3]] l)

     foreach bl m
     [
       file-type ?
       file-type comma
     ]

     file-print last m
  file-close
end 

to output-write-main-header-csv20[f]
  let l (map[? * 100] sort remove-duplicates[patch-depth]of voxels-peat) ;;cm
  let h["Day"]
  set h lput "Average surface water (mm/pixel)" h
  set h lput "Average acrotelm peat water (mm/pixel)" h
  set h lput "Average catotelm peat water (mm/pixel)" h

  foreach l
  [
       set h lput (word "Peat water layer " word ? "-cm depth (mm/pixel)") h
  ]
  let comma ","
  file-open f
     file-print "P E A T   V O X E L   M O D E L"
     file-print ""
     file-print "Developed by: Aryo Adhi Condro - 2017"
     file-print ""
     file-print "Supervisors: Idung Risdiyanto S. Si M.Sc. & Prof. Dr. Hidayat Pawitan"
     file-print ""
     file-print "GFM - IPB - Bogor - Indonesia"
     file-print ""
     file-print word "Simulation results reported on: " date-and-time
     file-print ""
     file-print "Pixel average - land use/cover type: bareland"
     file-print ""

     foreach bl h
     [
       file-type ?
       file-type comma
     ]

     file-print last h

  file-close
end 

to output-report-csv20[f]
  let day ifelse-value(run-day < 10) [word "0" run-day][run-day]
  let comma ","
  let l sort remove-duplicates[patch-depth]of voxels-peat with[patch-land = 4];;m
  let al sort remove-duplicates[patch-depth]of voxels-acrotelm with[patch-land = 4];;m
  let cl sort remove-duplicates[patch-depth]of voxels-catotelm with[patch-land = 4];;m


  file-open f
     file-type day
     file-type comma
     file-type mean[patch-water-height]of voxels-surface with[patch-land = 4]
     file-type comma

     let a (map[mean[patch-water-height-peat]of voxels-acrotelm with[patch-depth = ?]] al)
     let c (map[mean[patch-water-height-peat]of voxels-catotelm with[patch-depth = ?]] cl)

     file-type sum a
     file-type comma
     file-type sum a + sum c
     file-type comma

     let m (map[mean[patch-water-height-peat]of voxels-peat with[patch-depth = ? and patch-land = 4]] l)

     foreach bl m
     [
       file-type ?
       file-type comma
     ]

     file-print last m
  file-close
end