diff --git a/._ef5.ico b/._ef5.ico new file mode 100644 index 0000000..4c0326b Binary files /dev/null and b/._ef5.ico differ diff --git a/.gitignore b/.gitignore index 278c536..17884ed 100755 --- a/.gitignore +++ b/.gitignore @@ -21,3 +21,4 @@ depcomp stamp-h1 *.cache/* src/.* +EF5-US-Parameters-1.0.0/* diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..c320594 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,68 @@ +{ + "files.associations": { + "*.tbl": "cpp", + "defines.h": "c", + "configsection.h": "c", + "model.h": "c", + "modelbase.h": "c", + "basicgrids.h": "c", + "cctype": "c", + "clocale": "c", + "cmath": "c", + "cstdarg": "c", + "cstddef": "c", + "cstdio": "c", + "cstdlib": "c", + "cstring": "c", + "ctime": "c", + "cwchar": "c", + "cwctype": "c", + "array": "c", + "atomic": "c", + "strstream": "c", + "bit": "c", + "*.tcc": "c", + "bitset": "c", + "chrono": "c", + "condition_variable": "c", + "cstdint": "c", + "deque": "c", + "list": "c", + "map": "c", + "set": "c", + "unordered_map": "c", + "vector": "c", + "exception": "c", + "algorithm": "c", + "functional": "c", + "iterator": "c", + "memory": "c", + "memory_resource": "c", + "numeric": "c", + "optional": "c", + "random": "c", + "ratio": "c", + "string": "c", + "string_view": "c", + "system_error": "c", + "tuple": "c", + "type_traits": "c", + "utility": "c", + "fstream": "c", + "initializer_list": "c", + "iosfwd": "c", + "iostream": "c", + "istream": "c", + "limits": "c", + "mutex": "c", + "new": "c", + "ostream": "c", + "shared_mutex": "c", + "sstream": "c", + "stdexcept": "c", + "streambuf": "c", + "thread": "c", + "cinttypes": "c", + "typeinfo": "c" + } +} \ No newline at end of file diff --git a/Makefile.am b/Makefile.am index 6afc2c3..d19fc7f 100755 --- a/Makefile.am +++ b/Makefile.am @@ -9,7 +9,7 @@ unit_FILES = src/LAEAProjection.cpp src/GeographicProjection.cpp src/DistanceUni type_FILES = src/DatedName.cpp src/PETType.cpp src/PrecipType.cpp src/TempType.cpp src/GaugeMap.cpp config_FILES = src/BasicConfigSection.cpp src/PrecipConfigSection.cpp src/PETConfigSection.cpp src/TempConfigSection.cpp src/GaugeConfigSection.cpp src/BasinConfigSection.cpp src/CaliParamConfigSection.cpp src/ParamSetConfigSection.cpp src/RoutingCaliParamConfigSection.cpp src/RoutingParamSetConfigSection.cpp src/TaskConfigSection.cpp src/EnsTaskConfigSection.cpp src/ExecuteConfigSection.cpp src/Config.cpp src/SnowCaliParamConfigSection.cpp src/SnowParamSetConfigSection.cpp src/InundationCaliParamConfigSection.cpp src/InundationParamSetConfigSection.cpp input_FILES = src/RPSkewness.cpp src/TimeSeries.cpp src/PETReader.cpp src/PrecipReader.cpp src/TempReader.cpp src/TifGrid.cpp src/BifGrid.cpp src/AscGrid.cpp src/BasicGrids.cpp src/TRMMRTGrid.cpp src/MRMSGrid.cpp src/GridWriter.cpp src/GridWriterFull.cpp src/GriddedOutput.cpp -model_FILES = src/Model.cpp src/CRESTModel.cpp src/HyMOD.cpp src/SAC.cpp src/LinearRoute.cpp src/KinematicRoute.cpp src/ObjectiveFunc.cpp src/Simulator.cpp src/ARS.cpp src/DREAM.cpp src/dream_functions.cpp src/misc_functions.cpp src/Snow17Model.cpp src/HPModel.cpp src/SimpleInundation.cpp src/VCInundation.cpp +model_FILES = src/Model.cpp src/CRESTModel.cpp src/CRESTPhysModel.cpp src/HyMOD.cpp src/SAC.cpp src/LinearRoute.cpp src/KinematicRoute.cpp src/ObjectiveFunc.cpp src/Simulator.cpp src/ARS.cpp src/DREAM.cpp src/dream_functions.cpp src/misc_functions.cpp src/Snow17Model.cpp src/HPModel.cpp src/SimpleInundation.cpp src/VCInundation.cpp if WINDOWS AM_CXXFLAGS= -Wall -mwindows ${OPENMP_CFLAGS} __top_builddir__bin_ef5_SOURCES = $(unit_FILES) $(type_FILES) $(config_FILES) $(input_FILES) $(model_FILES) src/ExecutionController.cpp src/EF5Windows.cpp src/ef5.rc diff --git a/README.md b/README.md index 88b8a77..545338c 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ Ensemble Framework For Flash Flood Forecasting (EF5) === -![version](https://img.shields.io/badge/version-1.1-blue.svg?style=flat) [![Build Status](https://travis-ci.org/HyDROSLab/EF5.svg?branch=master)](https://travis-ci.org/HyDROSLab/EF5) +![version](https://img.shields.io/badge/version-v1.3.0-blue) [![Build Status](https://travis-ci.org/HyDROSLab/EF5.svg?branch=master)](https://travis-ci.org/HyDROSLab/EF5) ![lastCommit](https://img.shields.io/github/last-commit/chrimerss/EF5) EF5 was created by the Hydrometeorology and Remote Sensing Laboratory at the University of Oklahoma. The goal of EF5 is to have a framework for distributed hydrologic modeling that is user friendly, adaptable, expandable, all while being suitable for large scale (e.g. continental scale) modeling of flash floods with rapid forecast updates. Currently EF5 incorporates 3 water balance models including the Sacramento Soil Moisture Accouning Model (SAC-SMA), Coupled Routing and Excess Storage (CREST), and hydrophobic (HP). These water balance models can be coupled with either linear reservoir or kinematic wave routing. @@ -9,7 +9,7 @@ The goal of EF5 is to have a framework for distributed hydrologic modeling that EF5 has a homepage at [http://ef5.ou.edu](http://ef5.ou.edu). The training modules are found at [http://ef5.ou.edu/training/](http://ef5.ou.edu/training/) while the YouTube videos may be found at [https://www.youtube.com/channel/UCgoGJtdeqHgwoYIRhkgMwog](https://www.youtube.com/channel/UCgoGJtdeqHgwoYIRhkgMwog). The source code is found on GitHub at [https://github.com/HyDROSLab/EF5](https://github.com/HyDROSLab/EF5). -See [manual.html](manual.html) for the EF5 operating manual which describes configuration options. +See [manual.html](https://chrimerss.github.io/EF5/docs/) for the EF5 operating manual which describes configuration options. ## Compiling @@ -52,3 +52,4 @@ JJ Gourley Yang Hong +Zhi Li diff --git a/compile b/compile new file mode 100755 index 0000000..99e5052 --- /dev/null +++ b/compile @@ -0,0 +1,348 @@ +#! /bin/sh +# Wrapper for compilers which do not understand '-c -o'. + +scriptversion=2018-03-07.03; # UTC + +# Copyright (C) 1999-2018 Free Software Foundation, Inc. +# Written by Tom Tromey . +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2, or (at your option) +# any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +# As a special exception to the GNU General Public License, if you +# distribute this file as part of a program that contains a +# configuration script generated by Autoconf, you may include it under +# the same distribution terms that you use for the rest of that program. + +# This file is maintained in Automake, please report +# bugs to or send patches to +# . + +nl=' +' + +# We need space, tab and new line, in precisely that order. Quoting is +# there to prevent tools from complaining about whitespace usage. +IFS=" "" $nl" + +file_conv= + +# func_file_conv build_file lazy +# Convert a $build file to $host form and store it in $file +# Currently only supports Windows hosts. If the determined conversion +# type is listed in (the comma separated) LAZY, no conversion will +# take place. +func_file_conv () +{ + file=$1 + case $file in + / | /[!/]*) # absolute file, and not a UNC file + if test -z "$file_conv"; then + # lazily determine how to convert abs files + case `uname -s` in + MINGW*) + file_conv=mingw + ;; + CYGWIN*) + file_conv=cygwin + ;; + *) + file_conv=wine + ;; + esac + fi + case $file_conv/,$2, in + *,$file_conv,*) + ;; + mingw/*) + file=`cmd //C echo "$file " | sed -e 's/"\(.*\) " *$/\1/'` + ;; + cygwin/*) + file=`cygpath -m "$file" || echo "$file"` + ;; + wine/*) + file=`winepath -w "$file" || echo "$file"` + ;; + esac + ;; + esac +} + +# func_cl_dashL linkdir +# Make cl look for libraries in LINKDIR +func_cl_dashL () +{ + func_file_conv "$1" + if test -z "$lib_path"; then + lib_path=$file + else + lib_path="$lib_path;$file" + fi + linker_opts="$linker_opts -LIBPATH:$file" +} + +# func_cl_dashl library +# Do a library search-path lookup for cl +func_cl_dashl () +{ + lib=$1 + found=no + save_IFS=$IFS + IFS=';' + for dir in $lib_path $LIB + do + IFS=$save_IFS + if $shared && test -f "$dir/$lib.dll.lib"; then + found=yes + lib=$dir/$lib.dll.lib + break + fi + if test -f "$dir/$lib.lib"; then + found=yes + lib=$dir/$lib.lib + break + fi + if test -f "$dir/lib$lib.a"; then + found=yes + lib=$dir/lib$lib.a + break + fi + done + IFS=$save_IFS + + if test "$found" != yes; then + lib=$lib.lib + fi +} + +# func_cl_wrapper cl arg... +# Adjust compile command to suit cl +func_cl_wrapper () +{ + # Assume a capable shell + lib_path= + shared=: + linker_opts= + for arg + do + if test -n "$eat"; then + eat= + else + case $1 in + -o) + # configure might choose to run compile as 'compile cc -o foo foo.c'. + eat=1 + case $2 in + *.o | *.[oO][bB][jJ]) + func_file_conv "$2" + set x "$@" -Fo"$file" + shift + ;; + *) + func_file_conv "$2" + set x "$@" -Fe"$file" + shift + ;; + esac + ;; + -I) + eat=1 + func_file_conv "$2" mingw + set x "$@" -I"$file" + shift + ;; + -I*) + func_file_conv "${1#-I}" mingw + set x "$@" -I"$file" + shift + ;; + -l) + eat=1 + func_cl_dashl "$2" + set x "$@" "$lib" + shift + ;; + -l*) + func_cl_dashl "${1#-l}" + set x "$@" "$lib" + shift + ;; + -L) + eat=1 + func_cl_dashL "$2" + ;; + -L*) + func_cl_dashL "${1#-L}" + ;; + -static) + shared=false + ;; + -Wl,*) + arg=${1#-Wl,} + save_ifs="$IFS"; IFS=',' + for flag in $arg; do + IFS="$save_ifs" + linker_opts="$linker_opts $flag" + done + IFS="$save_ifs" + ;; + -Xlinker) + eat=1 + linker_opts="$linker_opts $2" + ;; + -*) + set x "$@" "$1" + shift + ;; + *.cc | *.CC | *.cxx | *.CXX | *.[cC]++) + func_file_conv "$1" + set x "$@" -Tp"$file" + shift + ;; + *.c | *.cpp | *.CPP | *.lib | *.LIB | *.Lib | *.OBJ | *.obj | *.[oO]) + func_file_conv "$1" mingw + set x "$@" "$file" + shift + ;; + *) + set x "$@" "$1" + shift + ;; + esac + fi + shift + done + if test -n "$linker_opts"; then + linker_opts="-link$linker_opts" + fi + exec "$@" $linker_opts + exit 1 +} + +eat= + +case $1 in + '') + echo "$0: No command. Try '$0 --help' for more information." 1>&2 + exit 1; + ;; + -h | --h*) + cat <<\EOF +Usage: compile [--help] [--version] PROGRAM [ARGS] + +Wrapper for compilers which do not understand '-c -o'. +Remove '-o dest.o' from ARGS, run PROGRAM with the remaining +arguments, and rename the output as expected. + +If you are trying to build a whole package this is not the +right script to run: please start by reading the file 'INSTALL'. + +Report bugs to . +EOF + exit $? + ;; + -v | --v*) + echo "compile $scriptversion" + exit $? + ;; + cl | *[/\\]cl | cl.exe | *[/\\]cl.exe | \ + icl | *[/\\]icl | icl.exe | *[/\\]icl.exe ) + func_cl_wrapper "$@" # Doesn't return... + ;; +esac + +ofile= +cfile= + +for arg +do + if test -n "$eat"; then + eat= + else + case $1 in + -o) + # configure might choose to run compile as 'compile cc -o foo foo.c'. + # So we strip '-o arg' only if arg is an object. + eat=1 + case $2 in + *.o | *.obj) + ofile=$2 + ;; + *) + set x "$@" -o "$2" + shift + ;; + esac + ;; + *.c) + cfile=$1 + set x "$@" "$1" + shift + ;; + *) + set x "$@" "$1" + shift + ;; + esac + fi + shift +done + +if test -z "$ofile" || test -z "$cfile"; then + # If no '-o' option was seen then we might have been invoked from a + # pattern rule where we don't need one. That is ok -- this is a + # normal compilation that the losing compiler can handle. If no + # '.c' file was seen then we are probably linking. That is also + # ok. + exec "$@" +fi + +# Name of file we expect compiler to create. +cofile=`echo "$cfile" | sed 's|^.*[\\/]||; s|^[a-zA-Z]:||; s/\.c$/.o/'` + +# Create the lock directory. +# Note: use '[/\\:.-]' here to ensure that we don't use the same name +# that we are using for the .o file. Also, base the name on the expected +# object file name, since that is what matters with a parallel build. +lockdir=`echo "$cofile" | sed -e 's|[/\\:.-]|_|g'`.d +while true; do + if mkdir "$lockdir" >/dev/null 2>&1; then + break + fi + sleep 1 +done +# FIXME: race condition here if user kills between mkdir and trap. +trap "rmdir '$lockdir'; exit 1" 1 2 15 + +# Run the compile. +"$@" +ret=$? + +if test -f "$cofile"; then + test "$cofile" = "$ofile" || mv "$cofile" "$ofile" +elif test -f "${cofile}bj"; then + test "${cofile}bj" = "$ofile" || mv "${cofile}bj" "$ofile" +fi + +rmdir "$lockdir" +exit $ret + +# Local Variables: +# mode: shell-script +# sh-indentation: 2 +# eval: (add-hook 'before-save-hook 'time-stamp) +# time-stamp-start: "scriptversion=" +# time-stamp-format: "%:y-%02m-%02d.%02H" +# time-stamp-time-zone: "UTC0" +# time-stamp-end: "; # UTC" +# End: diff --git a/manual.html b/docs/index.html similarity index 89% rename from manual.html rename to docs/index.html index 3814b35..a36b14a 100755 --- a/manual.html +++ b/docs/index.html @@ -1,572 +1,624 @@ - - - EF5 User Manual - - - - -
-

Table of Contents

-
    -
  1. About EF5
    2. -
  2. Hydrologic Water Balance Models -
      -
    1. CREST
      2. -
    2. SAC-SMA
      3. -
    3. HP
      4. -
    -
    3. -
  3. Routing Models -
      -
    1. Linear Reservior
      2. -
    2. Kinematic Wave
      3. -
    -
    4. -
  4. Snow Melt Models -
      -
    1. Snow-17
      2. -
    -
    5. -
  5. Inundation Models -
      -
    1. Simple Inundation
      2. -
    -
    6. -
  6. Compiling EF5
    7. -
  7. Configuration File -
      -
    1. Basic Information
      2. -
    2. Precipitation Information
      3. -
    3. Potentional Evapotranspiration Information
      4. -
    4. Gauge Locations
      5. -
    5. Basins
      6. -
    6. Parameter Sets -
        -
      1. CREST
        2. -
      2. SAC-SMA
        3. -
      3. HP
        4. -
      4. Linear Reservoir
        5. -
      5. Kinematic Wave
        6. -
      6. Snow-17
        7. -
      7. Simple Inundation
        8. -
      -
      7. -
    7. Tasks
      8. -
    8. Execute Block
      9. -
    -
    8. -
  8. Running EF5
    9. -
  9. Calibrating the Models
    10. -
  10. Appendix -
      -
    1. Complete Sample Configuration File
      2. -
    -
    11. -
-
-
-
    -
  1. About EF5
    2. -

    EF5 is designed to facilitate creation of ensemble forecasts for flash flood prediction. As such it will incorporate multi-model support while maintaining a single set of input data. Currently the only supported model is the Coupled Routing and Excess Storage (CREST) hydrologic model. Additionally, EF5 was designed with utilization of parallel computing in mind. Presently portions of CREST are optimized to take advantage of multi-core computing through the use of OpenMP.

    -
  2. Hydrologic Water Balance Models
    3. -
      -
    1. CREST
      2. -

      The Coupled Routing and Excess STorage (CREST) distributed hydrological model is a hybrid modeling strategy that was recently developed by the University of Oklahoma (http://hydro.ou.edu) and NASA SERVIR Project Team (www.servir.net). CREST simulates the spatiotemporal variation of water and energy fluxes and storages on a regular grid with the grid cell resolution being user-defined, thereby enabling global- and regional-scale applications. The scalability of CREST simulations is accomplished through sub-grid scale representation of soil moisture storage capacity (using a variable infiltration curve) and runoff generation processes (using linear reservoirs). The CREST model was initially developed to provide online global flood predictions with relatively coarse resolution, but it is also applicable at small scales, such as single basins. The CREST Model can be forced by gridded potential evapotranspiration and precipitation datasets such as, satellite-based precipitation estimates, gridded rain gauge observations, remote sensing platforms such as weather radar, and quantitative precipitation forecasts from numerical weather prediction models. The representation of the primary water fluxes such as infiltration and routing are closely related to the spatially variable land surface characteristics (i.e., vegetation, soil type, and topography). The runoff generation component and routing scheme are coupled, thus providing realistic interactions between atmospheric, land surface, and subsurface water. -

      More detailed information about CREST can be found in the following publication:
      Wang, J., Y. Hong, L. Li, J. J. Gourley, S. I. Khan, K. K. Yilmaz, R. F. Adler, F. S. Policelli, S. Habib, D. Irwin, A. S. Limaye, T. Korme, and L. Okello, 2011: The coupled routing and excess storage (CREST) distributed hydrological model. Hydrol. Sci. Journal, 56, 84-98, doi: 10.1080/02626667.2010.543087. -

      -
    2. SAC-SMA
      3. -

      The Sacramento Soil Moisture Accounting (SAC-SMA) Model was developed by the U.S. National Weather Service with the goal of parameterizing soil moisture characteristics in a fashion that would 1) logically distribute applied moisture in various depths and energy states in the soil 2) have rational percolation characteristics 3) allow an effective simulation of streamflow. -

      More detailed information about SAC-SMA can be found online at U.S. NWS SAC-SMA algorithm description

      -
    3. HP
      4. -

      The Hydrophobic (HP) water balance model features an entirely impervious surface where all rainfall is transformed into surface runoff.

      -
    -
  3. Routing Models
    4. -
      -
    1. Linear Reservoir
      2. -

      - The linear reservoir routing is adapted from the CREST implementation. There are two reservoirs, one for overland/surface runoff and one for subsurface runoff. -

      -
    2. Kinematic Wave
      3. -

      Kinematic wave routing is a simplified approximation of the Barré de Saint-Venant equations developed in 1871. Kinematic wave routing assumes the gravity force and friction force are equal and cancel while neglecting the acceleration terms.

      -
    -
  4. Snow Melt Models
    5. -
      -
    1. Snow-17
      2. -

      Snow-17 is a temperature index based snow melt module widely used by the U.S. NWS. -

      -
    -
  5. Inundation Models
    6. -
      -
    1. Simple Inundation
      2. -

      -

      -
    -
  6. Compiling EF5
    7. -

    - EF5 makes use of the TIFF and GEOTIFF libraries for use as a raster format. You can obtain binaries or source code and compilation documentation from:
    - Compiling EF5 can be accomplished in 3 steps:
    - 1. This step is only needed if you are a developer and are adding files to the Makefile.am 

    autoreconf --force --install

    - 2. This sets up the system, if you have a path where you would like to install the files then use ./configure --prefix=/path/to/someplace 
    ./configure

    - 3. This step actually compiles ef5 and generates the binary 
    make CXXFLAGS="-O3 -fopenmp"

    - Upon successful compilation there will be a binary called "ef5" in the "bin" directory. -

    -
  7. Configuration File
    8. -

    The configuration file specifies all of the user changeable settings for EF5. Information in this file controls the input forcings, output options, run methods, etc.

    -

    In general the configuration file is case insensitive. The only exception is file paths when working on case sensitive file systems such as those typically found in Linux/Unix.

    -

    The configuration file supports three different styles of comments. The three styles are bash (#), C (/**/) and C++ (//).

    -
    #All variables/names/etc with the exception of
-#file paths are case insensitive
-//Multiple comment types are supported
-/*
-   Including multi-line C-style comments
- */
-

    Example of the different comment styles.

    -
      -
    1. Basic Information
      2. -

      The basic section of the configuration file specifies the digital elevation model (DEM), drainage direction map (DDM), and flow accumulation map (FAM) files.
      -

      -[Basic]
-DEM=/EF5Demo/FF/basic/DEM.asc
-DDM=/EF5Demo/FF/basic/DDM.asc
-FAM=/EF5Demo/FF/basic/FAM.asc
-PROJ=laea
-ESRIDDM=true
-SELFFAM=true
-

      Example Basic Block

      -DEM: Specifies the location and file name of the DEM grid in ESRI ascii or float32 geotiff format.
      -DDM: Specifies the location and file name of the DDM grid in ESRI ascii or float32 geotiff format.
      -FAM: Specifies the location and file name of the FAM grid in ESRI ascii or float32 geotiff format.
      -PROJ: Specifies the projection that the model will be expecting for input files. Possible values are:
      -
      -Geographic: Standard geographic projection
-LAEA: Lambert Azimuthal Equal Area projection with the standard parallel at 45.0 and the central meridian at -100.0. 
-
      -ESRIDDM: Specifies if the DDM is in ESRI format or TauDEM format. Possible values are:
      -
      -true: The directions in the DDM are specified as 
      3264128
      16 1
      842
       
-false: The directions in the DDM are specified as 
      432
      5 1
      678
      
-
      -SELFFAM: Specifies if the flow accumulation map includes the current grid cell in the flow count. Possible values are:
      -
      -true: The lowest flow accumulation value for any grid cell will be 1.
-false: The lowest flow accumulation value for any grid cell will be 0.
-
      -

      -
    2. Precipitation Information
      3. -

      The precipitation forcing section specifies the information necessary to adequately describe the precipitation product that the model will ingest.
      -

      -[PrecipForcing Q2Precip]
-TYPE=BIF
-UNIT=mm/h
-FREQ=5u
-LOC=/EF5Demo/FF/precip
-NAME=Q2_YYYYMMDDHHUU.bif
-

      Example PrecipForcing Block

      -PrecipForcing Block Name: This specifies in the name of the precipitation forcing block as it will be referred to later in the configuration file. In the above example the block name is "Q2Precip".
      -TYPE: Specifies the file type of the precipitation. Possible type values are:
      -
      -ASC: An ESRI ASCII grid.
-BIF: A binary file version of an ESRI ASCII grid.
-TIF: A float32 geotiff grid.
-TRMMRT: TRMM Multisatellite Precipitation Analysis realtime binary grid. Can be gzip compressed.
-TRMMV7: TRMM Multisatellite Precipitation Analysis 3B42V7 HDF5 grid.
-MRMS: Multi-Radar Multi-Sensor binary grid.
-
      -UNIT: Specifies the units of the precipitation in the file. Supported length units are meters (m), centimeters (cm) and millimeters (mm). Supported time units are year (y), month (m), day (d), hour (h), minute (u) and second (s). Modifiers in front of the time portion are also supported. For example if your precipitation forcing file has units of millimeters per three hours then your "UNIT" line would appear as "UNIT=mm/3h".
      -FREQ: Specifies the frequency at which precipitation files should be ingested by the model. Supported time units are year (y), month (m), day (d), hour (h), minute (u) and second (s).
      -LOC: Specifies the directory location of the precipitation forcing files.
      -NAME: Specifies the naming convention of the precipitation forcing files. These files can (and should) contain valid time dates. The name can be of any format. YYYY will be replaced with the year, MM replaced with the month, DD replaced with the day, HH replaced with the hour, UU replaced with the minute and SS replaced with the second.

      -
    3. Potential Evapotranspiration (PET) Information
      4. -

      The potential evapotranspiration forcing section specifies the information necessary to adequately describe the PET product that the model will ingest.

      -
      -[PETForcing PET]
-TYPE=BIF
-UNIT=mm/h
-FREQ=m
-LOC=/EF5Demo/FF/pet
-NAME=PET_MM.bif
-

      Example PETForcing Block

      -PETForcing Block Name: This specifies in the name of the PET forcing block as it will be referred to later in the configuration file. In the above example the block name is "PET".
      -TYPE: Specifies the file type of the PET. Possible type values are:
      -
      -ASC: An ESRI ASCII grid.
-BIF: A binary version of an ESRI ASCII grid.
-TIF: A float32 geotiff grid.
-
      -UNIT: Specifies the units of the PET in the file. Supported length units are meters (m), centimeters (cm) and millimeters (mm). Support -ed time units are year (y), month (m), day (d), hour (h), minute (u) and second (s). Modifiers in front of the time portion are also supported. For example if your PET forcing file has units of millimeters per three hours then your "UNIT" line would appear as "UNIT=mm/3h".
      PET data may also be given as temperate data in degrees Celsius with unit "C". The temperature data is converted into PET.
      -FREQ: Specifies the frequency at which PET files should be ingested by the model. Supported time units are year (y), month (m), day (d -), hour (h), minute (u) and second (s).
      -LOC: Specifies the directory location of the PET forcing files.
      -NAME: Specifies the naming convention of the PET forcing files. These files can (and should) contain valid time dates. The name can be - of any format. YYYY will be replaced with the year, MM replaced with the month, DD replaced with the day, HH replaced with the hour, UU replaced with the minute and SS replaced with the second.

      -
    4. Gauge Locations
      5. -

      These blocks specify the location of the gauges to the model. This is useful if you want time series output at a point and also to specify parameters. Basins are treated as collections of gauges with the outlet gauge being the independent gauge and all other gauges inside a basin being dependent gauges.

      -
      [Gauge OKC]
-LON=-97.01
-LAT=35.68
-OBS=/EF5Demo/obs/okc.csv
-BASINAREA=341.88
-OUTPUTTS=TRUE
-
-[Gauge AR]
-LON=-93.62
-LAT=34.37
-

      Example Gauge Location Blocks

      -Gauge Location Block Name: This specifies in the name of the gauge location block as it will be referred to later in the configuration file. In the above example the block names are "OKC" and "AR".
      -LON: This is the longitude of the gauge in an unprojected coordinate system. EF5 will reproject the gauge point into the appropriate coordinate system.
      -LAT: This is the latitude of the gauge in an unprojected coordinate system. EF5 will reproject the gauge point into the appropriate coordinate system.
      -CELLX: (optional) This is the x-coordinate of the gauge in the basic file, used in place of LAT-LON.
      -CELLY: (optional) This is the y-coordinate of the gauge in the basic file, used in place of LAT-LON.
      -OBS: (optional) This is a file containing a time series specifying the observed values at the gauge. This is only used during model calibration.
      -BASINAREA: (optional) Observed contributing basin area for this gauge in km2. EF5 will search nearby grid cells to match flow accumulations as closely as possible.
      -OUPUTTS: (optional) Output a time series file for this gauge. Possible values are TRUE, YES, FALSE and NO. Default value is YES
      -WANTDA: (optional) If we are doing data assimilation, do we want it done for this gauge? Possible values are TRUE, YES, FALSE and NO. Default value is YES
      -WANTCO: (optional) Do we want this gauge included in the combined output file? Possible values are TRUE, YES, FALSE and NO. Default value is NO

      - -
    5. Basins
      6. -

      These blocks do not represent actual physical basins but rather a collection of specified gauge locations that may or may not form a single independent basin.

      -
      [Basin FF]
-GAUGE=OKC
-GAUGE=AR
-

      Example Basin Block

      -Basin Block Name: This specifies in the name of the basin block as it will be referred to later in the configuration file. In the above example the block name is "FF".
      -GAUGE: This is the name of the gauge location block to include as part of this basin.

      - -
    6. Parameter Sets
      7. -

      These blocks control the distributed model parameter settings. Parameters are specified per gauge and propogated upstream. Therefore each independent gauge must have a parameter set specified. Gridded parameters can be specified using the parameter name combined with "_grid" with a value specifying the grid file. When gridded parameters and uniform parameters are specified the uniform parameters act as scalar multipliers on the gridded parameter values.

      -
        -
      1. CREST Parameter Set
        2. -
        -[CrestParamSet ABRFC]
-wm_grid=/hydros/zflamig/EF5_Jess/params/crest_params/wm_usa.tif
-im_grid=/hydros/zflamig/EF5_Jess/params/crest_params/im_usa.tif
-fc_grid=/hydros/zflamig/EF5_Jess/params/crest_params/ksat_usa.tif
-b_grid=/hydros/zflamig/EF5_Jess/params/crest_params/b_usa.tif
-gauge=03455500
-wm=1.00
-b=1.0
-im=0.01
-ke=1.0
-fc=1.00
-iwu=50.0
-

        Example CrestParamSet Block

        - GAUGE: This is the name of the gauge location block for which the parameters that follow it will be applied.
        - WM: The maximum soil water capacity (depth integrated pore space) of the soil layer. (mm)
        - B: The exponent of the variable infiltration curve.
        - IM: The impervious area ratio. (% between 0 and 100)
        - KE: The multiplier to convert between input PET and local actual ET.
        - FC: The soil saturated hydraulic conductivity (Ksat). (mm hr-1)
        - IWU: The initial value of soil water. This is a percentage of the WM. (% between 0 and 100)

        - -
      2. SAC-SMA Parameter Set
        3. - 
        -[SacParamSet ABRFC]
-UZTWM_grid=/data/FLASH/v2/sac_params/uztwm_usa.tif
-UZFWM_grid=/data/FLASH/v2/sac_params/uzfwm_usa.tif
-UZK_grid=/data/FLASH/v2/sac_params/uzk_usa.tif
-ZPERC_grid=/data/FLASH/v2/sac_params/zperc_usa.tif
-REXP_grid=/data/FLASH/v2/sac_params/rexp_usa.tif
-LZTWM_grid=/data/FLASH/v2/sac_params/lztwm_usa.tif
-LZFSM_grid=/data/FLASH/v2/sac_params/lzfsm_usa.tif
-LZFPM_grid=/data/FLASH/v2/sac_params/lzfpm_usa.tif
-LZSK_grid=/data/FLASH/v2/sac_params/lzsk_usa.tif
-LZPK_grid=/data/FLASH/v2/sac_params/lzpk_usa.tif
-PFREE_grid=/data/FLASH/v2/sac_params/pfree_usa.tif
-gauge=01055000
-UZTWM=1.0
-UZFWM=1.0
-UZK=1.0
-PCTIM=0.101
-ADIMP=0.10
-RIVA=1.001
-ZPERC=1.0
-REXP=1.0
-LZTWM=1.0
-LZFSM=1.0
-LZFPM=1.0
-LZSK=1.0
-LZPK=1.0
-PFREE=1.0
-SIDE=0.0
-RSERV=0.3
-ADIMC=1.0
-UZTWC=0.55
-UZFWC=0.14
-LZTWC=0.56
-LZFSC=0.11
-LZFPC=0.46
-

        Example SacParamSet Block

        - GAUGE: This is the name of the gauge location block for which the parameters that follow it will be applied.
        - UZTWM: Upper zone tension water maximum (mm)
        - UZFWM: Upper zone free water maximum (mm)
        - UZK: Fractional daily upper zone free water withdrawal rate
        - PCTIM: Minimum impervious area (%)
        - ADIMP: Additional impervious area (%)
        - RIVA: Riparian vegetation area (%)
        - ZPERC: Maximum percolation rate
        - REXP: Exponent for the percolation equation
        - LZTWM: Lower zone tension water capacity (mm)
        - LZFSM: Lower zone supplemental free water capacity (mm)
        - LZFPM: Lower zone primary free water capacity (mm)
        - LZSK: Fractional daily supplemental withdrawal rate
        - LZPK: Fractional daily primary withdrawal rate
        - PFREE: Percent of percolated water which always goes directly to lower zone free water storages (%)
        - SIDE: Ratio of non-channel baseflow (deep recharge) to channel (visible) baseflow
        - RSERV: Percent of lower zone free water which cannot be transferred to lower zone tension water (%)
        - ADIMC: Tension water contents of the ADIMP area (mm)
        - UZTWC: Upper zone tension water contents (mm)
        - UZFWC: Upper zone free water contents (mm)
        - LZTWC: Lower zone tension water contents (mm)
        - LZFSC: Lower zone free supplemental contents (mm)
        - LZFPC: Lower zone free primary contents (mm)

        -
      3. HP Parameter Set
        4. - To be completed in a future revision.

        -
      4. Linear Reservoir Parameter Set
        5. - 
        -[lrparamset rundu]
-gauge=rundu
-coem=1611.115479
-river=307.980042
-under=2531.556641
-leako=0.918236
-leaki=0.017568
-th=8.140809
-iso=0.000040
-isu=0.000073
-
        - To be completed in a future revision

        -
      5. Kinematic Wave Parameter Set
        6. - 
        -[KWParamSet rundu]
-GAUGE=rundu
-UNDER=1.673110
-LEAKI=0.043105
-TH=6.658569
-ISU=0.000000
-ALPHA=2.991570
-BETA=0.932080
-ALPHA0=4.603945
-

        Example KWParamSet Block

        - GAUGE: This is the name of the gauge location block for which the parameters that follow it will be applied.
        - TH: The number of grid cells needed to flow into a cell for it to be declared a channel grid cell.
        - UNDER: The interflow flow speed multiplier
        - LEAKI: The amount of water leaked out of the interflow reservoir at each time step. Ranges between 0 and 1.
        - ISU: The initial value of the interflow reservoir.
        - ALPHA: The multiplier in Q = alpha*Abeta used for channel routing
        - BETA: The exponent in Q = alpha*Abeta used for channel routing
        - ALPHA0: The Alpha value used for overland routing.

        -
      6. Snow-17 Parameter Set
        7. - 
        -[snow17paramset tarbela]
-GAUGE=tarbela
-UADJ=0.184653
-MBASE=0.047224
-MFMAX=1.068658
-MFMIN=0.516059
-TIPM=0.911706
-NMF=0.077336
-PLWHC=0.093812
-SCF=2.219492
-

        Example Snow17ParamSet Block

        - GAUGE: This is the name of the gauge location block for which the parameters that follow it will be applied.
        - UADJ: The number of grid cells needed to flow into a cell for it to be declared a channel grid cell.
        - MBASE: The interflow flow speed multiplier
        - MFMAX: The amount of water leaked out of the interflow reservoir at each time step. Ranges between 0 and 1.
        - MFMIN: The initial value of the interflow reservoir.
        - TIPM: The multiplier in Q = alpha*Abeta used for channel routing
        - NMF: The exponent in Q = alpha*Abeta used for channel routing
        - PLWHC: The Alpha value used for overland routing.
        - SCF: The Alpha value used for overland routing.

        -
      7. Simple Inundation Parameter Set
        8. - 
        -[simpleinundationparamset rundu]
-gauge=rundu
-alpha=2.991570
-beta=0.932080
-

        Example SimpleInundationParamSet Block

        - GAUGE: This is the name of the gauge location block for which the parameters that follow it will be applied.
        - ALPHA: The multiplier in A = (Q/alpha)1/beta used for computing cross-sectional area
        - BETA: The exponent in A = (Q/alpha)1/beta used for computing cross-sectional area

        -
      -
    7. Tasks
      8. -

      Tasks define the necessary information about which model to run, over what time period to run, with what time step to run, etc.

      -
      [Task RunFF]
-STYLE=SIMU
-MODEL=CREST
-BASIN=FF
-PRECIP=Q2_PRECIP
-PET=PET
-OUTPUT=/EF5Demo/FF/output/
-PARAM_SET=FF
-TIMESTEP=5u
-TIME_BEGIN=201006010000
-TIME_END=201006010030
-

      Example Task Block

      - - STYLE: The style of task that this is. Possible values are:
      - 
      -SIMU: A simulation run.
-SIMU_RP: A simulation run that will produce standard deviation, average and skew coefficient grids for estimating return period.
-CALI_DREAM: A calibration run using DREAM.
-CLIP_BASIN: Clips the basic files to the specified BASIN.
-CLIP_GAUGE: Clips the basic files to the first specified gauge.
      - MODEL: The water balance model that this task should use. Possible values are:
      - 
      -CREST: The Coupled Routing and Excess STorage distributed hydrological model.
-SAC: The SAC-SMA model.
-HyMOD: HyMod.
      - ROUTING: The routing method that this task should use. Possible values are:
      - 
      -KW: Kinematic Wave Routing.
-LR: Linear Reservoir Routing.
      - SNOW: (Optional) The snow melt model that this task should use. Possible values are:
      - 
      SNOW17: The Snow-17 snow melt model.
      - INUNDATION: (Optional) The inundation model that this task should use. Possible values are:
      - 
      SIMPLEINUNDATIOn: The Simple Inundation model.
      - BASIN: The basin block name which defines the area over which the model should be run.
      - DEFAULTPARAMSGAUGE: (Optional) The gauge for which parameters are specified and used for gauges which did not have parameters specified. Useful when modeling large areas with gridded parameters.
      - PRECIP: The precipitation block name which defines the precipitation.
      - PRECIPFORECAST: (Optional) The precipitation block name which defines the forecast precipitation. Used if the precipiation specified in PRECIP is unavailable for the current time.
      - TEMP: (Required if using SNOW) The temperature block name which defines the temperature data to be used.
      - TEMPFORECAST: (Optional if using SNOW) The temperature block name which defines the forecast temperatures. Used if the temperature specified in TEMP is unavailable for the current time.
      - PET: The PET block name which defines the PET to use.
      - PARAM_SET: The parameter set block name which defines which set of water balance parameters to use.
      - ROUTING_PARAM_SET: The parameter set block name which defines which set of routing parameters to use.
      - SNOW_PARAM_SET: (Required if using SNOW) The parameter set block name which defines which set of snow parameters to use.
      - INUNDATION_PARAM_SET: (Required if using INUNDATION) The parameter set block name which defines which set of inundation parameters to use.
      - CALI_PARAM: (Required if using CALI_DREAM) The parameter set block name which defines which set of water balance parameters settings for calibration.
      - ROUTING_CALI_PARAM: (Required if using CALI_DREAM) The parameter set block name which defines which set of routing parameters to use for calibration.
      - SNOW_CALI_PARAM: (Required if using SNOW, CALI_DREAM) The parameter set block name which defines which set of snow parameters to use for calibration.
      - INUNDATION_CALI_PARAM: (Required if using INUNDATION, CALI_DREAM) The parameter set block name which defines which set of inundation parameters to use for calibration.
      - PRELOAD_FILE: (Optional) The file path and name where for the preload file. The preload file contains the forcings (Precip, PET, Temp) defined for the current time period and basin extent. Generated by EF5 if it does not exist. Useful for faster runs when forcings are not changing such as with manual calibration.
      - STATES: (Optional) The location where output files should be written.
      - TIMESTEP: The time step to use when running the model. Supported time units are year (y), month (m), day (d), hour (h), minute (u) and second (s).
      - TIME_BEGIN: The initialization time for the model run. YYYYMMDDHHUUSS format.
      - TIME_END: The ending time for the model run. YYYYMMDDHHUUSS format.
      - TIME_WARMEND: The end of the warm-up period. YYYYMMDDHHUUSS format.
      - TIME_STATE: (Optional) The time at which to output model states. YYYYMMDDHHUUSS format.
      - OUTPUT: The location where output files should be written.
      - OUTPUT_GRIDS: (Optional) Which grids should be output, combine together with | :
      - 
      -NONE: No gridded output.
-STREAMFLOW: Streamflow (cms)
-SOILMOISTURE: Soil moisture (%)
-RETURNPERIOD: Streamflow return period (years)
-PRECIP: Precipitaiton (mm)
-PET: Potential evapotranspiration (mm)
-SNOWWATER: Snow water equivalent from the snow melt model (mm)
-TEMPERATURE: Temperature (C)
-INUNDATION: Water depth (m)
-MAXSTREAMFLOW: Maximum streamflow during run (cms)
-MAXSOILMOISTURE: Maximum soil moisture during run (%)
-MAXRETURNPERIOD: Maximum return period during run (years)
-MAXSNOWWATER: Maximum snow water equivalent (mm)
      - DA_FILE: (Optional) The input observations to be added for use with streamflow data assimilation,
      - CO_FILE: (Optional) The combined output file if one is desired.
      - RP_STDGRID: (Optional) The geotiff file representing the standard deviation for the Log-Pearson Type III return period distribution.
      - RP_AVGGRID: (Optional) The geotiff file representing the average for the Log-Pearson Type III return period distribution.
      - RP_CSGRID: (Optional) The geotiff file representing the skew coefficient for the Log-Pearson Type III return period distribution.
      -
    8. Execute Block
      9. -

      The execute block defines which tasks to run when the ef5 program is executed

      -
      [Execute]
-TASK=RunFF
-

      Example Execute Block

      - TASK: The task block name which should be executed.

      -
    -
  8. Running EF5
    9. - 
    ef5 [controlfile]
    -

    Running ef5 is straight forwarded, you can called the executable with no arguments and it will assume a default configuration file name of "control.txt" in the current working directory or you can pass in a configuration file name as an argument.

    -
  9. Calibrating the Models
    10. -

    This section will be filled in once calibration methods are implemented.

    -
  10. Appendix
    11. -
      -
    1. Complete Sample Configuration File
      2. - 
      -/*
- * This is an example configuration file for EF5
- */
-
-[Basic]
-DEM=/EF5Demo/FF/basic/DEM.asc
-DDM=/EF5Demo/FF/basic/DDM.asc
-FAM=/EF5Demo/FF/basic/FAM.asc
-PROJ=laea
-ESRIDDM=true
-
-[PrecipForcing Q2Precip]
-TYPE=BIF
-UNIT=mm/h
-FREQ=5u
-LOC=/EF5Demo/FF/precip
-NAME=Q2_YYYYMMDDHHUU.bif
-
-[PETForcing PET]
-TYPE=BIF
-UNIT=mm/h
-FREQ=m
-LOC=/EF5Demo/FF/pet
-NAME=PET_MM.bif
-
-[Gauge OKC]
-LON=-97.01
-LAT=35.68
-OBS=/EF5Demo/obs/okc.csv
-
-[Gauge AR]
-LON=-93.62
-LAT=34.37
-
-[Basin FF]
-GAUGE=OKC
-GAUGE=AR
-
-[CrestParamSet FF]
-GAUGE=AR
-COEM=24.230076 EXPM=0.502391 RIVER=1.73056
-UNDER=0.291339 LEAKO=0.56668 LEAKI=0.251648
-TH=63.20205 GM=1.364364 PWM=71.96465
-PB=0.964355 PIM=6.508687 PKE=0.19952
-PFC=2.578529 IWU=53.52593 ISO=5.899539
-ISU=17.31128
-GAUGE=OKC
-COEM=24.230076 EXPM=0.502391 RIVER=1.73056
-UNDER=0.291339 LEAKO=0.56668 LEAKI=0.251648
-TH=63.20205 GM=1.364364 PWM=71.96465
-PB=0.964355 PIM=6.508687 PKE=0.19952
-PFC=2.578529 IWU=53.52593 ISO=5.899539
-ISU=17.31128
-
-[Task RunFF]
-STYLE=SIMU
-MODEL=CREST
-BASIN=FF
-PRECIP=Q2_PRECIP
-PET=PET
-OUTPUT=/EF5Demo/FF/output/
-PARAM_SET=FF
-TIMESTEP=5u
-TIME_BEGIN=201006010000
-TIME_END=201006010030
-
-[Execute]
-TASK=RunFF

      Example of a full EF5 configuration file.

      -
    -
-
- - + + + EF5 User Manual + + + + + +
+
    +
  1. About EF5
    2. +

    EF5 is designed to facilitate creation of ensemble forecasts for flash flood prediction. As such it will incorporate multi-model support while maintaining a single set of input data. Currently the only supported model is the Coupled Routing and Excess Storage (CREST) hydrologic model. Additionally, EF5 was designed with utilization of parallel computing in mind. Presently portions of CREST are optimized to take advantage of multi-core computing through the use of OpenMP.

    +
  2. Hydrologic Water Balance Models
    3. +
      +
    1. CREST
      2. +

      The Coupled Routing and Excess STorage (CREST) distributed hydrological model is a hybrid modeling strategy that was recently developed by the University of Oklahoma (http://hydro.ou.edu) and NASA SERVIR Project Team (www.servir.net). CREST simulates the spatiotemporal variation of water and energy fluxes and storages on a regular grid with the grid cell resolution being user-defined, thereby enabling global- and regional-scale applications. The scalability of CREST simulations is accomplished through sub-grid scale representation of soil moisture storage capacity (using a variable infiltration curve) and runoff generation processes (using linear reservoirs). The CREST model was initially developed to provide online global flood predictions with relatively coarse resolution, but it is also applicable at small scales, such as single basins. The CREST Model can be forced by gridded potential evapotranspiration and precipitation datasets such as, satellite-based precipitation estimates, gridded rain gauge observations, remote sensing platforms such as weather radar, and quantitative precipitation forecasts from numerical weather prediction models. The representation of the primary water fluxes such as infiltration and routing are closely related to the spatially variable land surface characteristics (i.e., vegetation, soil type, and topography). The runoff generation component and routing scheme are coupled, thus providing realistic interactions between atmospheric, land surface, and subsurface water. +

      More detailed information about CREST can be found in the following publication:
      Wang, J., Y. Hong, L. Li, J. J. Gourley, S. I. Khan, K. K. Yilmaz, R. F. Adler, F. S. Policelli, S. Habib, D. Irwin, A. S. Limaye, T. Korme, and L. Okello, 2011: The coupled routing and excess storage (CREST) distributed hydrological model. Hydrol. Sci. Journal, 56, 84-98, doi: 10.1080/02626667.2010.543087. +

      +
    2. CRESTPHYS
      3. +

      CRESTPHYS builds upon the CREST model. Differently, it separates interflow with baseflow by using a fill-spill bucket to represent conceptual groundwater reservior. This approach is applied also in the National Water Model. We hope this enhanced module can bring benefits in representing baseflow and flow recession limb. +

      More detailed information about CRESTPHYS will follow... +

      Contact Allen Zhi Li (li1995@ou.edu) if you have questions regarding this module +

      +
    3. SAC-SMA
      4. +

      The Sacramento Soil Moisture Accounting (SAC-SMA) Model was developed by the U.S. National Weather Service with the goal of parameterizing soil moisture characteristics in a fashion that would 1) logically distribute applied moisture in various depths and energy states in the soil 2) have rational percolation characteristics 3) allow an effective simulation of streamflow. +

      More detailed information about SAC-SMA can be found online at U.S. NWS SAC-SMA algorithm description

      +
    4. HP
      5. +

      The Hydrophobic (HP) water balance model features an entirely impervious surface where all rainfall is transformed into surface runoff.

      +
    +
  3. Routing Models
    4. +
      +
    1. Linear Reservoir
      2. +

      + The linear reservoir routing is adapted from the CREST implementation. There are two reservoirs, one for overland/surface runoff and one for subsurface runoff. +

      +
    2. Kinematic Wave
      3. +

      Kinematic wave routing is a simplified approximation of the Barré de Saint-Venant equations developed in 1871. Kinematic wave routing assumes the gravity force and friction force are equal and cancel while neglecting the acceleration terms.

      +
    +
  4. Snow Melt Models
    5. +
      +
    1. Snow-17
      2. +

      Snow-17 is a temperature index based snow melt module widely used by the U.S. NWS. +

      +
    +
  5. Inundation Models
    6. +
      +
    1. Simple Inundation
      2. +

      +

      +
    +
  6. Compiling EF5
    7. +

    + EF5 makes use of the TIFF and GEOTIFF libraries for use as a raster format. You can obtain binaries or source code and compilation documentation from:
    + Compiling EF5 can be accomplished in 3 steps:
    + 1. This step is only needed if you are a developer and are adding files to the Makefile.am 

    autoreconf --force --install

    + 2. This sets up the system, if you have a path where you would like to install the files then use ./configure --prefix=/path/to/someplace 
    ./configure

    + 3. This step actually compiles ef5 and generates the binary 
    make CXXFLAGS="-O3 -fopenmp"

    + Upon successful compilation there will be a binary called "ef5" in the "bin" directory. +

    +
  7. Configuration File
    8. +

    The configuration file specifies all of the user changeable settings for EF5. Information in this file controls the input forcings, output options, run methods, etc.

    +

    In general the configuration file is case insensitive. The only exception is file paths when working on case sensitive file systems such as those typically found in Linux/Unix.

    +

    The configuration file supports three different styles of comments. The three styles are bash (#), C (/**/) and C++ (//).

    +
    #All variables/names/etc with the exception of
+#file paths are case insensitive
+//Multiple comment types are supported
+/*
+   Including multi-line C-style comments
+ */
+

    Example of the different comment styles.

    +
      +
    1. Basic Information
      2. +

      The basic section of the configuration file specifies the digital elevation model (DEM), drainage direction map (DDM), and flow accumulation map (FAM) files.
      +

      +[Basic]
+DEM=/EF5Demo/FF/basic/DEM.asc
+DDM=/EF5Demo/FF/basic/DDM.asc
+FAM=/EF5Demo/FF/basic/FAM.asc
+PROJ=laea
+ESRIDDM=true
+SELFFAM=true
+

      Example Basic Block

      +DEM: Specifies the location and file name of the DEM grid in ESRI ascii or float32 geotiff format.
      +DDM: Specifies the location and file name of the DDM grid in ESRI ascii or float32 geotiff format.
      +FAM: Specifies the location and file name of the FAM grid in ESRI ascii or float32 geotiff format.
      +PROJ: Specifies the projection that the model will be expecting for input files. Possible values are:
      +
      +Geographic: Standard geographic projection
+LAEA: Lambert Azimuthal Equal Area projection with the standard parallel at 45.0 and the central meridian at -100.0. 
+
      +ESRIDDM: Specifies if the DDM is in ESRI format or TauDEM format. Possible values are:
      +
      +true: The directions in the DDM are specified as 
      3264128
      16 1
      842
       
+false: The directions in the DDM are specified as 
      432
      5 1
      678
      
+
      +SELFFAM: Specifies if the flow accumulation map includes the current grid cell in the flow count. Possible values are:
      +
      +true: The lowest flow accumulation value for any grid cell will be 1.
+false: The lowest flow accumulation value for any grid cell will be 0.
+
      +

      +
    2. Precipitation Information
      3. +

      The precipitation forcing section specifies the information necessary to adequately describe the precipitation product that the model will ingest.
      +

      +[PrecipForcing Q2Precip]
+TYPE=BIF
+UNIT=mm/h
+FREQ=5u
+LOC=/EF5Demo/FF/precip
+NAME=Q2_YYYYMMDDHHUU.bif
+

      Example PrecipForcing Block

      +PrecipForcing Block Name: This specifies in the name of the precipitation forcing block as it will be referred to later in the configuration file. In the above example the block name is "Q2Precip".
      +TYPE: Specifies the file type of the precipitation. Possible type values are:
      +
      +ASC: An ESRI ASCII grid.
+BIF: A binary file version of an ESRI ASCII grid.
+TIF: A float32 geotiff grid.
+TRMMRT: TRMM Multisatellite Precipitation Analysis realtime binary grid. Can be gzip compressed.
+TRMMV7: TRMM Multisatellite Precipitation Analysis 3B42V7 HDF5 grid.
+MRMS: Multi-Radar Multi-Sensor binary grid.
+
      +UNIT: Specifies the units of the precipitation in the file. Supported length units are meters (m), centimeters (cm) and millimeters (mm). Supported time units are year (y), month (m), day (d), hour (h), minute (u) and second (s). Modifiers in front of the time portion are also supported. For example if your precipitation forcing file has units of millimeters per three hours then your "UNIT" line would appear as "UNIT=mm/3h".
      +FREQ: Specifies the frequency at which precipitation files should be ingested by the model. Supported time units are year (y), month (m), day (d), hour (h), minute (u) and second (s).
      +LOC: Specifies the directory location of the precipitation forcing files.
      +NAME: Specifies the naming convention of the precipitation forcing files. These files can (and should) contain valid time dates. The name can be of any format. YYYY will be replaced with the year, MM replaced with the month, DD replaced with the day, HH replaced with the hour, UU replaced with the minute and SS replaced with the second.

      +
    3. Potential Evapotranspiration (PET) Information
      4. +

      The potential evapotranspiration forcing section specifies the information necessary to adequately describe the PET product that the model will ingest.

      +
      +[PETForcing PET]
+TYPE=BIF
+UNIT=mm/h
+FREQ=m
+LOC=/EF5Demo/FF/pet
+NAME=PET_MM.bif
+

      Example PETForcing Block

      +PETForcing Block Name: This specifies in the name of the PET forcing block as it will be referred to later in the configuration file. In the above example the block name is "PET".
      +TYPE: Specifies the file type of the PET. Possible type values are:
      +
      +ASC: An ESRI ASCII grid.
+BIF: A binary version of an ESRI ASCII grid.
+TIF: A float32 geotiff grid.
+
      +UNIT: Specifies the units of the PET in the file. Supported length units are meters (m), centimeters (cm) and millimeters (mm). Support +ed time units are year (y), month (m), day (d), hour (h), minute (u) and second (s). Modifiers in front of the time portion are also supported. For example if your PET forcing file has units of millimeters per three hours then your "UNIT" line would appear as "UNIT=mm/3h".
      PET data may also be given as temperate data in degrees Celsius with unit "C". The temperature data is converted into PET.
      +FREQ: Specifies the frequency at which PET files should be ingested by the model. Supported time units are year (y), month (m), day (d +), hour (h), minute (u) and second (s).
      +LOC: Specifies the directory location of the PET forcing files.
      +NAME: Specifies the naming convention of the PET forcing files. These files can (and should) contain valid time dates. The name can be + of any format. YYYY will be replaced with the year, MM replaced with the month, DD replaced with the day, HH replaced with the hour, UU replaced with the minute and SS replaced with the second.

      +
    4. Gauge Locations
      5. +

      These blocks specify the location of the gauges to the model. This is useful if you want time series output at a point and also to specify parameters. Basins are treated as collections of gauges with the outlet gauge being the independent gauge and all other gauges inside a basin being dependent gauges.

      +
      [Gauge OKC]
+LON=-97.01
+LAT=35.68
+OBS=/EF5Demo/obs/okc.csv
+BASINAREA=341.88
+OUTPUTTS=TRUE
+
+[Gauge AR]
+LON=-93.62
+LAT=34.37
+

      Example Gauge Location Blocks

      +Gauge Location Block Name: This specifies in the name of the gauge location block as it will be referred to later in the configuration file. In the above example the block names are "OKC" and "AR".
      +LON: This is the longitude of the gauge in an unprojected coordinate system. EF5 will reproject the gauge point into the appropriate coordinate system.
      +LAT: This is the latitude of the gauge in an unprojected coordinate system. EF5 will reproject the gauge point into the appropriate coordinate system.
      +CELLX: (optional) This is the x-coordinate of the gauge in the basic file, used in place of LAT-LON.
      +CELLY: (optional) This is the y-coordinate of the gauge in the basic file, used in place of LAT-LON.
      +OBS: (optional) This is a file containing a time series specifying the observed values at the gauge. This is only used during model calibration.
      +BASINAREA: (optional) Observed contributing basin area for this gauge in km2. EF5 will search nearby grid cells to match flow accumulations as closely as possible.
      +OUPUTTS: (optional) Output a time series file for this gauge. Possible values are TRUE, YES, FALSE and NO. Default value is YES
      +WANTDA: (optional) If we are doing data assimilation, do we want it done for this gauge? Possible values are TRUE, YES, FALSE and NO. Default value is YES
      +WANTCO: (optional) Do we want this gauge included in the combined output file? Possible values are TRUE, YES, FALSE and NO. Default value is NO

      + +
    5. Basins
      6. +

      These blocks do not represent actual physical basins but rather a collection of specified gauge locations that may or may not form a single independent basin.

      +
      [Basin FF]
+GAUGE=OKC
+GAUGE=AR
+

      Example Basin Block

      +Basin Block Name: This specifies in the name of the basin block as it will be referred to later in the configuration file. In the above example the block name is "FF".
      +GAUGE: This is the name of the gauge location block to include as part of this basin.

      + +
    6. Parameter Sets
      7. +

      These blocks control the distributed model parameter settings. Parameters are specified per gauge and propogated upstream. Therefore each independent gauge must have a parameter set specified. Gridded parameters can be specified using the parameter name combined with "_grid" with a value specifying the grid file. When gridded parameters and uniform parameters are specified the uniform parameters act as scalar multipliers on the gridded parameter values.

      +
        +
      1. CREST Parameter Set
        2. +
        +[CrestParamSet ABRFC]
+wm_grid=/hydros/zflamig/EF5_Jess/params/crest_params/wm_usa.tif
+im_grid=/hydros/zflamig/EF5_Jess/params/crest_params/im_usa.tif
+fc_grid=/hydros/zflamig/EF5_Jess/params/crest_params/ksat_usa.tif
+b_grid=/hydros/zflamig/EF5_Jess/params/crest_params/b_usa.tif
+gauge=03455500
+wm=1.00
+b=1.0
+im=0.01
+ke=1.0
+fc=1.00
+iwu=50.0
+

        Example CrestParamSet Block

        + GAUGE: This is the name of the gauge location block for which the parameters that follow it will be applied.
        + WM: The maximum soil water capacity (depth integrated pore space) of the soil layer. (mm)
        + B: The exponent of the variable infiltration curve.
        + IM: The impervious area ratio. (% between 0 and 100)
        + KE: The multiplier to convert between input PET and local actual ET.
        + FC: The soil saturated hydraulic conductivity (Ksat). (mm hr-1)
        + IWU: The initial value of soil water. This is a percentage of the WM. (% between 0 and 100)

        +
      2. CRESTPHYS Parameter Set
        3. +
        +  [CrestphysParamSet ABRFC]
+  gauge=03455500
+  wm_grid=/hydros/wm_usa.tif
+  im_grid=/hydros/im_usa.tif
+  fc_grid=/hydros/ksat_usa.tif
+  b_grid=/hydros/b_usa.tif
+  hmaxaq_grid=/hydros//params/hmaxaq.tif
+  gwe_grid=/hydros/params/gwe.tif
+  gwc_grid=/hydros/gwc.tif
+  wm=1.00
+  b=1.0
+  im=0.01
+  ke=1.0
+  ksoil=0.1
+  fc=1.00
+  iwu=50.0
+  igw=25.0
+  hmaxaq=1.0
+  gwc=1.0
+  gwe=1.0
+

        Example CrestPhysParamSet Block

        +GAUGE: This is the name of the gauge location block for which the parameters that follow it will be applied.
        +WM: The maximum soil water capacity (depth integrated pore space) of the soil layer. (mm)
        +B: The exponent of the variable infiltration curve.
        +IM: The impervious area ratio. (% between 0 and 100)
        +KE: The multiplier to convert between input PET and local actual ET.
        +KSOIL: The coefficient to drain soil moisture into groundwater (0-1)
        +FC: The soil saturated hydraulic conductivity (Ksat). (mm hr-1)
        +IWU: The initial value of soil water. This is a percentage of the WM. (% between 0 and 100)
        +IGW: The initial value of groundwater. This is a percentage of the HMAXAQ. (% between 0 and 100)
        +HMAXAQ: The maximum groundwater depth (m)
        +GWC: The coefficient (multiplier) to partition groundwater and baseflow
        +GWE: The exponent to partition groundwater and baseflow
        +
        + + +
      3. SAC-SMA Parameter Set
        4. + 
        +
+[SacParamSet ABRFC]
+UZTWM_grid=/data/FLASH/v2/sac_params/uztwm_usa.tif
+UZFWM_grid=/data/FLASH/v2/sac_params/uzfwm_usa.tif
+UZK_grid=/data/FLASH/v2/sac_params/uzk_usa.tif
+ZPERC_grid=/data/FLASH/v2/sac_params/zperc_usa.tif
+REXP_grid=/data/FLASH/v2/sac_params/rexp_usa.tif
+LZTWM_grid=/data/FLASH/v2/sac_params/lztwm_usa.tif
+LZFSM_grid=/data/FLASH/v2/sac_params/lzfsm_usa.tif
+LZFPM_grid=/data/FLASH/v2/sac_params/lzfpm_usa.tif
+LZSK_grid=/data/FLASH/v2/sac_params/lzsk_usa.tif
+LZPK_grid=/data/FLASH/v2/sac_params/lzpk_usa.tif
+PFREE_grid=/data/FLASH/v2/sac_params/pfree_usa.tif
+gauge=01055000
+UZTWM=1.0
+UZFWM=1.0
+UZK=1.0
+PCTIM=0.101
+ADIMP=0.10
+RIVA=1.001
+ZPERC=1.0
+REXP=1.0
+LZTWM=1.0
+LZFSM=1.0
+LZFPM=1.0
+LZSK=1.0
+LZPK=1.0
+PFREE=1.0
+SIDE=0.0
+RSERV=0.3
+ADIMC=1.0
+UZTWC=0.55
+UZFWC=0.14
+LZTWC=0.56
+LZFSC=0.11
+LZFPC=0.46
+

        Example SacParamSet Block

        + GAUGE: This is the name of the gauge location block for which the parameters that follow it will be applied.
        + UZTWM: Upper zone tension water maximum (mm)
        + UZFWM: Upper zone free water maximum (mm)
        + UZK: Fractional daily upper zone free water withdrawal rate
        + PCTIM: Minimum impervious area (%)
        + ADIMP: Additional impervious area (%)
        + RIVA: Riparian vegetation area (%)
        + ZPERC: Maximum percolation rate
        + REXP: Exponent for the percolation equation
        + LZTWM: Lower zone tension water capacity (mm)
        + LZFSM: Lower zone supplemental free water capacity (mm)
        + LZFPM: Lower zone primary free water capacity (mm)
        + LZSK: Fractional daily supplemental withdrawal rate
        + LZPK: Fractional daily primary withdrawal rate
        + PFREE: Percent of percolated water which always goes directly to lower zone free water storages (%)
        + SIDE: Ratio of non-channel baseflow (deep recharge) to channel (visible) baseflow
        + RSERV: Percent of lower zone free water which cannot be transferred to lower zone tension water (%)
        + ADIMC: Tension water contents of the ADIMP area (mm)
        + UZTWC: Upper zone tension water contents (mm)
        + UZFWC: Upper zone free water contents (mm)
        + LZTWC: Lower zone tension water contents (mm)
        + LZFSC: Lower zone free supplemental contents (mm)
        + LZFPC: Lower zone free primary contents (mm)

        +
      4. HP Parameter Set
        5. + To be completed in a future revision.

        +
      5. Linear Reservoir Parameter Set
        6. + 
        +[lrparamset rundu]
+gauge=rundu
+coem=1611.115479
+river=307.980042
+under=2531.556641
+leako=0.918236
+leaki=0.017568
+th=8.140809
+iso=0.000040
+isu=0.000073
+
        + To be completed in a future revision

        +
      6. Kinematic Wave Parameter Set
        7. + 
        +[KWParamSet rundu]
+GAUGE=rundu
+UNDER=1.673110
+LEAKI=0.043105
+TH=6.658569
+ISU=0.000000
+ALPHA=2.991570
+BETA=0.932080
+ALPHA0=4.603945
+

        Example KWParamSet Block

        + GAUGE: This is the name of the gauge location block for which the parameters that follow it will be applied.
        + TH: The number of grid cells needed to flow into a cell for it to be declared a channel grid cell.
        + UNDER: The interflow flow speed multiplier
        + LEAKI: The amount of water leaked out of the interflow reservoir at each time step. Ranges between 0 and 1.
        + ISU: The initial value of the interflow reservoir.
        + ALPHA: The multiplier in Q = alpha*Abeta used for channel routing
        + BETA: The exponent in Q = alpha*Abeta used for channel routing
        + ALPHA0: The Alpha value used for overland routing.

        +
      7. Snow-17 Parameter Set
        8. + 
        +[snow17paramset tarbela]
+GAUGE=tarbela
+UADJ=0.184653
+MBASE=0.047224
+MFMAX=1.068658
+MFMIN=0.516059
+TIPM=0.911706
+NMF=0.077336
+PLWHC=0.093812
+SCF=2.219492
+

        Example Snow17ParamSet Block

        + GAUGE: This is the name of the gauge location block for which the parameters that follow it will be applied.
        + UADJ: The number of grid cells needed to flow into a cell for it to be declared a channel grid cell.
        + MBASE: The interflow flow speed multiplier
        + MFMAX: The amount of water leaked out of the interflow reservoir at each time step. Ranges between 0 and 1.
        + MFMIN: The initial value of the interflow reservoir.
        + TIPM: The multiplier in Q = alpha*Abeta used for channel routing
        + NMF: The exponent in Q = alpha*Abeta used for channel routing
        + PLWHC: The Alpha value used for overland routing.
        + SCF: The Alpha value used for overland routing.

        +
      8. Simple Inundation Parameter Set
        9. + 
        +[simpleinundationparamset rundu]
+gauge=rundu
+alpha=2.991570
+beta=0.932080
+

        Example SimpleInundationParamSet Block

        + GAUGE: This is the name of the gauge location block for which the parameters that follow it will be applied.
        + ALPHA: The multiplier in A = (Q/alpha)1/beta used for computing cross-sectional area
        + BETA: The exponent in A = (Q/alpha)1/beta used for computing cross-sectional area

        +
      +
    7. Tasks
      8. +

      Tasks define the necessary information about which model to run, over what time period to run, with what time step to run, etc.

      +
      [Task RunFF]
+STYLE=SIMU
+MODEL=CREST
+BASIN=FF
+PRECIP=Q2_PRECIP
+PET=PET
+OUTPUT=/EF5Demo/FF/output/
+PARAM_SET=FF
+TIMESTEP=5u
+TIME_BEGIN=201006010000
+TIME_END=201006010030
+

      Example Task Block

      + + STYLE: The style of task that this is. Possible values are:
      + 
      +SIMU: A simulation run.
+SIMU_RP: A simulation run that will produce standard deviation, average and skew coefficient grids for estimating return period.
+CALI_DREAM: A calibration run using DREAM.
+CLIP_BASIN: Clips the basic files to the specified BASIN.
+CLIP_GAUGE: Clips the basic files to the first specified gauge.
      + MODEL: The water balance model that this task should use. Possible values are:
      + 
      +CREST: The Coupled Routing and Excess STorage distributed hydrological model.
+CRESTPHYS: CREST model + baseflow bucket
+SAC: The SAC-SMA model.
+HyMOD: HyMod.
      + ROUTING: The routing method that this task should use. Possible values are:
      + 
      +KW: Kinematic Wave Routing.
+LR: Linear Reservoir Routing.
      + SNOW: (Optional) The snow melt model that this task should use. Possible values are:
      + 
      SNOW17: The Snow-17 snow melt model.
      + INUNDATION: (Optional) The inundation model that this task should use. Possible values are:
      + 
      SIMPLEINUNDATIOn: The Simple Inundation model.
      + BASIN: The basin block name which defines the area over which the model should be run.
      + DEFAULTPARAMSGAUGE: (Optional) The gauge for which parameters are specified and used for gauges which did not have parameters specified. Useful when modeling large areas with gridded parameters.
      + PRECIP: The precipitation block name which defines the precipitation.
      + PRECIPFORECAST: (Optional) The precipitation block name which defines the forecast precipitation. Used if the precipiation specified in PRECIP is unavailable for the current time.
      + TEMP: (Required if using SNOW) The temperature block name which defines the temperature data to be used.
      + TEMPFORECAST: (Optional if using SNOW) The temperature block name which defines the forecast temperatures. Used if the temperature specified in TEMP is unavailable for the current time.
      + PET: The PET block name which defines the PET to use.
      + PARAM_SET: The parameter set block name which defines which set of water balance parameters to use.
      + ROUTING_PARAM_SET: The parameter set block name which defines which set of routing parameters to use.
      + SNOW_PARAM_SET: (Required if using SNOW) The parameter set block name which defines which set of snow parameters to use.
      + INUNDATION_PARAM_SET: (Required if using INUNDATION) The parameter set block name which defines which set of inundation parameters to use.
      + CALI_PARAM: (Required if using CALI_DREAM) The parameter set block name which defines which set of water balance parameters settings for calibration.
      + ROUTING_CALI_PARAM: (Required if using CALI_DREAM) The parameter set block name which defines which set of routing parameters to use for calibration.
      + SNOW_CALI_PARAM: (Required if using SNOW, CALI_DREAM) The parameter set block name which defines which set of snow parameters to use for calibration.
      + INUNDATION_CALI_PARAM: (Required if using INUNDATION, CALI_DREAM) The parameter set block name which defines which set of inundation parameters to use for calibration.
      + PRELOAD_FILE: (Optional) The file path and name where for the preload file. The preload file contains the forcings (Precip, PET, Temp) defined for the current time period and basin extent. Generated by EF5 if it does not exist. Useful for faster runs when forcings are not changing such as with manual calibration.
      + STATES: (Optional) The location where output files should be written.
      + TIMESTEP: The time step to use when running the model. Supported time units are year (y), month (m), day (d), hour (h), minute (u) and second (s).
      + TIME_BEGIN: The initialization time for the model run. YYYYMMDDHHUUSS format.
      + TIME_END: The ending time for the model run. YYYYMMDDHHUUSS format.
      + TIME_WARMEND: The end of the warm-up period. YYYYMMDDHHUUSS format.
      + TIME_STATE: (Optional) The time at which to output model states. YYYYMMDDHHUUSS format.
      + OUTPUT: The location where output files should be written.
      + OUTPUT_GRIDS: (Optional) Which grids should be output, combine together with | :
      + 
      +NONE: No gridded output.
+streamflow: Streamflow (cms)
+unitstreamflow: unit discharge (q/fac)
+soilmoisture: Soil moisture (%)
+
+returnperiod: Streamflow return period (years)
+precip: Precipitaiton (mm)
+pet: Potential evapotranspiration (mm)
+snowwater: Snow water equivalent from the snow melt model (mm)
+temperature: Temperature (C)
+inundation: Water depth (m)
+maxstreamflow: Maximum streamflow during run (cms)
+maxsoilmoisture: Maximum soil moisture during run (%)
+maxreturnperiod: Maximum return period during run (years)
+maxsnowwater: Maximum snow water equivalent (mm)
+runoff: Surface rain  (mm/hr)
+groundwater: Conceptual groundwater table (mm)
+subrunoff: Subsurface runoff (mm/hr)
      + DA_FILE: (Optional) The input observations to be added for use with streamflow data assimilation,
      + CO_FILE: (Optional) The combined output file if one is desired.
      + RP_STDGRID: (Optional) The geotiff file representing the standard deviation for the Log-Pearson Type III return period distribution.
      + RP_AVGGRID: (Optional) The geotiff file representing the average for the Log-Pearson Type III return period distribution.
      + RP_CSGRID: (Optional) The geotiff file representing the skew coefficient for the Log-Pearson Type III return period distribution.
      +
    8. Execute Block
      9. +

      The execute block defines which tasks to run when the ef5 program is executed

      +
      [Execute]
+TASK=RunFF
+

      Example Execute Block

      + TASK: The task block name which should be executed.

      +
    +
  8. Running EF5
    9. + 
    ef5 [controlfile]
    +

    Running ef5 is straight forwarded, you can called the executable with no arguments and it will assume a default configuration file name of "control.txt" in the current working directory or you can pass in a configuration file name as an argument.

    +
  9. Calibrating the Models
    10. +

    This section will be filled in once calibration methods are implemented.

    +
  10. Appendix
    11. +
      +
    1. Complete Sample Configuration File
      2. + 
      +/*
+ * This is an example configuration file for EF5
+ */
+
+[Basic]
+DEM=/EF5Demo/FF/basic/DEM.asc
+DDM=/EF5Demo/FF/basic/DDM.asc
+FAM=/EF5Demo/FF/basic/FAM.asc
+PROJ=laea
+ESRIDDM=true
+
+[PrecipForcing Q2Precip]
+TYPE=BIF
+UNIT=mm/h
+FREQ=5u
+LOC=/EF5Demo/FF/precip
+NAME=Q2_YYYYMMDDHHUU.bif
+
+[PETForcing PET]
+TYPE=BIF
+UNIT=mm/h
+FREQ=m
+LOC=/EF5Demo/FF/pet
+NAME=PET_MM.bif
+
+[Gauge OKC]
+LON=-97.01
+LAT=35.68
+OBS=/EF5Demo/obs/okc.csv
+
+[Gauge AR]
+LON=-93.62
+LAT=34.37
+
+[Basin FF]
+GAUGE=OKC
+GAUGE=AR
+
+[CrestParamSet FF]
+GAUGE=AR
+COEM=24.230076 EXPM=0.502391 RIVER=1.73056
+UNDER=0.291339 LEAKO=0.56668 LEAKI=0.251648
+TH=63.20205 GM=1.364364 PWM=71.96465
+PB=0.964355 PIM=6.508687 PKE=0.19952
+PFC=2.578529 IWU=53.52593 ISO=5.899539
+ISU=17.31128
+GAUGE=OKC
+COEM=24.230076 EXPM=0.502391 RIVER=1.73056
+UNDER=0.291339 LEAKO=0.56668 LEAKI=0.251648
+TH=63.20205 GM=1.364364 PWM=71.96465
+PB=0.964355 PIM=6.508687 PKE=0.19952
+PFC=2.578529 IWU=53.52593 ISO=5.899539
+ISU=17.31128
+
+[Task RunFF]
+STYLE=SIMU
+MODEL=CREST
+BASIN=FF
+PRECIP=Q2_PRECIP
+PET=PET
+OUTPUT=/EF5Demo/FF/output/
+PARAM_SET=FF
+TIMESTEP=5u
+TIME_BEGIN=201006010000
+TIME_END=201006010030
+
+[Execute]
+TASK=RunFF

      Example of a full EF5 configuration file.

      +
    +
+
+ + \ No newline at end of file diff --git a/install/Mac OS X user manual.docx b/install/Mac OS X user manual.docx new file mode 100644 index 0000000..9233cf5 Binary files /dev/null and b/install/Mac OS X user manual.docx differ diff --git a/install/install.md b/install/install.md new file mode 100644 index 0000000..6b44ff9 --- /dev/null +++ b/install/install.md @@ -0,0 +1,3 @@ +# Mac OS install + +Please find the attachment. diff --git a/src/AscGrid.cpp b/src/AscGrid.cpp index a14727f..dfcbf88 100644 --- a/src/AscGrid.cpp +++ b/src/AscGrid.cpp @@ -205,7 +205,9 @@ void WriteFloatAscGrid(const char *file, FloatGrid *grid) { for (long row = 0; row < grid->numRows; row++) { long lastCol = grid->numCols - 1; for (long col = 0; col < grid->numCols; col++) { - fprintf(fileH, "%.05f%s", grid->data[row][col], + // 2019-04: output gridded surface runoff --------------------------------- + //fprintf(fileH, "%.05f%s", grid->data[row][col], + fprintf(fileH, "%.02f%s", grid->data[row][col], (col == lastCol) ? "\n" : " "); } } diff --git a/src/CRESTModel.cpp b/src/CRESTModel.cpp index b666938..88e53cd 100644 --- a/src/CRESTModel.cpp +++ b/src/CRESTModel.cpp @@ -10,6 +10,7 @@ static const char *stateStrings[] = { "SM", + "GW", }; CRESTModel::CRESTModel() {} @@ -100,22 +101,27 @@ void CRESTModel::SaveStates(TimeVar *currentTime, char *statePath, } } -bool CRESTModel::WaterBalance(float stepHours, std::vector *precip, +bool CRESTModel::WaterBalance(float stepHours, + std::vector *precip, std::vector *pet, std::vector *fastFlow, std::vector *slowFlow, - std::vector *soilMoisture) { + std::vector *baseFlow, + std::vector *soilMoisture, + std::vector *groundwater) { size_t numNodes = nodes->size(); #if _OPENMP //#pragma omp parallel for #endif + for (size_t i = 0; i < numNodes; i++) { GridNode *node = &nodes->at(i); CRESTGridNode *cNode = &(crestNodes[i]); WaterBalanceInt(node, cNode, stepHours, precip->at(i), pet->at(i), - &(fastFlow->at(i)), &(slowFlow->at(i))); + &(fastFlow->at(i)), &(slowFlow->at(i)), + &(baseFlow->at(i))); soilMoisture->at(i) = cNode->states[STATE_CREST_SM] * 100.0 / cNode->params[PARAM_CREST_WM]; } @@ -125,7 +131,7 @@ bool CRESTModel::WaterBalance(float stepHours, std::vector *precip, void CRESTModel::WaterBalanceInt(GridNode *node, CRESTGridNode *cNode, float stepHours, float precipIn, float petIn, - float *fastFlow, float *slowFlow) { + float *fastFlow, float *slowFlow, float *baseFlow) { double precip = precipIn * stepHours; // precipIn is mm/hr, precip is mm double pet = petIn * stepHours; // petIn in mm/hr, pet is mm double R = 0.0, Wo = 0.0; diff --git a/src/CRESTModel.h b/src/CRESTModel.h index b5d3420..1721c1a 100644 --- a/src/CRESTModel.h +++ b/src/CRESTModel.h @@ -8,6 +8,7 @@ enum STATES_CREST { STATE_CREST_SM, STATE_CREST_QTY }; enum CREST_LAYER { CREST_LAYER_OVERLAND, CREST_LAYER_INTERFLOW, + CREST_LAYER_BASEFLOW, CREST_LAYER_QTY, }; @@ -37,14 +38,16 @@ class CRESTModel : public WaterBalanceModel { bool WaterBalance(float stepHours, std::vector *precip, std::vector *pet, std::vector *fastFlow, std::vector *slowFlow, - std::vector *soilMoisture); + std::vector *baseFlow, + std::vector *soilMoisture, + std::vector *groundwater); bool IsLumped() { return false; } const char *GetName() { return "crest"; } private: void WaterBalanceInt(GridNode *node, CRESTGridNode *cNode, float stepHours, float precipIn, float petIn, float *fastFlow, - float *slowFlow); + float *slowFlow, float *baseFlow); void InitializeParameters(std::map *paramSettings, std::vector *paramGrids); diff --git a/src/CRESTPhysModel.cpp b/src/CRESTPhysModel.cpp new file mode 100644 index 0000000..121bb30 --- /dev/null +++ b/src/CRESTPhysModel.cpp @@ -0,0 +1,458 @@ +#include +#include +#include +#if _OPENMP +#include +#endif +#include "AscGrid.h" +#include "CRESTPhysModel.h" +#include "DatedName.h" + +static const char *stateStrings[] = { + "SM", + "GW" +}; + +CRESTPHYSModel::CRESTPHYSModel() {} + +CRESTPHYSModel::~CRESTPHYSModel() {} + +bool CRESTPHYSModel::InitializeModel( + std::vector *newNodes, + std::map *paramSettings, + std::vector *paramGrids) { + + nodes = newNodes; + if (crestphysNodes.size() != nodes->size()) { + crestphysNodes.resize(nodes->size()); + } + + // Fill in modelIndex in the gridNodes + size_t numNodes = nodes->size(); + for (size_t i = 0; i < numNodes; i++) { + GridNode *node = &nodes->at(i); + node->modelIndex = i; + } + + InitializeParameters(paramSettings, paramGrids); + + return true; +} + +void CRESTPHYSModel::InitializeStates(TimeVar *beginTime, char *statePath) { + DatedName timeStr; + timeStr.SetNameStr("YYYYMMDD_HHUU"); + timeStr.ProcessNameLoose(NULL); + timeStr.UpdateName(beginTime->GetTM()); + + char buffer[300]; + for (int p = 0; p < STATE_CRESTPHYS_QTY; p++) { + sprintf(buffer, "%s/crestphys_%s_%s.tif", statePath, stateStrings[p], + timeStr.GetName()); + + FloatGrid *sGrid = ReadFloatTifGrid(buffer); + if (sGrid) { + printf("Using CRESTPHYS %s State Grid %s\n", stateStrings[p], buffer); + if (g_DEM->IsSpatialMatch(sGrid)) { + for (size_t i = 0; i < nodes->size(); i++) { + GridNode *node = &nodes->at(i); + CRESTPHYSGridNode *cNode = &(crestphysNodes[i]); + if (sGrid->data[node->y][node->x] != sGrid->noData) { + cNode->states[p] = sGrid->data[node->y][node->x]; + } + } + } else { + GridLoc pt; + for (size_t i = 0; i < nodes->size(); i++) { + GridNode *node = &(nodes->at(i)); + CRESTPHYSGridNode *cNode = &(crestphysNodes[i]); + if (sGrid->GetGridLoc(node->refLoc.x, node->refLoc.y, &pt) && + sGrid->data[pt.y][pt.x] != sGrid->noData) { + cNode->states[p] = sGrid->data[pt.y][pt.x]; + } + } + } + delete sGrid; + } else { + printf("CRESTPHYS %s State Grid %s not found!\n", stateStrings[p], buffer); + } + } +} + +void CRESTPHYSModel::SaveStates(TimeVar *currentTime, char *statePath, + GridWriterFull *gridWriter) { + DatedName timeStr; + timeStr.SetNameStr("YYYYMMDD_HHUU"); + timeStr.ProcessNameLoose(NULL); + timeStr.UpdateName(currentTime->GetTM()); + + std::vector dataVals; + dataVals.resize(nodes->size()); + + char buffer[300]; + for (int p = 0; p < STATE_CRESTPHYS_QTY; p++) { + sprintf(buffer, "%s/crestphys_%s_%s.tif", statePath, stateStrings[p], + timeStr.GetName()); + for (size_t i = 0; i < nodes->size(); i++) { + CRESTPHYSGridNode *cNode = &(crestphysNodes[i]); + dataVals[i] = cNode->states[p]; + } + gridWriter->WriteGrid(nodes, &dataVals, buffer, false); + } +} + +bool CRESTPHYSModel::WaterBalance(float stepHours, + std::vector *precip, + std::vector *pet, + std::vector *fastFlow, + std::vector *interFlow, + std::vector *baseFlow, + std::vector *soilMoisture, + std::vector *groundwater) { + + size_t numNodes = nodes->size(); + +#if _OPENMP + //#pragma omp parallel for +#endif + for (size_t i = 0; i < numNodes; i++) { + GridNode *node = &nodes->at(i); + CRESTPHYSGridNode *cNode = &(crestphysNodes[i]); + WaterBalanceInt(node, + cNode, + stepHours, + precip->at(i), + pet->at(i), + &(fastFlow->at(i)), + &(interFlow->at(i)), + &(baseFlow->at(i))); + soilMoisture->at(i) = + cNode->states[STATE_CRESTPHYS_SM] * 100.0 / cNode->params[PARAM_CREST_WM]; + groundwater ->at(i) = cNode->states[STATE_CRESTPHYS_GW]; + } + + + return true; +} + +void CRESTPHYSModel::WaterBalanceInt(GridNode *node, CRESTPHYSGridNode *cNode, + float stepHours, float precipIn, float petIn, + float *fastFlow, float *interFlow, float *baseFlow) { + double precip = precipIn * stepHours; // precipIn is mm/hr, precip is mm + double pet = petIn * stepHours; // petIn in mm/hr, pet is mm + double R = 0.0, Wo = 0.0; + float baseflowExcess=0.0; + + double adjPET = pet * cNode->params[PARAM_CRESTPHYS_KE]; + double temX = 0.0; + + // If we aren't a channel cell, add routed in overland to precip + /*if (!node->channelGridCell) { + precip += *fastFlow; + *fastFlow = 0.0; + }*/ + + // Base flow continuity + cNode->states[STATE_CRESTPHYS_GW]+= *baseFlow; + + // We have more water coming in than leaving via ET. + if (precip > adjPET) { + double precipSoil = + (precip - adjPET) * (1 - cNode->params[PARAM_CRESTPHYS_IM]); // This is the + // precip that + // makes it to + // the soil + double precipImperv = + precip - adjPET - precipSoil; // Portion of precip on impervious area + + // cNode->states[STATE_CRESTPHYS_SM] += *slowFlow; + //*slowFlow = 0.0; + + double interflowExcess = + cNode->states[STATE_CRESTPHYS_SM] - cNode->params[PARAM_CRESTPHYS_WM]; + if (interflowExcess < 0.0) { + interflowExcess = 0.0; + } + + if (cNode->states[STATE_CRESTPHYS_SM] > cNode->params[PARAM_CRESTPHYS_WM]) { + cNode->states[STATE_CRESTPHYS_SM] = cNode->params[PARAM_CRESTPHYS_WM]; + } + + if (cNode->states[STATE_CRESTPHYS_SM] < cNode->params[PARAM_CRESTPHYS_WM]) { + double Wmaxm = + cNode->params[PARAM_CRESTPHYS_WM] * (1 + cNode->params[PARAM_CRESTPHYS_B]); + double A = Wmaxm * (1 - pow(1 - cNode->states[STATE_CRESTPHYS_SM] / + cNode->params[PARAM_CRESTPHYS_WM], + 1 / (1 + cNode->params[PARAM_CRESTPHYS_B]))); + if (precipSoil + A >= Wmaxm) { + R = precipSoil - + (cNode->params[PARAM_CRESTPHYS_WM] - + cNode->states[STATE_CRESTPHYS_SM]); // Leftovers after filling SM + + if (R < 0) { + printf("R to %f, [%f, %f, %f, %f, %f, %f]\n", R, + cNode->params[PARAM_CRESTPHYS_WM], cNode->params[PARAM_CRESTPHYS_B], + cNode->states[STATE_CRESTPHYS_SM], A, Wmaxm, precipSoil); + R = 0.0; + } + + Wo = cNode->params[PARAM_CRESTPHYS_WM]; + + } else { + double infiltration = + cNode->params[PARAM_CRESTPHYS_WM] * + (pow(1 - A / Wmaxm, 1 + cNode->params[PARAM_CRESTPHYS_B]) - + pow(1 - (A + precipSoil) / Wmaxm, + 1 + cNode->params[PARAM_CRESTPHYS_B])); + if (infiltration > precipSoil) { + infiltration = precipSoil; + } else if (infiltration < 0.0) { + printf("Infiltration went to %f, [%f, %f, %f, %f, %f, %f]\n", + infiltration, cNode->params[PARAM_CRESTPHYS_WM], + cNode->params[PARAM_CRESTPHYS_B], cNode->states[STATE_CRESTPHYS_SM], A, + Wmaxm, precipSoil); + } + + R = precipSoil - infiltration; + + if (R < 0) { + printf("R (infil) to %f, [%f, %f, %f, %f, %f, %f]\n", R, + cNode->params[PARAM_CRESTPHYS_WM], cNode->params[PARAM_CRESTPHYS_B], + cNode->states[STATE_CRESTPHYS_SM], A, Wmaxm, precipSoil); + R = 0.0; + } + Wo = cNode->states[STATE_CRESTPHYS_SM] + infiltration; + } + } else { + R = precipSoil; + Wo = cNode->params[PARAM_CRESTPHYS_WM]; + } + + // Now R is excess water, split it between overland & interflow + + temX = (cNode->states[STATE_CRESTPHYS_SM] + Wo) / + cNode->params[PARAM_CRESTPHYS_WM] / 2 * + (cNode->params[PARAM_CRESTPHYS_FC] * + stepHours); // Calculate how much water can infiltrate + + if (R <= temX) { + cNode->excess[CRESTPHYS_LAYER_BASEFLOW] = R; + } else { + cNode->excess[CRESTPHYS_LAYER_BASEFLOW] = temX; + } + cNode->excess[CRESTPHYS_LAYER_INTERFLOW] = + R - cNode->excess[CRESTPHYS_LAYER_BASEFLOW]; + + cNode->excess[CRESTPHYS_LAYER_OVERLAND]= precipImperv; + + cNode->actET = adjPET; + + cNode->excess[CRESTPHYS_LAYER_OVERLAND] += + interflowExcess; // Extra interflow that got routed. + } else { // All the incoming precip goes straight to ET + cNode->excess[CRESTPHYS_LAYER_OVERLAND] = 0.0; + cNode->excess[CRESTPHYS_LAYER_BASEFLOW] = 0.0; + cNode->excess[CRESTPHYS_LAYER_INTERFLOW] = 0.0; + + // cNode->states[STATE_CRESTPHYS_SM] += *slowFlow; + //*slowFlow = 0.0; + + double interflowExcess = + cNode->states[STATE_CRESTPHYS_SM] - cNode->params[PARAM_CRESTPHYS_WM]; + if (interflowExcess < 0.0) { + interflowExcess = 0.0; + } + cNode->excess[CRESTPHYS_LAYER_OVERLAND] = interflowExcess; + + if (cNode->states[STATE_CRESTPHYS_SM] > cNode->params[PARAM_CRESTPHYS_WM]) { + cNode->states[STATE_CRESTPHYS_SM] = cNode->params[PARAM_CRESTPHYS_WM]; + } + + double ExcessET = (adjPET - precip) * cNode->states[STATE_CRESTPHYS_SM] / + cNode->params[PARAM_CRESTPHYS_WM]; + if (ExcessET < cNode->states[STATE_CRESTPHYS_SM]) { + Wo = cNode->states[STATE_CRESTPHYS_SM] - + ExcessET; // We can evaporate away ExcessET too. + } else { + Wo = 0.0; // We don't have enough to evaporate ExcessET. + ExcessET = cNode->states[STATE_CRESTPHYS_SM]; + } + double resET = (adjPET - ExcessET - precip)*cNode->states[STATE_CRESTPHYS_GW]/cNode->params[PARAM_CRESTPHYS_HMAXAQ]; + + if (resETstates[STATE_CRESTPHYS_GW]){ + cNode->states[STATE_CRESTPHYS_GW]-= resET; + } + else { + resET= cNode->states[STATE_CRESTPHYS_GW]; + cNode->states[STATE_CRESTPHYS_GW]-= resET; + } + // printf("groundwater 0: %f\n", cNode->states[STATE_CRESTPHYS_GW]); + + // printf("groundwater 1: %f\n", cNode->states[STATE_CRESTPHYS_GW]); + cNode->actET = resET+ExcessET + precip; + } + + // Modified by Allen 02/14/23, we drain soil moisture based on a rate factor KSOIL + double soil2gw=Wo*cNode->params[PARAM_CRESTPHYS_KSOIL]; + if (Wo>=soil2gw){ + Wo-= soil2gw; + } + else{ + Wo = 0.0; + soil2gw=Wo; + } + + cNode->states[STATE_CRESTPHYS_SM] = Wo; + cNode->states[STATE_CRESTPHYS_GW] += (cNode->excess[CRESTPHYS_LAYER_BASEFLOW]+soil2gw); + // printf("groundwater 2: %f\n", cNode->states[STATE_CRESTPHYS_GW]); + // Add Overland Excess Water to fastFlow + *fastFlow += (cNode->excess[CRESTPHYS_LAYER_OVERLAND] / (stepHours * 3600.0f)); + + // Add Interflow Excess Water to slowFlow + *interFlow += (cNode->excess[CRESTPHYS_LAYER_INTERFLOW] / (stepHours * 3600.0f)); + + // Base flow we apply a fill-spill bucket + if (cNode->states[STATE_CRESTPHYS_GW]>cNode->params[PARAM_CRESTPHYS_HMAXAQ]) { + // spill + baseflowExcess= cNode->states[STATE_CRESTPHYS_GW]-cNode->params[PARAM_CRESTPHYS_HMAXAQ]; + cNode->states[STATE_CRESTPHYS_GW]=cNode->params[PARAM_CRESTPHYS_HMAXAQ]; + // printf("set groundwater level %.2f to maximum aquifer depth %.2f...\n",cNode->states[STATE_CRESTPHYS_GW], cNode->params[PARAM_CRESTPHYS_HMAXAQ]); + } + else{ + baseflowExcess=0.0; + } + double baseflowExp= cNode->params[PARAM_CRESTPHYS_GWC]*(exp(cNode->params[PARAM_CRESTPHYS_GWE]* + cNode->states[STATE_CRESTPHYS_GW]/cNode->params[PARAM_CRESTPHYS_HMAXAQ])-1); + // printf("GWC: %f, GWE: %f, GW: %f, HMAX: %f, baseflow: %f\n", cNode->params[PARAM_CRESTPHYS_GWC], + // cNode->params[PARAM_CRESTPHYS_GWE], cNode->states[STATE_CRESTPHYS_GW],cNode->params[PARAM_CRESTPHYS_HMAXAQ], baseflowExp); + // printf("groundwater 3: %f\n", cNode->states[STATE_CRESTPHYS_GW]); + if (baseflowExp>cNode->states[STATE_CRESTPHYS_GW]){ + baseflowExp= cNode->states[STATE_CRESTPHYS_GW]; + cNode->states[STATE_CRESTPHYS_GW]=0; + } + else{ + cNode->states[STATE_CRESTPHYS_GW]-=baseflowExp; + } + // printf("groundwater flow : %f\n", baseflowExp); + // printf("groundwater 4: %f\n", cNode->states[STATE_CRESTPHYS_GW]); + *baseFlow = ((baseflowExcess+baseflowExp) / (stepHours * 3600.0f)); + // printf("baseflow: %f\n", *baseFlow); + + // Calculate Discharge as the sum of the leaks + //*discharge = (overlandLeak + interflowLeak) * node->area / 3.6; +} + +void CRESTPHYSModel::InitializeParameters( + std::map *paramSettings, + std::vector *paramGrids) { + + // This pass distributes parameters + size_t numNodes = nodes->size(); + size_t unused = 0; + for (size_t i = 0; i < numNodes; i++) { + GridNode *node = &nodes->at(i); + CRESTPHYSGridNode *cNode = &(crestphysNodes[i]); + if (!node->gauge) { + unused++; + continue; + } + /*if (i == 0) { + float *paramsBlah = (*paramSettings)[node->gauge]; + printf("Thread %i, %f\n", omp_get_thread_num(), + paramsBlah[PARAM_CREST_QTY-1]); + }*/ + // Copy all of the parameters over + memcpy(cNode->params, (*paramSettings)[node->gauge], + sizeof(float) * PARAM_CRESTPHYS_QTY); + + // Some of the parameters are special, deal with that here + if (!paramGrids->at(PARAM_CRESTPHYS_IM)) { + cNode->params[PARAM_CRESTPHYS_IM] /= 100.0; + } + + // Deal with the distributed parameters here + GridLoc pt; + for (size_t paramI = 0; paramI < PARAM_CRESTPHYS_QTY; paramI++) { + FloatGrid *grid = paramGrids->at(paramI); + if (grid && g_DEM->IsSpatialMatch(grid)) { + if (grid->data[node->y][node->x] == 0) { + grid->data[node->y][node->x] = 0.01; + } + cNode->params[paramI] *= grid->data[node->y][node->x]; + } else if (grid && + grid->GetGridLoc(node->refLoc.x, node->refLoc.y, &pt)) { + if (grid->data[pt.y][pt.x] == 0) { + grid->data[pt.y][pt.x] = 0.01; + // printf("Using nodata value in param %s\n", + // modelParamStrings[MODEL_CREST][paramI]); + } + cNode->params[paramI] *= grid->data[pt.y][pt.x]; + } + } + + if (!paramGrids->at(PARAM_CRESTPHYS_IWU)) { + cNode->states[STATE_CRESTPHYS_SM] = cNode->params[PARAM_CRESTPHYS_IWU] * + cNode->params[PARAM_CRESTPHYS_WM] / 100.0; + } + + if (!paramGrids->at(PARAM_CRESTPHYS_IGW)) { + cNode->states[STATE_CRESTPHYS_GW] = cNode->params[PARAM_CRESTPHYS_IGW] * + cNode->params[PARAM_CRESTPHYS_HMAXAQ] / 100.0; + } + + if (cNode->params[PARAM_CRESTPHYS_WM] < 0.0) { + cNode->params[PARAM_CRESTPHYS_WM] = 100.0; + } + // printf("Maximum aquifer depth: %f\n",cNode->params[PARAM_CRESTPHYS_HMAXAQ]); + if (cNode->states[STATE_CRESTPHYS_SM] < 0.0) { + printf("Node Soil Moisture(%f) is less than 0, setting to 0.\n", + cNode->states[STATE_CRESTPHYS_SM]); + cNode->states[STATE_CRESTPHYS_SM] = 0.0; + } else if (cNode->states[STATE_CRESTPHYS_SM] > cNode->params[PARAM_CREST_WM]) { + printf("Node Soil Moisture(%f) is greater than WM, setting to %f.\n", + cNode->states[STATE_CRESTPHYS_SM], cNode->params[PARAM_CREST_WM]); + } + + if (cNode->params[PARAM_CRESTPHYS_IM] < 0.0) { + // printf("Node Impervious Area(%f) is less than 0, setting to 0.\n", + // cNode->params[PARAM_CREST_IM]); + cNode->params[PARAM_CRESTPHYS_IM] = 0.0; + } else if (cNode->params[PARAM_CRESTPHYS_IM] > 1.0) { + // printf("Node Impervious Area(%f) is greater than 1, setting to 1.\n", + // cNode->params[PARAM_CREST_IM]); + cNode->params[PARAM_CRESTPHYS_IM] = 1.0; + } + + if (cNode->params[PARAM_CRESTPHYS_B] < 0.0) { + // printf("Node B (%f) is less than 0, setting to 0.\n", + // cNode->params[PARAM_CREST_B]); + cNode->params[PARAM_CRESTPHYS_B] = 1.0; + } else if (cNode->params[PARAM_CRESTPHYS_B] != cNode->params[PARAM_CRESTPHYS_B]) { + // printf("Node B (%f) NaN, setting to %f.\n", + // cNode->params[PARAM_CREST_B], 0.0); + cNode->params[PARAM_CRESTPHYS_B] = 0.0; + } + + if (cNode->params[PARAM_CRESTPHYS_FC] < 0.0) { + // printf("Node B (%f) is less than 0, setting to 0.\n", + // cNode->params[PARAM_CREST_B]); + cNode->params[PARAM_CRESTPHYS_FC] = 1.0; + } + if (cNode->params[PARAM_CRESTPHYS_IGW]<0.0){ + cNode->params[PARAM_CRESTPHYS_IGW]=0.0; + } + + if (cNode->params[PARAM_CRESTPHYS_HMAXAQ]<0.0){ + cNode->params[PARAM_CRESTPHYS_HMAXAQ]=0.1; + } + + if (cNode->params[PARAM_CRESTPHYS_GWE]<0.0){ + cNode->params[PARAM_CRESTPHYS_GWE]=0; + } + + if (cNode->params[PARAM_CRESTPHYS_GWC]<0.0){ + cNode->params[PARAM_CRESTPHYS_GWC]=0; + } + + } +} diff --git a/src/CRESTPhysModel.h b/src/CRESTPhysModel.h new file mode 100644 index 0000000..1f424e1 --- /dev/null +++ b/src/CRESTPhysModel.h @@ -0,0 +1,64 @@ +#ifndef CRESTPHYS_MODEL_H +#define CRESTPHYS_MODEL_H + +#include "ModelBase.h" + +enum STATES_CRESTPHYS { STATE_CRESTPHYS_SM, + STATE_CRESTPHYS_GW, + STATE_CRESTPHYS_QTY }; + +enum CRESTPHYS_LAYER { + CRESTPHYS_LAYER_OVERLAND, + CRESTPHYS_LAYER_INTERFLOW, + CRESTPHYS_LAYER_BASEFLOW, + CRESTPHYS_LAYER_QTY, +}; + +struct CRESTPHYSGridNode : BasicGridNode { + float params[PARAM_CRESTPHYS_QTY]; + float states[STATE_CRESTPHYS_QTY]; + + // These guys are the state variables + // double soilMoisture; // Current depth of water filling available pore space + + // These are results of the water balance + double excess[CRESTPHYS_LAYER_QTY]; + double actET; +}; + +class CRESTPHYSModel : public WaterBalanceModel { + +public: + CRESTPHYSModel(); + ~CRESTPHYSModel(); + bool InitializeModel(std::vector *newNodes, + std::map *paramSettings, + std::vector *paramGrids); + void InitializeStates(TimeVar *beginTime, char *statePath); + void SaveStates(TimeVar *currentTime, char *statePath, + GridWriterFull *gridWriter); + bool WaterBalance(float stepHours, + std::vector *precip, + std::vector *pet, + std::vector *fastFlow, + std::vector *interFlow, + std::vector *baseFlow, + std::vector *soilMoisture, + std::vector *groundwater); + + bool IsLumped() { return false; } + const char *GetName() { return "crestphys"; } + +private: + void WaterBalanceInt(GridNode *node, CRESTPHYSGridNode *cNode, float stepHours, + float precipIn, float petIn, float *fastFlow, + float *interFlow, float *baseflow); + void + InitializeParameters(std::map *paramSettings, + std::vector *paramGrids); + + std::vector *nodes; + std::vector crestphysNodes; +}; + +#endif diff --git a/src/CaliParamConfigSection.cpp b/src/CaliParamConfigSection.cpp index f02021f..286a699 100644 --- a/src/CaliParamConfigSection.cpp +++ b/src/CaliParamConfigSection.cpp @@ -63,7 +63,7 @@ CONFIG_SEC_RET CaliParamConfigSection::ProcessKeyValue(char *name, } } ERROR_LOGF("Unknown objective function option \"%s\"!", value); - INFO_LOGF("Valid objective function options are \"%s\"", "NSCE, CC, SSE"); + INFO_LOGF("Valid objective function options are \"%s\"", "NSCE, LOGNSCE, CC, SSE, KGE"); return INVALID_RESULT; } else if (!strcasecmp(name, "ars_topnum")) { ars_topNum = atoi(value); diff --git a/src/DREAM.cpp b/src/DREAM.cpp index 334ff77..9220677 100644 --- a/src/DREAM.cpp +++ b/src/DREAM.cpp @@ -204,6 +204,7 @@ void DREAM::CalibrateParams() { MEMORYCHECK(pointerRUNvar->Table_JumpRate, "at dream.c: Memory Allocation for DREAM run variable " "Table_JumpRate not successfull\n"); + //------Step 1: Sample s points in the parameter //space---------------------------// if Extra.InitPopulation = 'LHS_BASED' // Latin hypercube sampling when indicated @@ -213,16 +214,20 @@ void DREAM::CalibrateParams() { // Step 2: Calculate posterior density associated with each value in x allocate2D(&p, pointerMCMC->seq, 2); + // printf(" 0 allocating memory ...\n"); MEMORYCHECK( p, "at dream.c: Memory Allocation for DREAM variable p not successfull\n"); log_p = (float *)malloc(pointerMCMC->seq * sizeof(float)); + // printf(" 1 allocating memory ...\n"); MEMORYCHECK(log_p, "at dream.c: Memory Allocation for DREAM variable log_p " "not successfull\n"); + // printf(" 2 allocating memory ...\n"); CompDensity(p, log_p, x, pointerMCMC, pointerInput, 3); - + // printf(" 3 allocating memory ...\n"); // Save the initial population, density and log density in one matrix X allocate2D(&X, pointerMCMC->seq, pointerMCMC->n + 2); + // printf(" 4 allocating memory ...\n"); MEMORYCHECK( X, "at dream.c: Memory Allocation for DREAM variable X not successfull\n"); @@ -233,9 +238,10 @@ void DREAM::CalibrateParams() { X[i][pointerInput->nPar] = p[i][0]; X[i][pointerInput->nPar + 1] = log_p[i]; } - + // printf(" 5 allocating memory ...\n"); // Then initialize the sequences // if save in memory = Yes + // printf("Start initializing sequences...\n"); InitSequences(X, pointerRUNvar->Sequences, pointerMCMC); // Reduced sample collection if reduced_sample_collection = Yes @@ -263,6 +269,7 @@ void DREAM::CalibrateParams() { Gelman(pointerOutput->R_stat, 0, pointerRUNvar->Sequences, pointerRUNvar->iloc, pointerMCMC->n, pointerMCMC->seq, 0); + // printf("preparing allocate variables \n"); // Allocate Variables allocate2D(&x_old, pointerMCMC->seq, pointerMCMC->n); MEMORYCHECK(x_old, "at dream.c: Memory Allocation for DREAM variable x_old " @@ -311,7 +318,7 @@ void DREAM::CalibrateParams() { #ifdef _WIN32 setIteration(0); #endif - + // printf("Start iteration...\n"); // Now start iteration ... while ((pointerRUNvar->Iter < pointerMCMC->ndraw) && !converged) { // Loop a number of times diff --git a/src/Defines.h b/src/Defines.h index de53852..53ea102 100644 --- a/src/Defines.h +++ b/src/Defines.h @@ -1,7 +1,7 @@ #ifndef DEFINES_H #define DEFINES_H -#define EF5_VERSION "1.2.3" +#define EF5_VERSION "1.3.3" #define CONFIG_MAX_LEN 256 diff --git a/src/EF5.cpp b/src/EF5.cpp index 705dc93..671ef14 100644 --- a/src/EF5.cpp +++ b/src/EF5.cpp @@ -55,7 +55,8 @@ int main(int argc, char *argv[]) { void PrintStartupMessage() { printf("%s", "********************************************************\n"); printf("%s", "** Ensemble Framework For Flash Flood Forecasting **\n"); - printf("** Version %s **\n", + printf("** Version %s **\n", EF5_VERSION); + printf("** Last modified date: 2023-02-14 **\n"); printf("%s", "********************************************************\n"); } diff --git a/src/GridWriterFull.cpp b/src/GridWriterFull.cpp index ac37342..f29d7f7 100644 --- a/src/GridWriterFull.cpp +++ b/src/GridWriterFull.cpp @@ -14,7 +14,10 @@ void GridWriterFull::Initialize() { grid.numCols = g_DEM->numCols; grid.numRows = g_DEM->numRows; grid.cellSize = g_DEM->cellSize; - grid.noData = -9999.0f; + // 2019-04: output gridded surface runoff --------------------------------- + //grid.noData = -9999.0f; + grid.noData = -9.0f; + // --------------------------------- grid.geoSet = g_DEM->geoSet; grid.modelType = g_DEM->modelType; grid.geographicType = g_DEM->geographicType; diff --git a/src/GriddedOutput.cpp b/src/GriddedOutput.cpp index 4a8420f..404f3b0 100644 --- a/src/GriddedOutput.cpp +++ b/src/GriddedOutput.cpp @@ -23,6 +23,9 @@ const char *GriddedOutputText[] = { "maxthresholdexceedance", "maxthresholdexceedancep", "precipaccum", + "runoff", + "groundwater", + "subrunoff", }; const int GriddedOutputFlags[] = { @@ -45,4 +48,7 @@ const int GriddedOutputFlags[] = { OG_MAXTHRES, OG_MAXTHRESP, OG_PRECIPACCUM, + OG_RUNOFF, + OG_GW, + OG_SUBSURF }; diff --git a/src/GriddedOutput.h b/src/GriddedOutput.h index 8621118..eb45d2c 100644 --- a/src/GriddedOutput.h +++ b/src/GriddedOutput.h @@ -21,9 +21,12 @@ enum SUPPORTED_OUTPUT_GRIDS { OG_MAXTHRES = 32768, OG_MAXTHRESP = 65536, OG_PRECIPACCUM = 131072, + OG_RUNOFF=262144, //surface runoff (overland) + OG_GW=524288, //Groundwater state + OG_SUBSURF=1048576, //subsurface runoff }; -#define OG_QTY 19 +#define OG_QTY 22 extern const char *GriddedOutputText[]; extern const int GriddedOutputFlags[]; diff --git a/src/HPModel.cpp b/src/HPModel.cpp index f89b072..a05d92d 100644 --- a/src/HPModel.cpp +++ b/src/HPModel.cpp @@ -39,7 +39,9 @@ bool HPModel::WaterBalance(float stepHours, std::vector *precip, std::vector *pet, std::vector *fastFlow, std::vector *slowFlow, - std::vector *soilMoisture) { + std::vector *baseFlow, + std::vector *soilMoisture, + std::vector *groundwater) { size_t numNodes = nodes->size(); diff --git a/src/HPModel.h b/src/HPModel.h index b14d48a..96e88d9 100644 --- a/src/HPModel.h +++ b/src/HPModel.h @@ -21,7 +21,9 @@ class HPModel : public WaterBalanceModel { bool WaterBalance(float stepHours, std::vector *precip, std::vector *pet, std::vector *fastFlow, std::vector *slowFlow, - std::vector *soilMoisture); + std::vector *baseFlow, + std::vector *soilMoisture, + std::vector *groundwater); bool IsLumped() { return false; } const char *GetName() { return "hp"; } diff --git a/src/HyMOD.cpp b/src/HyMOD.cpp index aa09fb4..42fe54e 100644 --- a/src/HyMOD.cpp +++ b/src/HyMOD.cpp @@ -43,7 +43,9 @@ void HyMOD::SaveStates(TimeVar *currentTime, char *statePath, bool HyMOD::WaterBalance(float stepHours, std::vector *precip, std::vector *pet, std::vector *fastFlow, std::vector *slowFlow, - std::vector *soilMoisture) { + std::vector *baseFlow, + std::vector *soilMoisture, + std::vector *groundwater) { size_t numNodes = nodes->size(); size_t i = 0; diff --git a/src/HyMOD.h b/src/HyMOD.h index 1c93644..49f63df 100644 --- a/src/HyMOD.h +++ b/src/HyMOD.h @@ -36,7 +36,9 @@ class HyMOD : public WaterBalanceModel { bool WaterBalance(float stepHours, std::vector *precip, std::vector *pet, std::vector *fastFlow, std::vector *slowFlow, - std::vector *soilMoisture); + std::vector *baseFlow, + std::vector *soilMoisture, + std::vector *groundwater); bool IsLumped() { return true; } const char *GetName() { return "hymod"; } diff --git a/src/KinematicRoute.cpp b/src/KinematicRoute.cpp index 3c096d4..2950898 100644 --- a/src/KinematicRoute.cpp +++ b/src/KinematicRoute.cpp @@ -76,6 +76,7 @@ bool KWRoute::InitializeModel( KWGridNode *cNode = &(kwNodes[i]); cNode->slopeSqrt = pow(node->slope, 0.5f); cNode->hillSlopeSqrt = pow(node->slope * 0.5, 0.5f); + cNode->incomingWater[KW_LAYER_BASEFLOW]=0.0; cNode->incomingWater[KW_LAYER_INTERFLOW] = 0.0; cNode->incomingWater[KW_LAYER_FASTFLOW] = 0.0; cNode->incomingWaterOverland = 0.0; @@ -94,7 +95,8 @@ bool KWRoute::InitializeModel( void KWRoute::InitializeStates(TimeVar *beginTime, char *statePath, std::vector *fastFlow, - std::vector *slowFlow) { + std::vector *interFlow, + std::vector *baseFlow) { DatedName timeStr; timeStr.SetNameStr("YYYYMMDD_HHUU"); timeStr.ProcessNameLoose(NULL); @@ -159,7 +161,8 @@ void KWRoute::SaveStates(TimeVar *currentTime, char *statePath, } bool KWRoute::Route(float stepHours, std::vector *fastFlow, - std::vector *slowFlow, + std::vector *interFlow, + std::vector *baseFlow, std::vector *discharge) { if (!initialized) { @@ -172,28 +175,32 @@ bool KWRoute::Route(float stepHours, std::vector *fastFlow, for (long i = numNodes - 1; i >= 0; i--) { KWGridNode *cNode = &(kwNodes[i]); RouteInt(stepHours * 3600.0f, &(nodes->at(i)), cNode, fastFlow->at(i), - slowFlow->at(i)); + interFlow->at(i), baseFlow->at(i)); } for (size_t i = 0; i < numNodes; i++) { KWGridNode *cNode = &(kwNodes[i]); - slowFlow->at(i) = 0.0; // cNode->incomingWater[KW_LAYER_INTERFLOW]; + interFlow->at(i) = 0.0; // cNode->incomingWater[KW_LAYER_INTERFLOW]; + baseFlow ->at(i)= 0.0; // cNode->incomingWater[KW_LAYER_BASEFLOW]; fastFlow->at(i) = 0.0; // cNode->incomingWater[KW_LAYER_FASTFLOW]; cNode->incomingWaterOverland = 0.0; cNode->incomingWaterChannel = 0.0; if (!cNode->channelGridCell) { float q = cNode->incomingWater[KW_LAYER_FASTFLOW] * nodes->at(i).horLen; q += (cNode->incomingWater[KW_LAYER_INTERFLOW] * nodes->at(i).area / 3.6); + q += (cNode -> incomingWater[KW_LAYER_BASEFLOW]*nodes->at(i).area/3.6); discharge->at(i) = q; // * (stepHours * 3600.0f); } else { discharge->at(i) = cNode->incomingWater[KW_LAYER_FASTFLOW]; } cNode->states[STATE_KW_IR] = - cNode->states[STATE_KW_IR] + cNode->incomingWater[KW_LAYER_INTERFLOW]; + cNode->states[STATE_KW_IR] + cNode->incomingWater[KW_LAYER_INTERFLOW] + cNode->incomingWater[KW_LAYER_BASEFLOW]; cNode->incomingWater[KW_LAYER_INTERFLOW] = 0.0; // Zero here so we can save states cNode->incomingWater[KW_LAYER_FASTFLOW] = 0.0; // Zero here so we can save states + cNode->incomingWater[KW_LAYER_BASEFLOW] = + 0.0; // Zero here so we can save states } // InitializeRouting(stepHours * 3600.0f); @@ -202,9 +209,11 @@ bool KWRoute::Route(float stepHours, std::vector *fastFlow, } void KWRoute::RouteInt(float stepSeconds, GridNode *node, KWGridNode *cNode, - float fastFlow, float slowFlow) { + float fastFlow, float interFlow, float baseFlow) { if (!cNode->channelGridCell) { + // overland routing + // printf("%f\n",baseFlow); float beta = 0.6; float alpha = cNode->params[PARAM_KINEMATIC_ALPHA0]; @@ -264,7 +273,7 @@ void KWRoute::RouteInt(float stepSeconds, GridNode *node, KWGridNode *cNode, cNode->incomingWater[KW_LAYER_FASTFLOW] = newq; // Add Interflow Excess Water to Reservoir - cNode->states[STATE_KW_IR] += slowFlow; + cNode->states[STATE_KW_IR] += (interFlow+baseFlow); double interflowLeak = cNode->states[STATE_KW_IR] * cNode->params[PARAM_KINEMATIC_LEAKI]; // printf(" %f ", interflowLeak); @@ -282,33 +291,36 @@ void KWRoute::RouteInt(float stepSeconds, GridNode *node, KWGridNode *cNode, printf(" 0 got %f ", cNode->routeCNode[0][KW_LAYER_INTERFLOW]->incomingWater[KW_LAYER_INTERFLOW]); }*/ - cNode->routeCNode[0][KW_LAYER_INTERFLOW] - ->incomingWater[KW_LAYER_INTERFLOW] += leakAmount; + // printf("0 got here... incoming base flow: %f, leakAmount: %f\n",cNode->routeCNode[0][KW_LAYER_INTERFLOW]->incomingWater[KW_LAYER_INTERFLOW], leakAmount); + cNode->routeCNode[0][KW_LAYER_BASEFLOW]->incomingWater[KW_LAYER_BASEFLOW] += leakAmount; //*res += leakAmount; // Make this an atomic add for parallelization } + if (cNode->routeCNode[1][KW_LAYER_INTERFLOW]) { double interflowLeak1 = interflowLeak * cNode->routeAmount[1][KW_LAYER_INTERFLOW] * node->area / cNode->routeNode[1][KW_LAYER_INTERFLOW]->area; double leakAmount = interflowLeak1; // printf(" 1 got %f ", leakAmount); - double *res = &(cNode->routeCNode[1][KW_LAYER_INTERFLOW] - ->incomingWater[KW_LAYER_INTERFLOW]); + double *res = &(cNode->routeCNode[1][KW_LAYER_BASEFLOW]->incomingWater[KW_LAYER_BASEFLOW]); *res += leakAmount; // Make this an atomic add for parallelization } + } else { // First do overland routing float beta = 0.6; float alpha = cNode->params[PARAM_KINEMATIC_ALPHA0]; - slowFlow += cNode->incomingWater[KW_LAYER_INTERFLOW]; + interFlow += cNode->incomingWater[KW_LAYER_INTERFLOW]; // printf(" %f ", cNode->incomingWater[KW_LAYER_INTERFLOW]); fastFlow /= 1000.0; // mm to m - slowFlow /= 1000.0; // mm to m - float newInWater = (fastFlow + slowFlow); + baseFlow /= 1000.0; + interFlow /= 1000.0; // mm to m + float newInWater = (fastFlow + interFlow+ baseFlow); + // printf("fast flow: %f, interflow: %f, baseflow: %f\n",fastFlow,interFlow,baseFlow); float A, B, C, D, E; float backDiffq = 0.0; @@ -415,6 +427,7 @@ void KWRoute::RouteInt(float stepSeconds, GridNode *node, KWGridNode *cNode, cNode->incomingWater[KW_LAYER_FASTFLOW] = newWater; cNode->incomingWater[KW_LAYER_INTERFLOW] = 0.0; + cNode->incomingWater[KW_LAYER_BASEFLOW] = 0.0; } } @@ -440,6 +453,7 @@ void KWRoute::InitializeParameters( cNode->states[STATE_KW_IR] = cNode->params[PARAM_KINEMATIC_ISU]; } cNode->incomingWater[KW_LAYER_INTERFLOW] = 0.0; + cNode->incomingWater[KW_LAYER_BASEFLOW] = 0.0; cNode->incomingWater[KW_LAYER_FASTFLOW] = 0.0; // Deal with the distributed parameters here @@ -514,6 +528,7 @@ void KWRoute::InitializeRouting(float timeSeconds) { float nexTimeUnder = node->horLen / speedUnder; cNode->nexTime[KW_LAYER_INTERFLOW] = nexTimeUnder; + cNode->nexTime[KW_LAYER_BASEFLOW] = nexTimeUnder; } // This pass figures out which cell water is routed to @@ -555,6 +570,12 @@ void KWRoute::InitializeRouting(float timeSeconds) { cNode->routeNode[1][KW_LAYER_INTERFLOW] = previousNode; cNode->routeCNode[1][KW_LAYER_INTERFLOW] = (previousNode) ? &(kwNodes[previousNode->modelIndex]) : NULL; + cNode->routeNode[0][KW_LAYER_BASEFLOW] = currentNode; + cNode->routeCNode[0][KW_LAYER_BASEFLOW] = + (currentNode) ? &(kwNodes[currentNode->modelIndex]) : NULL; + cNode->routeNode[1][KW_LAYER_BASEFLOW] = previousNode; + cNode->routeCNode[1][KW_LAYER_BASEFLOW] = + (previousNode) ? &(kwNodes[previousNode->modelIndex]) : NULL; if (currentNode && !kwNodes[currentNode->modelIndex].channelGridCell) { if ((currentSeconds - previousSeconds) > 0) { cNode->routeAmount[0][KW_LAYER_INTERFLOW] = @@ -562,10 +583,18 @@ void KWRoute::InitializeRouting(float timeSeconds) { (currentSeconds - previousSeconds); cNode->routeAmount[1][KW_LAYER_INTERFLOW] = 1.0 - cNode->routeAmount[0][KW_LAYER_INTERFLOW]; - } + + cNode->routeAmount[0][KW_LAYER_BASEFLOW] = + (timeSeconds - previousSeconds) / + (currentSeconds - previousSeconds); + cNode->routeAmount[1][KW_LAYER_BASEFLOW] = + 1.0 - cNode->routeAmount[0][KW_LAYER_BASEFLOW]; + } } else { cNode->routeAmount[0][KW_LAYER_INTERFLOW] = 1.0; cNode->routeAmount[1][KW_LAYER_INTERFLOW] = 0.0; + cNode->routeAmount[0][KW_LAYER_BASEFLOW] = 1.0; + cNode->routeAmount[1][KW_LAYER_BASEFLOW] = 0.0; } } } diff --git a/src/KinematicRoute.h b/src/KinematicRoute.h index 13b2d23..b4e6511 100644 --- a/src/KinematicRoute.h +++ b/src/KinematicRoute.h @@ -6,6 +6,7 @@ enum KW_LAYER { KW_LAYER_FASTFLOW, KW_LAYER_INTERFLOW, + KW_LAYER_BASEFLOW, KW_LAYER_QTY, }; @@ -45,16 +46,17 @@ class KWRoute : public RoutingModel { std::vector *paramGrids); void InitializeStates(TimeVar *beginTime, char *statePath, std::vector *fastFlow, - std::vector *slowFlow); + std::vector *interFlow, + std::vector *baseFlow); void SaveStates(TimeVar *currentTime, char *statePath, GridWriterFull *gridWriter); bool Route(float stepHours, std::vector *fastFlow, - std::vector *slowFlow, std::vector *discharge); + std::vector *interFlow, std::vector *baseFlow, std::vector *discharge); float GetMaxSpeed() { return maxSpeed; } private: void RouteInt(float stepSeconds, GridNode *node, KWGridNode *cNode, - float fastFlow, float slowFlow); + float fastFlow, float interFlow, float baseFlow); void InitializeParameters(std::map *paramSettings, std::vector *paramGrids); diff --git a/src/LinearRoute.cpp b/src/LinearRoute.cpp index f2005d0..7b8f46d 100644 --- a/src/LinearRoute.cpp +++ b/src/LinearRoute.cpp @@ -39,12 +39,14 @@ bool LRRoute::InitializeModel( void LRRoute::InitializeStates(TimeVar *beginTime, char *statePath, std::vector *fastFlow, - std::vector *slowFlow) {} + std::vector *interFlow, + std::vector *baseFlow) {} void LRRoute::SaveStates(TimeVar *currentTime, char *statePath, GridWriterFull *gridWriter) {} bool LRRoute::Route(float stepHours, std::vector *fastFlow, - std::vector *slowFlow, + std::vector *interFlow, + std::vector *baseFlow, std::vector *discharge) { if (!initialized) { @@ -56,14 +58,14 @@ bool LRRoute::Route(float stepHours, std::vector *fastFlow, for (size_t i = 0; i < numNodes; i++) { LRGridNode *cNode = &(lrNodes[i]); - RouteInt(&(nodes->at(i)), cNode, fastFlow->at(i), slowFlow->at(i)); + RouteInt(&(nodes->at(i)), cNode, fastFlow->at(i), interFlow->at(i), baseFlow->at(i)); } for (size_t i = 0; i < numNodes; i++) { LRGridNode *cNode = &(lrNodes[i]); fastFlow->at(i) = cNode->incomingWater[LR_LAYER_OVERLAND]; cNode->incomingWater[LR_LAYER_OVERLAND] = 0.0; - slowFlow->at(i) = cNode->incomingWater[LR_LAYER_INTERFLOW]; + interFlow->at(i) = cNode->incomingWater[LR_LAYER_INTERFLOW]; cNode->incomingWater[LR_LAYER_INTERFLOW] = 0.0; } @@ -73,7 +75,7 @@ bool LRRoute::Route(float stepHours, std::vector *fastFlow, } void LRRoute::RouteInt(GridNode *node, LRGridNode *cNode, float fastFlow, - float slowFlow) { + float interFlow, float baseFlow) { if (!node->channelGridCell) { cNode->reservoirs[LR_LAYER_OVERLAND] += fastFlow; @@ -91,7 +93,7 @@ void LRRoute::RouteInt(GridNode *node, LRGridNode *cNode, float fastFlow, } // Add Interflow Excess Water to Reservoir - cNode->reservoirs[LR_LAYER_INTERFLOW] += slowFlow; + cNode->reservoirs[LR_LAYER_INTERFLOW] += interFlow; double interflowLeak = cNode->reservoirs[LR_LAYER_INTERFLOW] * cNode->params[PARAM_LINEAR_LEAKI]; cNode->reservoirs[LR_LAYER_INTERFLOW] -= interflowLeak; diff --git a/src/LinearRoute.h b/src/LinearRoute.h index 409cbb8..f2308ef 100644 --- a/src/LinearRoute.h +++ b/src/LinearRoute.h @@ -35,16 +35,19 @@ class LRRoute : public RoutingModel { std::vector *paramGrids); void InitializeStates(TimeVar *beginTime, char *statePath, std::vector *fastFlow, - std::vector *slowFlow); + std::vector *baseFlow, + std::vector *interFlow); void SaveStates(TimeVar *currentTime, char *statePath, GridWriterFull *gridWriter); bool Route(float stepHours, std::vector *fastFlow, - std::vector *slowFlow, std::vector *discharge); + std::vector *interFlow, + std::vector *baseFlow, + std::vector *discharge); float GetMaxSpeed() { return maxSpeed; } private: void RouteInt(GridNode *node, LRGridNode *cNode, float fastFlow, - float slowFlow); + float interFlow, float baseFlow); void InitializeParameters(std::map *paramSettings, std::vector *paramGrids); diff --git a/src/Model.cpp b/src/Model.cpp index 92315fe..1d18f9a 100644 --- a/src/Model.cpp +++ b/src/Model.cpp @@ -48,6 +48,14 @@ const char *CREST[] = { #define ADDPARAMCREST(a, b) }; +const char *CRESTPHYS[] = { +#undef ADDPARAMCRESTPHYS +#define ADDPARAMCRESTPHYS(a, b) a, +#include "Models.tbl" +#undef ADDPARAMCRESTPHYS +#define ADDPARAMCRESTPHYS(a, b) +}; + const char *HYMOD[] = { #undef ADDPARAMHYMOD #define ADDPARAMHYMOD(a, b) a, @@ -69,6 +77,7 @@ const char *SAC[] = { const char **modelParamStrings[] = { paramStrings::HP, paramStrings::CREST, + paramStrings::CRESTPHYS, paramStrings::HYMOD, paramStrings::SAC, @@ -93,6 +102,14 @@ const char *CREST[] = { #define ADDPARAMCREST(a, b) }; +const char *CRESTPHYS[] = { +#undef ADDPARAMCRESTPHYS +#define ADDPARAMCRESTPHYS(a, b) a "_grid", +#include "Models.tbl" +#undef ADDPARAMCRESTPHYS +#define ADDPARAMCRESTPHYS(a, b) +}; + const char *HYMOD[] = { #undef ADDPARAMHYMOD #define ADDPARAMHYMOD(a, b) a "_grid", @@ -111,7 +128,10 @@ const char *SAC[] = { } // namespace paramGridStrings const char **modelParamGridStrings[] = { - paramGridStrings::HP, paramGridStrings::CREST, paramGridStrings::HYMOD, + paramGridStrings::HP, + paramGridStrings::CREST, + paramGridStrings::CRESTPHYS, + paramGridStrings::HYMOD, paramGridStrings::SAC}; const int numModelParams[] = { diff --git a/src/Model.h b/src/Model.h index b37e810..fbd4450 100644 --- a/src/Model.h +++ b/src/Model.h @@ -3,6 +3,7 @@ #define ADDMODEL(a, b) #define ADDPARAMCREST(a, b) +#define ADDPARAMCRESTPHYS(a, b) #define ADDPARAMHYMOD(a, b) #define ADDPARAMSAC(a, b) #define ADDPARAMHP(a, b) @@ -47,6 +48,16 @@ enum CREST_PARAMS { PARAM_CREST_QTY, }; +enum CRESTPHYS_PARAMS { +#undef ADDPARAMCRESTPHYS +#define ADDPARAMCRESTPHYS(a, b) PARAM_CRESTPHYS_##b, +#include "Models.tbl" +#undef ADDPARAMCRESTPHYS +#define ADDPARAMCRESTPHYS(a, b) + + PARAM_CRESTPHYS_QTY, +}; + enum HYMOD_PARAMS { #undef ADDPARAMHYMOD #define ADDPARAMHYMOD(a, b) PARAM_HYMOD_##b, diff --git a/src/ModelBase.h b/src/ModelBase.h index 232689f..19dec38 100644 --- a/src/ModelBase.h +++ b/src/ModelBase.h @@ -20,7 +20,8 @@ class Model { virtual bool RunTimeStep(long index, float stepHours, std::vector *precip, std::vector *pet, std::vector *discharge, - std::vector *soilMoisture) = 0; + std::vector *soilMoisture, + std::vector *groundwater) = 0; virtual bool IsLumped() = 0; virtual const char *GetName() = 0; }; @@ -35,11 +36,14 @@ class WaterBalanceModel { virtual void InitializeStates(TimeVar *beginTime, char *statePath) = 0; virtual void SaveStates(TimeVar *currentTime, char *statePath, GridWriterFull *gridWriter) = 0; - virtual bool WaterBalance(float stepHours, std::vector *precip, + virtual bool WaterBalance(float stepHours, + std::vector *precip, std::vector *pet, std::vector *fastFlow, - std::vector *slowFlow, - std::vector *soilMoisture) = 0; + std::vector *interFlow, + std::vector *baseFlow, + std::vector *soilMoisture, + std::vector *groundwater) = 0; virtual bool IsLumped() = 0; virtual const char *GetName() = 0; }; @@ -53,11 +57,13 @@ class RoutingModel { std::vector *paramGrids) = 0; virtual void InitializeStates(TimeVar *beginTime, char *statePath, std::vector *fastFlow, - std::vector *slowFlow) = 0; + std::vector *interFlow, + std::vector *baseFlow) = 0; virtual void SaveStates(TimeVar *currentTime, char *statePath, GridWriterFull *gridWriter) = 0; virtual bool Route(float stepHours, std::vector *fastFlow, - std::vector *slowFlow, + std::vector *interFlow, + std::vector *baseFlow, std::vector *discharge) = 0; virtual float GetMaxSpeed() = 0; virtual float SetObsInflow(long index, float inflow) = 0; diff --git a/src/Models.tbl b/src/Models.tbl index 75c3c07..ba80be3 100755 --- a/src/Models.tbl +++ b/src/Models.tbl @@ -10,6 +10,20 @@ ADDPARAMCREST("ke", KE) ADDPARAMCREST("fc", FC) ADDPARAMCREST("iwu", IWU) + +ADDMODEL("crestphys", CRESTPHYS) +ADDPARAMCRESTPHYS("wm", WM) +ADDPARAMCRESTPHYS("b", B) +ADDPARAMCRESTPHYS("im", IM) +ADDPARAMCRESTPHYS("ke", KE) +ADDPARAMCRESTPHYS("fc", FC) +ADDPARAMCRESTPHYS("ksoil", KSOIL) +ADDPARAMCRESTPHYS("iwu", IWU) +ADDPARAMCRESTPHYS("igw", IGW) +ADDPARAMCRESTPHYS("hmaxaq", HMAXAQ) +ADDPARAMCRESTPHYS("gwc", GWC) +ADDPARAMCRESTPHYS("gwe", GWE) + ADDMODEL("hymod", HYMOD) ADDPARAMHYMOD("huz", HUZ) ADDPARAMHYMOD("b", B) diff --git a/src/ObjectiveFunc.cpp b/src/ObjectiveFunc.cpp index 46ff3ad..8555356 100644 --- a/src/ObjectiveFunc.cpp +++ b/src/ObjectiveFunc.cpp @@ -3,18 +3,22 @@ #include #include -const char *objectiveStrings[] = {"nsce", "cc", "sse"}; +const char *objectiveStrings[] = {"nsce", "lognsce", "cc", "sse","kge"}; const OBJECTIVE_GOAL objectiveGoals[] = { + OBJECTIVE_GOAL_MAXIMIZE, OBJECTIVE_GOAL_MAXIMIZE, OBJECTIVE_GOAL_MAXIMIZE, OBJECTIVE_GOAL_MINIMIZE, + OBJECTIVE_GOAL_MAXIMIZE, }; // Internal functions for calculating actual objective function scores +static float CalcLOGNSCE(std::vector *obs, std::vector *sim); static float CalcNSCE(std::vector *obs, std::vector *sim); static float CalcCC(std::vector *obs, std::vector *sim); static float CalcSSE(std::vector *obs, std::vector *sim); +static float CalcKGE(std::vector *obs, std::vector *sim); // This is the main function for calculating objective functions, everything // passes through here first. @@ -23,11 +27,15 @@ float CalcObjFunc(std::vector *obs, std::vector *sim, switch (obj) { case OBJECTIVE_NSCE: - return CalcNSCE(obs, sim); + return CalcNSCE(obs,sim); + case OBJECTIVE_LOGNSCE: + return CalcLOGNSCE(obs, sim); case OBJECTIVE_CC: return CalcCC(obs, sim); case OBJECTIVE_SSE: return CalcSSE(obs, sim); + case OBJECTIVE_KGE: + return CalcKGE(obs, sim); default: return 0; } @@ -64,6 +72,41 @@ float CalcNSCE(std::vector *obs, std::vector *sim) { } } +float CalcLOGNSCE(std::vector *obs, std::vector *sim){ + float obsMean = 0, obsAcc = 0, simAcc = 0, validQs = 0; + size_t totalTimeSteps = obs->size(); + for (size_t tsIndex = 0; tsIndex < totalTimeSteps; tsIndex++) { + if ((*obs)[tsIndex] == (*obs)[tsIndex] && + (*sim)[tsIndex] == (*sim)[tsIndex]&& + (*sim)[tsIndex] !=0 && + (*obs)[tsIndex] !=0) { + obsMean += (*obs)[tsIndex]; + validQs++; + } + } + + obsMean /= validQs; + float logobsMean= log(obsMean); + + for (size_t tsIndex = 0; tsIndex < totalTimeSteps; tsIndex++) { + if ((*obs)[tsIndex] == (*obs)[tsIndex] && + (*sim)[tsIndex] == (*sim)[tsIndex] && + (*sim)[tsIndex] !=0 && + (*obs)[tsIndex] !=0) { + // printf("%f %f\n", (*obs)[tsIndex], (*sim)[tsIndex]); + obsAcc += powf( log((*obs)[tsIndex]) - logobsMean, 2.0); + simAcc += powf( log((*obs)[tsIndex]) - log((*sim)[tsIndex]), 2.0); + } + } + + float result = 1.0 - (simAcc / obsAcc); + if (result == result) { + return result; + } else { + return -10000000000.0; + } +} + float CalcCC(std::vector *obs, std::vector *sim) { float obsMean = 0, simMean = 0, obsAcc2 = 0, obsAcc = 0, simAcc = 0; @@ -109,3 +152,40 @@ float CalcSSE(std::vector *obs, std::vector *sim) { return sse; } +// TODO +float CalcKGE(std::vector *obs, std::vector *sim){ + float obsMean = 0,simMean = 0, obsAcc = 0, simAcc = 0, validQs = 0, simVar=0, obsAcc2=0; + size_t totalTimeSteps = obs->size(); + for (size_t tsIndex = 0; tsIndex < totalTimeSteps; tsIndex++) { + if ((*obs)[tsIndex] == (*obs)[tsIndex] && + (*sim)[tsIndex] == (*sim)[tsIndex]) { + obsMean += (*obs)[tsIndex]; + simMean += (*sim)[tsIndex]; + validQs++; + } + } + + obsMean /= validQs; + simMean /= validQs; + + for (size_t tsIndex = 0; tsIndex < totalTimeSteps; tsIndex++) { + if ((*obs)[tsIndex] == (*obs)[tsIndex] && + (*sim)[tsIndex] == (*sim)[tsIndex]) { + // printf("%f %f\n", (*obs)[tsIndex], (*sim)[tsIndex]); + obsAcc += powf((*obs)[tsIndex] - obsMean, 2.0); + obsAcc2 += (((*obs)[tsIndex] - obsMean) * ((*sim)[tsIndex] - simMean)); + simAcc += powf((*obs)[tsIndex] - (*sim)[tsIndex], 2.0); + simVar += powf((*sim)[tsIndex] - simMean, 2.0); + } + } + float r= obsAcc2 / (sqrt(obsAcc) * sqrt(simAcc)); + float beta= simMean/obsMean; + float gamma= sqrt(simVar)/simMean/(sqrt(obsAcc/obsMean)); + float result= 1 - sqrt(powf(r-1,2)+powf(beta-1,2)+powf(gamma-1,2)); + + if (result == result) { + return result; + } else { + return -10000000000.0; + } +} diff --git a/src/ObjectiveFunc.h b/src/ObjectiveFunc.h index cb21009..1776318 100644 --- a/src/ObjectiveFunc.h +++ b/src/ObjectiveFunc.h @@ -5,8 +5,10 @@ enum OBJECTIVES { OBJECTIVE_NSCE, + OBJECTIVE_LOGNSCE, OBJECTIVE_CC, OBJECTIVE_SSE, + OBJECTIVE_KGE, OBJECTIVE_QTY, }; diff --git a/src/ParamSetConfigSection.cpp b/src/ParamSetConfigSection.cpp index e2b49cc..7f10924 100644 --- a/src/ParamSetConfigSection.cpp +++ b/src/ParamSetConfigSection.cpp @@ -77,6 +77,7 @@ CONFIG_SEC_RET ParamSetConfigSection::ProcessKeyValue(char *name, char *value) { for (int i = 0; i < numParams; i++) { // printf("%i %i %s %s\n", model, i, modelParamStrings[model][i], name); if (strcasecmp(name, modelParamStrings[model][i]) == 0) { + INFO_LOGF("readed parameter \"%s\" !", name); if (currentParamsSet[i]) { ERROR_LOGF("Duplicate parameter \"%s\" in parameter set!", name); @@ -89,7 +90,6 @@ CONFIG_SEC_RET ParamSetConfigSection::ProcessKeyValue(char *name, char *value) { return VALID_RESULT; } } - // We got here so we must not know what this parameter is! ERROR_LOGF("Unknown parameter name \"%s\".", name); return INVALID_RESULT; diff --git a/src/SAC.cpp b/src/SAC.cpp index c2dc50a..37a887a 100644 --- a/src/SAC.cpp +++ b/src/SAC.cpp @@ -221,7 +221,9 @@ void SAC::SaveStates(TimeVar *currentTime, char *statePath, bool SAC::WaterBalance(float stepHours, std::vector *precip, std::vector *pet, std::vector *fastFlow, std::vector *slowFlow, - std::vector *soilMoisture) { + std::vector *baseFlow, + std::vector *soilMoisture, + std::vector *groundwater) { size_t numNodes = nodes->size(); size_t i = 0; @@ -238,6 +240,7 @@ bool SAC::WaterBalance(float stepHours, std::vector *precip, WaterBalanceInt(node, cNode, stepHours, precip->at(i), pet->at(i)); fastFlow->at(i) += (cNode->dischargeF / (stepHours * 3600.0f)); slowFlow->at(i) += (cNode->dischargeS / (stepHours * 3600.0f)); + baseFlow->at(i) += (cNode->dischargeS / (stepHours * 3600.0f)); soilMoisture->at(i) = 100.0 * (cNode->UZTWC + cNode->UZFWC) / (cNode->params[PARAM_SAC_UZTWM] + cNode->params[PARAM_SAC_UZFWM]); diff --git a/src/SAC.h b/src/SAC.h index 9575808..1785be3 100644 --- a/src/SAC.h +++ b/src/SAC.h @@ -25,7 +25,9 @@ class SAC : public WaterBalanceModel { bool WaterBalance(float stepHours, std::vector *precip, std::vector *pet, std::vector *fastFlow, std::vector *slowFlow, - std::vector *soilMoisture); + std::vector *baseFlow, + std::vector *soilMoisture, + std::vector *groundwater); bool IsLumped() { return false; } const char *GetName() { return "sac"; } diff --git a/src/Simulator.cpp b/src/Simulator.cpp index a8dd137..83afc4d 100644 --- a/src/Simulator.cpp +++ b/src/Simulator.cpp @@ -6,6 +6,7 @@ #endif #include "BasicGrids.h" #include "CRESTModel.h" +#include "CRESTPhysModel.h" #include "GridWriterFull.h" #include "GriddedOutput.h" #include "HPModel.h" @@ -80,6 +81,16 @@ bool Simulator::InitializeBasic(TaskConfigSection *task) { timeStepTemp = NULL; } + if (timeStepPrecip->GetTimeInSec() < timeStep->GetTimeInSec()) { + ERROR_LOG("The time step for precipitation must be greater or equal to the overall time step."); + return false; + } + + if (timeStepPET->GetTimeInSec() < timeStep->GetTimeInSec()) { + ERROR_LOG("The time step for PET must be greater or equal to the overall time step."); + return false; + } + // Initialize unit converters precipConvert = (3600.0 / (float)task->GetPrecipSec()->GetUnitTime()->GetTimeInSec()); @@ -102,13 +113,16 @@ bool Simulator::InitializeBasic(TaskConfigSection *task) { // Initialize time information currentTime = *(task->GetTimeBegin()); + currentTimes= *(task->GetTimeBegins()); // Leave for multi-event calibration if necessary currentTimePrecip = *(task->GetTimeBegin()); currentTimeQPF = *(task->GetTimeBegin()); currentTimePET = *(task->GetTimeBegin()); currentTimeTemp = *(task->GetTimeBegin()); currentTimeTempF = *(task->GetTimeBegin()); beginTime = *(task->GetTimeBegin()); - endTime = *(task->GetTimeEnd()); + beginTimes= *(task->GetTimeBegins()); // Leave for multi-event calibration if necessary + endTime = *(task->GetTimeEnd()); // Leave for multi-event calibration if necessary + endTimes= *(task->GetTimeEnds()); warmEndTime = *(task->GetTimeWarmEnd()); // Initialize file name information @@ -209,12 +223,14 @@ bool Simulator::InitializeBasic(TaskConfigSection *task) { "(Invalid gauge location?)"); return false; } - // Create the appropriate model switch (task->GetModel()) { case MODEL_CREST: wbModel = new CRESTModel(); break; + case MODEL_CRESTPHYS: + wbModel = new CRESTPHYSModel(); + break; case MODEL_HYMOD: wbModel = new HyMOD(); break; @@ -314,8 +330,10 @@ bool Simulator::InitializeSimu(TaskConfigSection *task) { avgSWE.resize(gauges->size()); avgT.resize(gauges->size()); avgSM.resize(gauges->size()); + avgGW.resize(gauges->size()); avgFF.resize(gauges->size()); avgSF.resize(gauges->size()); + avgBF.resize(gauges->size()); // Initialize forcing vectors & discharge storage vector currentPrecipSimu.resize(nodes.size()); @@ -343,6 +361,7 @@ bool Simulator::InitializeSimu(TaskConfigSection *task) { } currentFF.resize(nodes.size()); currentSF.resize(nodes.size()); + currentBF.resize(nodes.size()); currentQ.resize(nodes.size()); currentSWE.resize(nodes.size()); currentDepth.resize(nodes.size()); @@ -350,6 +369,7 @@ bool Simulator::InitializeSimu(TaskConfigSection *task) { outputRP = false; currentFF.resize(lumpedNodes.size()); currentSF.resize(lumpedNodes.size()); + currentBF.resize(lumpedNodes.size()); currentQ.resize(lumpedNodes.size()); currentSWE.resize(lumpedNodes.size()); } @@ -367,8 +387,8 @@ bool Simulator::InitializeSimu(TaskConfigSection *task) { // setvbuf(gaugeOutputs[i], NULL, _IONBF, 0); fprintf(gaugeOutputs[i], "%s", "Time,Discharge(m^3 s^-1),Observed(m^3 s^-1),Precip(mm " - "h^-1),PET(mm h^-1),SM(%),Fast Flow(mm*1000),Slow " - "Flow(mm*1000)"); + "h^-1),PET(mm h^-1),SM(%),Groundwater (mm),Fast Flow(mm*1000),Slow " + "Flow(mm*1000),Base Flow(mm*1000)"); if (sModel) { fprintf(gaugeOutputs[i], "%s", ",Temperature (C),SWE(mm)"); } @@ -503,11 +523,13 @@ bool Simulator::InitializeCali(TaskConfigSection *task) { if (!wbModel->IsLumped()) { currentFF.resize(nodes.size()); currentSF.resize(nodes.size()); + currentBF.resize(nodes.size()); currentQ.resize(nodes.size()); currentSWE.resize(nodes.size()); } else { currentFF.resize(lumpedNodes.size()); currentSF.resize(lumpedNodes.size()); + currentBF.resize(lumpedNodes.size()); currentQ.resize(lumpedNodes.size()); currentSWE.resize(lumpedNodes.size()); } @@ -541,6 +563,9 @@ bool Simulator::InitializeCali(TaskConfigSection *task) { case MODEL_CREST: caliWBModels[i] = new CRESTModel(); break; + case MODEL_CRESTPHYS: + caliWBModels[i] = new CRESTPHYSModel(); + break; case MODEL_HYMOD: caliWBModels[i] = new HyMOD(); break; @@ -960,15 +985,15 @@ void Simulator::SaveTSOutput() { GaugeConfigSection *gauge = gauges->at(i); if (gaugeOutputs[i]) { if (std::isfinite(currentQ[gauge->GetGridNodeIndex()])) { - fprintf(gaugeOutputs[i], "%s,%.2f,%.2f,%.2f,%.2f,%.2f,%.4f,%.4f", + fprintf(gaugeOutputs[i], "%s,%.2f,%.2f,%.2f,%.2f,%.2f,%.2f,%.4f,%.4f,%.4f", currentTimeText.GetName(), currentQ[gauge->GetGridNodeIndex()], gauge->GetObserved(¤tTime), avgPrecip[i], avgPET[i], - avgSM[i], avgFF[i] * 1000.0, avgSF[i] * 1000.0); + avgSM[i], avgGW[i], avgFF[i] * 1000.0, avgSF[i] * 1000.0, avgBF[i]*1000.0); } else { - fprintf(gaugeOutputs[i], "%s,%.2f,nan,%.2f,%.2f,%.2f,%.4f,%.4f", + fprintf(gaugeOutputs[i], "%s,%.2f,nan,%.2f,%.2f,%.2f,%.2f,%.4f,%.4f,%.4f", currentTimeText.GetName(), gauge->GetObserved(¤tTime), avgPrecip[i], avgPET[i], - avgSM[i], avgFF[i] * 1000.0, avgSF[i] * 1000.0); + avgSM[i], avgGW[i], avgFF[i] * 1000.0, avgSF[i] * 1000.0, avgBF[i]*1000.0); } if (sModel) { fprintf(gaugeOutputs[i], ",%.2f,%.2f", avgT[i], avgSWE[i]); @@ -1189,6 +1214,12 @@ void Simulator::SimulateDistributed(bool trackPeaks) { int currentYear = -1; int indexYear = -1; + // 2019-04: output gridded surface runoff --------------------------------- + std::vector currentFF_o; + + currentFF_o.resize(currentFF.size()); + // --------------------------------- + // Initialize TempReader if (sModel) { tempReader.ReadDEM(tempSec->GetDEM()); @@ -1244,7 +1275,7 @@ void Simulator::SimulateDistributed(bool trackPeaks) { qpfAccum.resize(nodes.size()); savePrecip = true; } - + // printf("Got here before initializing water balance models..."); // Initialize our models // NORMAL_LOGF("%s\n", "Got here!4"); wbModel->InitializeModel(&nodes, &fullParamSettings, ¶mGrids); @@ -1267,7 +1298,7 @@ void Simulator::SimulateDistributed(bool trackPeaks) { if (useStates) { wbModel->InitializeStates(¤tTime, statePath); if (rModel) { - rModel->InitializeStates(¤tTime, statePath, ¤tFF, ¤tSF); + rModel->InitializeStates(¤tTime, statePath, ¤tFF, ¤tSF, ¤tBF); } if (sModel) { sModel->InitializeStates(¤tTime, statePath); @@ -1276,10 +1307,10 @@ void Simulator::SimulateDistributed(bool trackPeaks) { for (size_t i = 0; i < currentFF.size(); i++) { currentFF[i] = 0.0; currentSF[i] = 0.0; + currentBF[i]= 0.0; currentSWE[i] = 0.0; } } - // Set up stuff for peak tracking here, if needed if (trackPeaks) { numYears = GetNumSimulatedYears(); @@ -1295,12 +1326,14 @@ void Simulator::SimulateDistributed(bool trackPeaks) { // Hard coded RP counting std::vector count2, rpGrid, rpMaxGrid, maxGrid; std::vector SM; + std::vector GW; std::vector dailyMaxQ, dailyMinSM, dailyMaxQHour; count2.resize(currentFF.size()); rpGrid.resize(currentFF.size()); rpMaxGrid.resize(currentFF.size()); maxGrid.resize(currentFF.size()); SM.resize(currentFF.size()); + GW.resize(currentFF.size()); for (size_t i = 0; i < currentFF.size(); i++) { count2[i] = 0.0; @@ -1310,7 +1343,9 @@ void Simulator::SimulateDistributed(bool trackPeaks) { currentFF[i] = 0.0; currentSF[i] = 0.0; currentQ[i] = 0.0; + currentBF[i]=0.0; } + // printf("After initialization...\n"); #if _OPENMP double timeTotal = 0.0, timeCount = 0.0; @@ -1321,6 +1356,7 @@ void Simulator::SimulateDistributed(bool trackPeaks) { // Here we load the input forcings & actually run the model for (currentTime.Increment(timeStep); currentTime <= endTime; currentTime.Increment(timeStep)) { + #if _OPENMP #ifndef _WIN32 double beginTime = omp_get_wtime(); @@ -1354,21 +1390,48 @@ void Simulator::SimulateDistributed(bool trackPeaks) { currentPrecip = ¤tPrecipSnow; } } - + // printf("Got here before executing water balance...\n"); done // Integrate the models for this timestep if (!preloadedForcings) { wbModel->WaterBalance(stepHoursReal, currentPrecip, ¤tPETSimu, - ¤tFF, ¤tSF, &SM); + ¤tFF, ¤tSF, ¤tBF, &SM, &GW); } else { wbModel->WaterBalance(stepHoursReal, &(currentPrecipCali[tsIndex]), - &(currentPETCali[tsIndex]), ¤tFF, ¤tSF, - &SM); + &(currentPETCali[tsIndex]), ¤tFF, ¤tSF, ¤tBF, + &SM, &GW); + } + if (griddedOutputs && ((griddedOutputs & OG_RUNOFF)==OG_RUNOFF)) { + // 2021-04 Allen: output gridded surface runoff --------------------------------- + std::vector _runoff; + _runoff.resize(currentFF.size()); + for (size_t i = 0; i < currentFF.size(); i++) { + float val = currentFF[i] * 3600.0f; + _runoff[i] = val; + } + sprintf(buffer, "%s/runoff.%s.%s.tif", outputPath, + currentTimeTextOutput.GetName(), wbModel->GetName()); + gridWriter.WriteGrid(&nodes, &_runoff, buffer, false); } + + if (griddedOutputs && ((griddedOutputs & OG_SUBSURF) == OG_SUBSURF)) { + std::vector _subsurface; + _subsurface.resize(currentSF.size()); + for (size_t i = 0; i < currentSF.size(); i++) { + float val = currentSF[i] * 3600.0f; + _subsurface[i] = val; + } + sprintf(buffer, "%s/subrunoff.%s.%s.tif", outputPath, + currentTimeTextOutput.GetName(), wbModel->GetName()); + gridWriter.WriteGrid(&nodes, &_subsurface, buffer, false); + } + if (outputTS) { gaugeMap.GaugeAverage(&nodes, ¤tFF, &avgFF); gaugeMap.GaugeAverage(&nodes, ¤tSF, &avgSF); + gaugeMap.GaugeAverage(&nodes, ¤tBF, &avgBF); } + if (rModel && wantsDA) { AssimilateData(); } @@ -1378,7 +1441,8 @@ void Simulator::SimulateDistributed(bool trackPeaks) { double beginTimeR = omp_get_wtime(); #endif #endif - rModel->Route(stepHoursReal, ¤tFF, ¤tSF, ¤tQ); + rModel->Route(stepHoursReal, ¤tFF, ¤tSF, ¤tBF, ¤tQ); + // printf("After routing...\n"); #if _OPENMP #ifndef _WIN32 double endTimeR = omp_get_wtime(); @@ -1389,6 +1453,7 @@ void Simulator::SimulateDistributed(bool trackPeaks) { for (size_t i = 0; i < currentFF.size(); i++) { currentFF[i] = 0.0; currentSF[i] = 0.0; + currentBF[i] = 0.0; } } if (saveStates && stateTime == currentTime) { @@ -1424,6 +1489,7 @@ void Simulator::SimulateDistributed(bool trackPeaks) { if (outputTS) { gaugeMap.GaugeAverage(&nodes, &SM, &avgSM); + gaugeMap.GaugeAverage(&nodes, &GW, &avgGW); if (!preloadedForcings) { gaugeMap.GaugeAverage(&nodes, currentPrecip, &avgPrecip); gaugeMap.GaugeAverage(&nodes, ¤tPETSimu, &avgPET); @@ -1487,12 +1553,18 @@ void Simulator::SimulateDistributed(bool trackPeaks) { currentDepth[i] = val; } gridWriter.WriteGrid(&nodes, ¤tDepth, buffer, false); - } + } if ((griddedOutputs & OG_SM) == OG_SM) { sprintf(buffer, "%s/sm.%s.%s.tif", outputPath, currentTimeTextOutput.GetName(), wbModel->GetName()); gridWriter.WriteGrid(&nodes, &SM, buffer, false); } + if ((griddedOutputs & OG_GW) == OG_GW) { + sprintf(buffer, "%s/gw.%s.%s.tif", outputPath, + currentTimeTextOutput.GetName(), wbModel->GetName()); + gridWriter.WriteGrid(&nodes, &GW, buffer, false); + } + if (outputRP && ((griddedOutputs & OG_QRP) == OG_QRP)) { sprintf(buffer, "%s/rp.%s.%s.tif", outputPath, currentTimeTextOutput.GetName(), wbModel->GetName()); @@ -1544,6 +1616,7 @@ void Simulator::SimulateDistributed(bool trackPeaks) { currentTimeTextOutput.GetName(), wbModel->GetName()); gridWriter.WriteGrid(&nodes, ¤tDepth, buffer, false); } + } #if _OPENMP @@ -1707,6 +1780,7 @@ void Simulator::SimulateLumped() { size_t tsIndex = 0; std::vector SM; + std::vector GW; SM.resize(currentFF.size()); // Initialize our model @@ -1750,12 +1824,12 @@ void Simulator::SimulateLumped() { gaugeMap.GaugeAverage(&nodes, ¤tPrecipSimu, &avgPrecip); gaugeMap.GaugeAverage(&nodes, ¤tPETSimu, &avgPET); - wbModel->WaterBalance(timeStepHours, &avgPrecip, &avgPET, ¤tFF, - ¤tSF, &SM); + wbModel->WaterBalance(timeStepHours, &avgPrecip, &avgPET, ¤tFF, ¤tBF, + ¤tSF, &SM, &GW); } else { wbModel->WaterBalance(timeStepHours, &(currentPrecipCali[tsIndex]), - &(currentPETCali[tsIndex]), ¤tFF, ¤tSF, - &SM); + &(currentPETCali[tsIndex]), ¤tFF, ¤tSF, ¤tBF, + &SM, &GW); } // We only output after the warmup period is over if (warmEndTime <= currentTime) { @@ -2042,8 +2116,8 @@ float Simulator::SimulateForCali(float *testParams) { WaterBalanceModel *runModel; RoutingModel *runRoutingModel; SnowModel *runSnowModel; - std::vector currentFFCali, currentSFCali, currentQCali, simQCali, - SMCali, currentSWECali, currentPrecipSnow; + std::vector currentFFCali, currentSFCali, currentBFCali, currentQCali, simQCali, + SMCali, GWCali, currentSWECali, currentPrecipSnow; TimeVar currentTimeCali; std::map *currentWBParamSettings; std::map *currentRParamSettings; @@ -2095,10 +2169,12 @@ float Simulator::SimulateForCali(float *testParams) { currentFFCali.resize(currentFF.size()); currentSFCali.resize(currentFF.size()); + currentBFCali.resize(currentFF.size()); currentQCali.resize(currentFF.size()); currentSWECali.resize(currentFF.size()); currentPrecipSnow.resize(currentFF.size()); SMCali.resize(currentFF.size()); + GWCali.resize(currentFF.size()); simQCali.resize(simQ.size()); avgPrecip.resize(gauges->size()); avgPET.resize(gauges->size()); @@ -2126,10 +2202,10 @@ float Simulator::SimulateForCali(float *testParams) { gaugeMap.GaugeAverage(&nodes, petVec, &avgPET); printf("%f %f\n", precipVec->at(300), petVec->at(300)); }*/ - runModel->WaterBalance(timeStepHours, precipVec, petVec, ¤tFFCali, - ¤tSFCali, &SMCali); + runModel->WaterBalance(timeStepHours, precipVec, petVec, ¤tFFCali, ¤tSFCali, + ¤tBFCali, &SMCali, &GWCali); - runRoutingModel->Route(timeStepHours, ¤tFFCali, ¤tSFCali, + runRoutingModel->Route(timeStepHours, ¤tFFCali, ¤tSFCali, ¤tBFCali, ¤tQCali); if (warmEndTime <= currentTimeCali) { @@ -2157,7 +2233,7 @@ float Simulator::SimulateForCali(float *testParams) { float *Simulator::SimulateForCaliTS(float *testParams) { WaterBalanceModel *runModel; - std::vector currentFFCali, currentSFCali, SMCali; + std::vector currentFFCali, currentSFCali, currentBFCali, SMCali, GWCali; float *simQCali; TimeVar currentTimeCali; std::map *currentParamSettings; @@ -2197,8 +2273,8 @@ float *Simulator::SimulateForCaliTS(float *testParams) { std::vector *precipVec = &(currentPrecipCali[tsIndex]); std::vector *petVec = &(currentPETCali[tsIndex]); - runModel->WaterBalance(timeStepHours, precipVec, petVec, ¤tFFCali, - ¤tSFCali, &SMCali); + runModel->WaterBalance(timeStepHours, precipVec, petVec, ¤tFFCali, ¤tBFCali, + ¤tSFCali, &SMCali, &GWCali); if (warmEndTime <= currentTimeCali) { if (!runModel->IsLumped()) { simQCali[tsIndexWarm] = currentFFCali[caliGauge->GetGridNodeIndex()]; diff --git a/src/Simulator.h b/src/Simulator.h index 3006216..4d7ab1c 100644 --- a/src/Simulator.h +++ b/src/Simulator.h @@ -84,10 +84,11 @@ class Simulator { TimeVar currentTime, currentTimePrecip, currentTimeQPF, currentTimePET, currentTimeTemp, currentTimeTempF, beginTime, endTime, warmEndTime, beginLRTime; + std::vector currentTimes, beginTimes, endTimes; DatedName *precipFile, *qpfFile, *petFile, *tempFile, *tempFFile, currentTimeText, currentTimeTextOutput; - std::vector currentFF, currentSF, currentQ, avgPrecip, avgPET, avgSWE, - currentSWE, avgT, avgSM, avgFF, avgSF, currentDepth; + std::vector currentFF, currentSF, currentBF, currentQ, avgPrecip, avgPET, avgSWE, + currentSWE, avgT, avgSM,avgGW, avgFF, avgSF, avgBF, currentDepth; std::vector paramGrids, paramGridsRoute, paramGridsSnow, paramGridsInundation; bool hasQPF, hasTempF, wantsDA; diff --git a/src/TaskConfigSection.cpp b/src/TaskConfigSection.cpp index 329c102..8975ab7 100644 --- a/src/TaskConfigSection.cpp +++ b/src/TaskConfigSection.cpp @@ -100,12 +100,17 @@ char *TaskConfigSection::GetCOFile() { return coFile; } TimeVar *TaskConfigSection::GetTimeBegin() { return &timeBegin; } +std::vector *TaskConfigSection::GetTimeBegins() { return &timeBegins; } + +std::vector *TaskConfigSection::GetTimeEnds() { return &timeEnds; } + TimeVar *TaskConfigSection::GetTimeState() { return &timeState; } TimeVar *TaskConfigSection::GetTimeWarmEnd() { return &timeWarmEnd; } TimeVar *TaskConfigSection::GetTimeEnd() { return &timeEnd; } + TimeVar *TaskConfigSection::GetTimeBeginLR() { if (!timestepLRSet || !timeBeginLRSet) { return NULL; @@ -195,15 +200,18 @@ CONFIG_SEC_RET TaskConfigSection::ProcessKeyValue(char *name, char *value) { "CLIP_GAUGE, BASIN_AVG"); return INVALID_RESULT; } else if (!strcasecmp(name, "model")) { + for (int i = 0; i < MODEL_QTY; i++) { + if (!strcasecmp(value, modelStrings[i])) { + modelSet = true; model = (MODELS)i; return VALID_RESULT; } } ERROR_LOGF("Unknown water balance model option \"%s\"!", value); - INFO_LOGF("Valid water balance options are \"%s\"", "CREST, SAC, HP"); + INFO_LOGF("Valid water balance options are \"%s\"", "CREST, CRESTPHYS, SAC, HP"); return INVALID_RESULT; } else if (!strcasecmp(name, "routing")) { for (int i = 0; i < ROUTE_QTY; i++) { @@ -428,10 +436,20 @@ CONFIG_SEC_RET TaskConfigSection::ProcessKeyValue(char *name, char *value) { } timeStateSet = true; } else if (!strcasecmp(name, "time_begin")) { - if (!timeBegin.LoadTime(value)) { + + char *part = strtok(value, "|"); + while (part != NULL) { + + if (!timeBegin.LoadTime(part)) { ERROR_LOGF("Unknown time begin option \"%s\"", value); return INVALID_RESULT; + } + else{ + timeBegins.push_back(&timeBegin); + } + part = strtok(NULL, "|"); } + timeBeginSet = true; } else if (!strcasecmp(name, "time_begin_lr")) { if (!timeBeginLR.LoadTime(value)) { @@ -446,10 +464,22 @@ CONFIG_SEC_RET TaskConfigSection::ProcessKeyValue(char *name, char *value) { } timeWarmEndSet = true; } else if (!strcasecmp(name, "time_end")) { - if (!timeEnd.LoadTime(value)) { + char *part = strtok(value, "|"); + while (part != NULL) { + + if (!timeEnd.LoadTime(part)) { ERROR_LOGF("Unknown time end option \"%s\"", value); return INVALID_RESULT; + } + else{ + timeEnds.push_back(&timeEnd); + } + part = strtok(NULL, "|"); } + // if (!timeEnd.LoadTime(value)) { + // ERROR_LOGF("Unknown time end option \"%s\"", value); + // return INVALID_RESULT; + // } timeEndSet = true; } else if (!strcasecmp(name, "cali_param")) { if (!modelSet) { diff --git a/src/TaskConfigSection.h b/src/TaskConfigSection.h index 246e0cc..a2375e6 100644 --- a/src/TaskConfigSection.h +++ b/src/TaskConfigSection.h @@ -20,6 +20,7 @@ #include "TimeUnit.h" #include "TimeVar.h" #include +#include class TaskConfigSection : public ConfigSection { @@ -45,8 +46,10 @@ class TaskConfigSection : public ConfigSection { char *GetDAFile(); char *GetCOFile(); TimeVar *GetTimeBegin(); + std::vector *GetTimeBegins(); TimeVar *GetTimeWarmEnd(); TimeVar *GetTimeEnd(); + std::vector *GetTimeEnds(); TimeVar *GetTimeState(); TimeUnit *GetTimeStep(); TimeUnit *GetTimeStepLR(); @@ -76,6 +79,7 @@ class TaskConfigSection : public ConfigSection { CONFIG_SEC_RET ProcessKeyValue(char *name, char *value); CONFIG_SEC_RET ValidateSection(); int GetGriddedOutputs() { return griddedOutputs; } + static bool IsDuplicate(char *name); @@ -125,6 +129,8 @@ class TaskConfigSection : public ConfigSection { int griddedOutputs; bool LoadGriddedOutputs(char *value); + std::vector timeBegins; + std::vector timeEnds; }; extern std::map g_taskConfigs;