/*
 * MATLAB Compiler: 3.0
 * Date: Mon Sep 20 09:40:04 2004
 * Arguments: "-B" "macro_default" "-O" "all" "-O" "fold_scalar_mxarrays:on"
 * "-O" "fold_non_scalar_mxarrays:on" "-O" "optimize_integer_for_loops:on" "-O"
 * "array_indexing:on" "-O" "optimize_conditionals:on" "-B" "sgl" "-m" "-W"
 * "main" "-L" "C" "-t" "-T" "link:exe" "-h" "libmmfile.mlib" "-W" "mainhg"
 * "libmwsglm.mlib" "recarga_multiple" 
 */
#include "recarga.h"
#include "dlmread.h"
#include "iga.h"
#include "libmatlbm.h"
#include "libmmfile.h"
#include "relcond.h"
#include "scsii.h"
#include "storageout.h"
#include "stormcheck.h"
#include "thetah.h"

static mxChar _array1_[4] = { 'l', 'o', 'n', 'g' };
static mxArray * _mxarray0_;
static mxArray * _mxarray2_;
static mxArray * _mxarray3_;
static mxArray * _mxarray4_;
static mxArray * _mxarray5_;

static mxChar _array7_[12] = { 's', 'o', 'i', 'l', 's', 'D',
                               'N', 'R', '.', 't', 'x', 't' };
static mxArray * _mxarray6_;

static mxChar _array9_[2] = { 0x005c, 't' };
static mxArray * _mxarray8_;

static double _array11_[4] = { 0.0, 1.0, 10.0, 4.0 };
static mxArray * _mxarray10_;
static mxArray * _mxarray12_;
static mxArray * _mxarray13_;

static double _array15_[3] = { .01, .01, .01 };
static mxArray * _mxarray14_;

static double _array17_[3] = { 1e-06, 1e-06, 1e-06 };
static mxArray * _mxarray16_;
static mxArray * _mxarray18_;
static mxArray * _mxarray19_;
static mxArray * _mxarray20_;
static mxArray * _mxarray21_;
static mxArray * _mxarray22_;
static mxArray * _mxarray23_;
static mxArray * _mxarray24_;
static mxArray * _mxarray25_;

static double _array27_[3] = { 1.6, 1.378, 1.242 };
static mxArray * _mxarray26_;

static double _array29_[3] = { -.00024246, -4.7992e-07, -1.2314e-10 };
static mxArray * _mxarray28_;
static mxArray * _mxarray30_;
static mxArray * _mxarray31_;

static double _array33_[3] = { 1.6, 1.242, 1.242 };
static mxArray * _mxarray32_;

static double _array35_[3] = { -.00024246, -1.2314e-10, -1.2314e-10 };
static mxArray * _mxarray34_;
static mxArray * _mxarray36_;
static mxArray * _mxarray37_;
static mxArray * _mxarray38_;
static mxArray * _mxarray39_;
static mxArray * _mxarray40_;
static mxArray * _mxarray41_;
static mxArray * _mxarray42_;
static mxArray * _mxarray43_;
static mxArray * _mxarray44_;
static mxArray * _mxarray45_;
static mxArray * _mxarray46_;

static double _array48_[11] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                                0.0, 0.0, 0.0, 0.0, 0.0 };
static mxArray * _mxarray47_;

static double _array50_[2] = { 0.0, 0.0 };
static mxArray * _mxarray49_;

static mxChar _array52_[6] = { '*', 'c', 'u', 'b', 'i', 'c' };
static mxArray * _mxarray51_;

static mxChar _array54_[7] = { '*', 'l', 'i', 'n', 'e', 'a', 'r' };
static mxArray * _mxarray53_;
static mxArray * _mxarray55_;
static mxArray * _mxarray56_;
static mxArray * _mxarray57_;

static mxChar _array59_[4] = { '.', 't', 'x', 't' };
static mxArray * _mxarray58_;

static mxChar _array61_[1] = { 'w' };
static mxArray * _mxarray60_;

static mxChar _array63_[2] = { 'a', 't' };
static mxArray * _mxarray62_;

static mxChar _array65_[13] = { 'R', 'E', 'C', 'A', 'R', 'G', 'A',
                                ' ', '2', '.', 'X', 0x005c, 'n' };
static mxArray * _mxarray64_;

static mxChar _array67_[12] = { 'D', 'a', 't', 'e', ':', 0x005c,
                                't', ' ', '%', 's', 0x005c, 'n' };
static mxArray * _mxarray66_;

static mxChar _array69_[19] = { 'O', 'u', 't', 'p', 'u', 't', ' ',
                                'F', 'i', 'l', 'e', ':', 0x005c,
                                't', ' ', '%', 's', 0x005c, 'n' };
static mxArray * _mxarray68_;

static mxChar _array71_[32] = { 'N', 'u', 'm', 'b', 'e', 'r', ' ', 'o', 'f',
                                ' ', 't', 'i', 'm', 'e', ' ', 's', 't', 'e',
                                'p', 's', '=', 0x005c, 't', ' ', '%', '-',
                                '4', '.', '0', 'f', 0x005c, 'n' };
static mxArray * _mxarray70_;

static mxChar _array73_[35] = { 'C', 'P', 'U', ' ', 'e', 'l', 'a', 'p', 's',
                                'e', 'd', ' ', 't', 'i', 'm', 'e', ' ', '(',
                                's', ')', '=', 0x005c, 't', ' ', '%', '-', '1',
                                '0', '.', '0', 'f', 0x005c, 'n', 0x005c, 'n' };
static mxArray * _mxarray72_;

static mxChar _array75_[13] = { 'I', 'N', 'P', 'U', 'T', ' ', 'T',
                                'E', 'R', 'M', 'S', 0x005c, 'n' };
static mxArray * _mxarray74_;

static mxChar _array77_[20] = { '%', '%', 'I', 'm', 'p', 'e', 'r',
                                'v', 'i', 'o', 'u', 's', '=', 0x005c,
                                't', ' ', '%', 'g', 0x005c, 'n' };
static mxArray * _mxarray76_;

static mxChar _array79_[19] = { 'P', 'e', 'r', 'v', 'i', 'o', 'u',
                                's', ' ', 'C', 'N', '=', 0x005c,
                                't', ' ', '%', 'g', 0x005c, 'n' };
static mxArray * _mxarray78_;

static mxChar _array81_[26] = { 'F', 'a', 'c', 'i', 'l', 'i', 't', 'y', ' ',
                                'A', 'r', 'e', 'a', ' ', '(', 'm', '2', ')',
                                '=', 0x005c, 't', ' ', '%', 'g', 0x005c, 'n' };
static mxArray * _mxarray80_;

static mxChar _array83_[23] = { 'T', 'r', 'i', 'b', '.', ' ', 'A', 'r',
                                'e', 'a', ' ', '(', 'm', '2', ')', '=',
                                0x005c, 't', ' ', '%', 'g', 0x005c, 'n' };
static mxArray * _mxarray82_;

static mxChar _array85_[19] = { 'R', 'A', 'T', 'I', 'O', 'i', 'm',
                                'p', '2', 'r', 'g', '=', 0x005c,
                                't', ' ', '%', 'g', 0x005c, 'n' };
static mxArray * _mxarray84_;

static mxChar _array87_[22] = { 'R', 'A', 'T', 'I', 'O', 'p', 'e', 'r',
                                'v', '2', 'r', 'g', '=', 0x005c, 't', ' ',
                                '%', 'g', 0x005c, 'n', 0x005c, 'n' };
static mxArray * _mxarray86_;

static mxChar _array89_[59] = { 'L', 'A', 'Y', 'E', 'R', ' ', ' ', ' ', ' ',
                                ' ', ' ', ' ', 0x005c, 't', ' ', 'D', 'E', 'P',
                                'T', 'H', ' ', '(', 'c', 'm', ')', ' ', ' ',
                                0x005c, 't', ' ', 'K', 's', 'a', 't', ' ', '(',
                                'c', 'm', '/', 'h', 'r', ')', 0x005c, 't', ' ',
                                'T', 'E', 'X', 'T', 'U', 'R', 'E', ' ', ' ',
                                ' ', ' ', ' ', 0x005c, 'n' };
static mxArray * _mxarray88_;

static mxChar _array91_[60] = { 'D', 'e', 'p', 'r', 'e', 's', 's', 'i', 'o',
                                'n', ' ', ' ', 0x005c, 't', ' ', ' ', ' ', ' ',
                                ' ', '%', 'g', ' ', ' ', ' ', ' ', ' ', ' ',
                                0x005c, 't', ' ', ' ', ' ', ' ', ' ', '-', '-',
                                '-', ' ', ' ', ' ', ' ', ' ', ' ', 0x005c, 't',
                                ' ', ' ', ' ', '-', '-', '-', ' ', ' ', ' ',
                                ' ', ' ', ' ', ' ', 0x005c, 'n' };
static mxArray * _mxarray90_;

static mxChar _array93_[29] = { 'R', 'o', 'o', 't', ' ', 'Z', 'o', 'n',
                                'e', ' ', ' ', ' ', 0x005c, 't', ' ',
                                ' ', ' ', ' ', ' ', '%', 'g', ' ', ' ',
                                ' ', ' ', ' ', ' ', 0x005c, 't' };
static mxArray * _mxarray92_;

static mxChar _array95_[15] = { ' ', ' ', ' ', ' ', '%', '4', '.', '2',
                                'f', ' ', ' ', ' ', ' ', 0x005c, 't' };
static mxArray * _mxarray94_;

static mxChar _array97_[14] = { ' ', ' ', ' ', ' ', '%', 'g', ' ', ' ',
                                ' ', ' ', ' ', ' ', 0x005c, 'n' };
static mxArray * _mxarray96_;

static mxChar _array99_[29] = { 'S', 't', 'o', 'r', 'a', 'g', 'e', ' ',
                                ' ', ' ', ' ', ' ', 0x005c, 't', ' ',
                                ' ', ' ', ' ', ' ', '%', 'g', ' ', ' ',
                                ' ', ' ', ' ', ' ', 0x005c, 't' };
static mxArray * _mxarray98_;

static mxChar _array101_[16] = { ' ', ' ', ' ', ' ', '%', '4', '.', '2', 'f',
                                 ' ', ' ', ' ', ' ', ' ', 0x005c, 't' };
static mxArray * _mxarray100_;

static mxChar _array103_[44] = { 'N', 'a', 't', 'i', 'v', 'e', ' ', ' ', ' ',
                                 ' ', ' ', ' ', 0x005c, 't', ' ', ' ', ' ', ' ',
                                 ' ', '-', '-', '-', ' ', ' ', ' ', ' ', ' ',
                                 0x005c, 't', ' ', ' ', ' ', ' ', ' ', '%', '4',
                                 '.', '2', 'f', ' ', ' ', ' ', 0x005c, 't' };
static mxArray * _mxarray102_;

static mxChar _array105_[43] = { 0x005c, 'n', 'I', 'n', 'i', 't', '.', ' ', 'C',
                                 'o', 'n', 'd', '.', ':', ' ', 's', 'u', 'c',
                                 't', 'i', 'o', 'n', ' ', 'h', 'e', 'a', 'd',
                                 ' ', '(', 'c', 'm', ')', '=', 0x005c, 't', ' ',
                                 '%', '5', '.', '0', 'f', 0x005c, 'n' };
static mxArray * _mxarray104_;

static mxChar _array107_[39] = { 'M', 'a', 'x', '.', ' ', 'U', 'n', 'd',
                                 'e', 'r', 'd', 'r', 'a', 'i', 'n', ' ',
                                 'f', 'l', 'o', 'w', ' ', '(', 'c', 'm',
                                 '/', 'h', 'r', ')', '=', 0x005c, 't', ' ',
                                 '%', '5', '.', '0', 'f', 0x005c, 'n' };
static mxArray * _mxarray106_;

static mxChar _array109_[31] = { 'U', 'n', 'd', 'e', 'r', 'd', 'r', 'a',
                                 'i', 'n', ' ', 'd', 'i', 'a', 'm', ' ',
                                 '(', 'm', 'm', ')', '=', 0x005c, 't', ' ',
                                 '%', '5', '.', '0', 'f', 0x005c, 'n' };
static mxArray * _mxarray108_;

static mxChar _array111_[25] = { 0x005c, 'n', 'W', 'A', 'T', 'E', 'R', ' ', 'B',
                                 'A', 'L', 'A', 'N', 'C', 'E', ' ', 'T', 'E',
                                 'R', 'M', 'S', 0x005c, 'n', 0x005c, 'n' };
static mxArray * _mxarray110_;

static mxChar _array113_[55] = { ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',
                                 ' ', ' ', ' ', 0x005c, 't', ' ', 'V', 'o',
                                 'l', 'u', 'm', 'e', '(', 'm', '3', ')',
                                 ' ', 0x005c, 't', ' ', 'D', 'e', 'p', 't',
                                 'h', ' ', '(', 'c', 'm', ')', 0x005c, 't',
                                 ' ', '%', '%', 'o', 'f', ' ', 'I', 'n',
                                 'f', 'l', 'o', 'w', ' ', 0x005c, 'n' };
static mxArray * _mxarray112_;

static mxChar _array115_[23] = { 'I', 'n', 'f', 'l', 'o', 'w', ' ', ' ',
                                 ' ', ' ', ' ', 0x005c, 't', ' ', '%', '1',
                                 '0', '.', '2', 'f', ' ', 0x005c, 't' };
static mxArray * _mxarray114_;

static mxChar _array117_[9] = { ' ', '%', '9', '.', '2',
                                'f', ' ', 0x005c, 't' };
static mxArray * _mxarray116_;

static mxChar _array119_[9] = { ' ', '%', '9', '.', '2',
                                'f', ' ', 0x005c, 'n' };
static mxArray * _mxarray118_;

static mxChar _array121_[23] = { 'R', 'u', 'n', 'o', 'f', 'f', ' ', ' ',
                                 ' ', ' ', ' ', 0x005c, 't', ' ', '%', '1',
                                 '0', '.', '2', 'f', ' ', 0x005c, 't' };
static mxArray * _mxarray120_;

static mxChar _array123_[23] = { 'R', 'e', 'c', 'h', 'a', 'r', 'g', 'e',
                                 ' ', ' ', ' ', 0x005c, 't', ' ', '%', '1',
                                 '0', '.', '2', 'f', ' ', 0x005c, 't' };
static mxArray * _mxarray122_;

static mxChar _array125_[23] = { 'E', 'v', 'a', 'p', 'o', 'r', 'a', 't',
                                 'i', 'o', 'n', 0x005c, 't', ' ', '%', '1',
                                 '0', '.', '2', 'f', ' ', 0x005c, 't' };
static mxArray * _mxarray124_;

static mxChar _array127_[23] = { 'U', 'n', 'd', 'e', 'r', 'd', 'r', 'a',
                                 'i', 'n', ' ', 0x005c, 't', ' ', '%', '1',
                                 '0', '.', '2', 'f', ' ', 0x005c, 't' };
static mxArray * _mxarray126_;

static mxChar _array129_[23] = { 'S', 't', 'o', 'r', 'a', 'g', 'e', ' ',
                                 ' ', ' ', ' ', 0x005c, 't', ' ', '%', '1',
                                 '0', '.', '2', 'f', ' ', 0x005c, 't' };
static mxArray * _mxarray128_;

static mxChar _array131_[23] = { 'I', 'n', 'f', 'i', 'l', 't', 'r', 'a',
                                 't', 'e', 'd', 0x005c, 't', ' ', '%', '1',
                                 '0', '.', '2', 'f', ' ', 0x005c, 'n' };
static mxArray * _mxarray130_;

static mxChar _array133_[35] = { 'P', 'r', 'e', 'c', 'i', 'p', '.', ' ', ' ',
                                 ' ', ' ', 0x005c, 't', ' ', ' ', ' ', ' ', ' ',
                                 ' ', ' ', ' ', ' ', ' ', ' ', ' ', 0x005c, 't',
                                 ' ', '%', '9', '.', '2', 'f', 0x005c, 'n' };
static mxArray * _mxarray132_;

static mxChar _array135_[23] = { 'I', 'm', 'p', '.', ' ', 'I', 'n', 'f',
                                 'l', 'o', 'w', 0x005c, 't', ' ', '%', '1',
                                 '0', '.', '2', 'f', ' ', 0x005c, 't' };
static mxArray * _mxarray134_;

static mxChar _array137_[23] = { 'P', 'e', 'r', 'v', '.', 'I', 'n', 'f',
                                 'l', 'o', 'w', 0x005c, 't', ' ', '%', '1',
                                 '0', '.', '2', 'f', ' ', 0x005c, 't' };
static mxArray * _mxarray136_;

static mxChar _array139_[30] = { 'B', 'a', 's', 'i', 'n', ' ', 'S', 't',
                                 'a', 'y', 'o', 'n', '(', 'c', 'm', ')',
                                 0x005c, 't', 0x005c, 't', ' ', '%', '9',
                                 '.', '2', 'f', 0x005c, 'n', 0x005c, 'n' };
static mxArray * _mxarray138_;

static mxChar _array141_[22] = { 'M', 'a', 's', 's', ' ', 'B', 'a', 'l',
                                 'a', 'n', 'c', 'e', ' ', 'C', 'h', 'e',
                                 'c', 'k', 's', ':', 0x005c, 'n' };
static mxArray * _mxarray140_;

static mxChar _array143_[24] = { 'V', 'i', 'n', 's', 'o', 'i', 'l', '=', 'i',
                                 'n', 'p', 'u', 't', '-', 'r', 'o', 'f', 'f',
                                 '-', 'd', 'v', 'h', 0x005c, 'n' };
static mxArray * _mxarray142_;

static mxChar _array145_[20] = { 'V', 'i', 'n', 's', 'o', 'i', 'l',
                                 ' ', '(', 'm', '3', ')', '=', 0x005c,
                                 't', ' ', '%', 'e', 0x005c, 'n' };
static mxArray * _mxarray144_;

static mxChar _array147_[31] = { 'V', 'i', 'n', 'f', ' ', 'd', 'i', 'f',
                                 'f', 'e', 'r', 'e', 'n', 'c', 'e', ' ',
                                 '(', '%', '%', ')', '=', 0x005c, 't', ' ',
                                 '%', '4', '.', '1', 'f', 0x005c, 'n' };
static mxArray * _mxarray146_;

static mxChar _array149_[23] = { 'V', 'r', 'e', 'c', 'h', 'a', 'r', 'g',
                                 '1', '=', 'i', 'n', 'f', '-', 'e', 't',
                                 '-', 'd', 's', 'w', 's', 0x005c, 'n' };
static mxArray * _mxarray148_;

static mxChar _array151_[22] = { 'V', 'r', 'e', 'c', 'h', 'a', 'r', 'g',
                                 '1', ' ', '(', 'm', '3', ')', '=', 0x005c,
                                 't', ' ', '%', 'e', 0x005c, 'n' };
static mxArray * _mxarray150_;

static mxChar _array153_[32] = { 'V', 'r', 'e', 'c', 'h', ' ', 'd', 'i', 'f',
                                 'f', 'e', 'r', 'e', 'n', 'c', 'e', ' ', '(',
                                 '%', '%', ')', '=', 0x005c, 't', ' ', '%',
                                 '4', '.', '1', 'f', 0x005c, 'n' };
static mxArray * _mxarray152_;

static mxChar _array155_[20] = { 'd', 'V', 's', 'o', 'i', 'l', '=',
                                 'i', 'n', 'f', '-', 'r', 'e', 'c',
                                 'h', '-', 'e', 't', 0x005c, 'n' };
static mxArray * _mxarray154_;

static mxChar _array157_[14] = { 'd', 'V', 's', 'o', 'i', 'l', '=', 0x005c,
                                 't', ' ', '%', 'e', 0x005c, 'n' };
static mxArray * _mxarray156_;

static mxChar _array159_[19] = { 'd', 'V', 's', 'f', 'c', 'e', ' ',
                                 '(', 'm', '3', ')', '=', 0x005c,
                                 't', ' ', '%', 'e', 0x005c, 'n' };
static mxArray * _mxarray158_;

static mxChar _array161_[21] = { 'd', 'V', 's', '2', '=', 'i', 'n', 'p',
                                 'u', 't', '-', 'r', 'o', 'f', 'f', '-',
                                 'i', 'n', 'f', 0x005c, 'n' };
static mxArray * _mxarray160_;

static mxChar _array163_[15] = { 'd', 'V', 's', '2', ' ', '(', 'm', '3',
                                 ')', '=', ' ', '%', 'e', 0x005c, 'n' };
static mxArray * _mxarray162_;

static mxChar _array165_[29] = { 0x005c, 'n', 'P', 'L', 'A', 'N', 'T',
                                 ' ', 'S', 'U', 'R', 'V', 'I', 'V', 'A',
                                 'B', 'I', 'L', 'I', 'T', 'Y', ' ', 'T',
                                 'E', 'R', 'M', 'S', 0x005c, 'n' };
static mxArray * _mxarray164_;

static mxChar _array167_[32] = { 't', 'o', 't', 'a', 'l', ' ', 't', 'i', 'm',
                                 'e', ' ', 'p', 'o', 'n', 'd', 'e', 'd', ' ',
                                 '(', 'h', ')', '=', 0x005c, 't', ' ', '%',
                                 '5', '.', '1', 'f', 0x005c, 'n' };
static mxArray * _mxarray166_;

static mxChar _array169_[30] = { 'm', 'a', 'x', ' ', 't', 'i', 'm', 'e',
                                 ' ', 'p', 'o', 'n', 'd', 'e', 'd', ' ',
                                 '(', 'h', ')', '=', 0x005c, 't', ' ',
                                 '%', '5', '.', '1', 'f', 0x005c, 'n' };
static mxArray * _mxarray168_;

static mxChar _array171_[38] = { 't', 'o', 't', 'a', 'l', ' ', 't', 'i',
                                 'm', 'e', ' ', 'R', 'Z', ' ', 's', 'a',
                                 't', 'u', 'r', 'a', 't', 'e', 'd', ' ',
                                 '(', 'h', ')', '=', 0x005c, 't', ' ',
                                 '%', '5', '.', '1', 'f', 0x005c, 'n' };
static mxArray * _mxarray170_;

static mxChar _array173_[45] = { 't', 'o', 't', 'a', 'l', ' ', 't', 'i',
                                 'm', 'e', ' ', 'R', 'Z', ' ', 'a', 't',
                                 ' ', 'w', 'i', 'l', 't', 'i', 'n', 'g',
                                 ' ', 'p', 'o', 'i', 'n', 't', ' ', '(',
                                 'h', ')', '=', 0x005c, 't', ' ', '%',
                                 '5', '.', '1', 'f', 0x005c, 'n' };
static mxArray * _mxarray172_;

static mxChar _array175_[38] = { 't', 'o', 't', 'a', 'l', ' ', 't', 'i',
                                 'm', 'e', ' ', 'T', 'Z', ' ', 's', 'a',
                                 't', 'u', 'r', 'a', 't', 'e', 'd', ' ',
                                 '(', 'h', ')', '=', 0x005c, 't', ' ',
                                 '%', '5', '.', '1', 'f', 0x005c, 'n' };
static mxArray * _mxarray174_;

static mxChar _array177_[36] = { 'm', 'a', 'x', ' ', 't', 'i', 'm', 'e',
                                 ' ', 'R', 'Z', ' ', 's', 'a', 't', 'u',
                                 'r', 'a', 't', 'e', 'd', ' ', '(', 'h',
                                 ')', '=', 0x005c, 't', ' ', '%', '5',
                                 '.', '1', 'f', 0x005c, 'n' };
static mxArray * _mxarray176_;

static mxChar _array179_[43] = { 'm', 'a', 'x', ' ', 't', 'i', 'm', 'e', ' ',
                                 'R', 'Z', ' ', 'a', 't', ' ', 'w', 'i', 'l',
                                 't', 'i', 'n', 'g', ' ', 'p', 'o', 'i', 'n',
                                 't', ' ', '(', 'h', ')', '=', 0x005c, 't',
                                 ' ', '%', '5', '.', '1', 'f', 0x005c, 'n' };
static mxArray * _mxarray178_;

static mxChar _array181_[36] = { 'm', 'a', 'x', ' ', 't', 'i', 'm', 'e',
                                 ' ', 'T', 'Z', ' ', 's', 'a', 't', 'u',
                                 'r', 'a', 't', 'e', 'd', ' ', '(', 'h',
                                 ')', '=', 0x005c, 't', ' ', '%', '5',
                                 '.', '1', 'f', 0x005c, 'n' };
static mxArray * _mxarray180_;

static mxChar _array183_[46] = { 'n', 'e', 'a', 'r', ' ', 's', 'a', 't',
                                 'u', 'r', 'a', 't', 'e', 'd', ' ', 't',
                                 'i', 'm', 'e', 's', ' ', 'f', 'o', 'r',
                                 ' ', 's', 'i', 'm', 'u', 'l', 'a', 't',
                                 'i', 'o', 'n', '=', 0x005c, 't', ' ',
                                 '%', '6', '.', '0', 'f', 0x005c, 'n' };
static mxArray * _mxarray182_;

static mxChar _array185_[15] = { '%', '6', '.', '1', 'f', 0x005c, 't', ' ',
                                 '%', '6', '.', '1', 'f', 0x005c, 'n' };
static mxArray * _mxarray184_;

static mxChar _array187_[45] = { 'w', 'i', 'l', 't', 'i', 'n', 'g', ' ',
                                 'p', 'o', 'i', 'n', 't', ' ', 't', 'i',
                                 'm', 'e', 's', ' ', 'f', 'o', 'r', ' ',
                                 's', 'i', 'm', 'u', 'l', 'a', 't', 'i',
                                 'o', 'n', '=', 0x005c, 't', ' ', '%',
                                 '6', '.', '0', 'f', 0x005c, 'n' };
static mxArray * _mxarray186_;

static mxChar _array189_[38] = { 'p', 'o', 'n', 'd', 'e', 'd', ' ', 't',
                                 'i', 'm', 'e', 's', ' ', 'f', 'o', 'r',
                                 ' ', 's', 'i', 'm', 'u', 'l', 'a', 't',
                                 'i', 'o', 'n', '=', 0x005c, 't', ' ',
                                 '%', '6', '.', '0', 'f', 0x005c, 'n' };
static mxArray * _mxarray188_;

static mxChar _array191_[40] = { 'o', 'v', 'e', 'r', 'f', 'l', 'o', 'w', ' ',
                                 't', 'i', 'm', 'e', 's', ' ', 'f', 'o', 'r',
                                 ' ', 's', 'i', 'm', 'u', 'l', 'a', 't', 'i',
                                 'o', 'n', '=', 0x005c, 't', ' ', '%', '6',
                                 '.', '0', 'f', 0x005c, 'n' };
static mxArray * _mxarray190_;

static mxChar _array193_[7] = { '_', 'W', 'B', '.', 't', 'x', 't' };
static mxArray * _mxarray192_;

static mxChar _array195_[127] = { 'T', 'i', 'm', 'e', '(', 'h', 'r', ')',
                                  0x005c, 't', ' ', 'R', 'u', 'n', 'o', 'n',
                                  '(', 'c', 'm', ')', 0x005c, 't', ' ', 'P',
                                  'o', 'n', 'd', 'i', 'n', 'g', '(', 'c',
                                  'm', ')', 0x005c, 't', ' ', 'I', 'n', 'f',
                                  'i', 'l', '(', 'c', 'm', ')', 0x005c, 't',
                                  ' ', 'R', 'u', 'n', 'o', 'f', 'f', '(',
                                  'c', 'm', ')', 0x005c, 't', ' ', 'D', 'r',
                                  'a', 'i', 'n', '(', 'c', 'm', ')', 0x005c,
                                  't', ' ', 'R', 'e', 'c', 'h', 'a', 'r',
                                  'g', 'e', '(', 'c', 'm', ')', 0x005c, 't',
                                  ' ', 'E', 'T', '(', 'c', 'm', ')', 0x005c,
                                  't', ' ', 'T', 'h', 'e', 't', 'a', 'R',
                                  'Z', 0x005c, 't', ' ', 'T', 'h', 'e', 't',
                                  'a', 'S', 'Z', 0x005c, 't', ' ', 'T', 'h',
                                  'e', 't', 'a', 'C', 'Z', 0x005c, 'n' };
static mxArray * _mxarray194_;

static mxChar _array197_[86] = { '%', '6', '.', '3', 'f', 0x005c, 't', ' ',
                                 '%', '6', '.', '3', 'f', 0x005c, 't', ' ',
                                 '%', '6', '.', '3', 'f', 0x005c, 't', ' ',
                                 '%', '6', '.', '3', 'f', 0x005c, 't', ' ',
                                 '%', '6', '.', '3', 'f', 0x005c, 't', ' ',
                                 '%', '6', '.', '3', 'f', 0x005c, 't', ' ',
                                 '%', '6', '.', '3', 'f', 0x005c, 't', ' ',
                                 '%', '6', '.', '3', 'f', 0x005c, 't', ' ',
                                 '%', '6', '.', '3', 'f', 0x005c, 't', ' ',
                                 '%', '6', '.', '3', 'f', 0x005c, 't', ' ',
                                 '%', '.', '3', 'f', 0x005c, 'n' };
static mxArray * _mxarray196_;

void InitializeModule_recarga(void) {
    _mxarray0_ = mclInitializeString(4, _array1_);
    _mxarray2_ = mclInitializeDouble(9.8);
    _mxarray3_ = mclInitializeDouble(.041666666666666664);
    _mxarray4_ = mclInitializeDouble(24.0);
    _mxarray5_ = mclInitializeDouble(.75);
    _mxarray6_ = mclInitializeString(12, _array7_);
    _mxarray8_ = mclInitializeString(2, _array9_);
    _mxarray10_ = mclInitializeDoubleVector(1, 4, _array11_);
    _mxarray12_ = mclInitializeDouble(1.0);
    _mxarray13_ = mclInitializeDouble(2.0);
    _mxarray14_ = mclInitializeDoubleVector(1, 3, _array15_);
    _mxarray16_ = mclInitializeDoubleVector(1, 3, _array17_);
    _mxarray18_ = mclInitializeDouble(-1.0);
    _mxarray19_ = mclInitializeDouble(3.1);
    _mxarray20_ = mclInitializeDouble(.1);
    _mxarray21_ = mclInitializeDouble(1.21);
    _mxarray22_ = mclInitializeDouble(.95);
    _mxarray23_ = mclInitializeDouble(60.0);
    _mxarray24_ = mclInitializeDouble(15.0);
    _mxarray25_ = mclInitializeDouble(10.0);
    _mxarray26_ = mclInitializeDoubleVector(1, 3, _array27_);
    _mxarray28_ = mclInitializeDoubleVector(1, 3, _array29_);
    _mxarray30_ = mclInitializeDouble(0.0);
    _mxarray31_ = mclInitializeDouble(5.0);
    _mxarray32_ = mclInitializeDoubleVector(1, 3, _array33_);
    _mxarray34_ = mclInitializeDoubleVector(1, 3, _array35_);
    _mxarray36_ = mclInitializeDouble(3.0);
    _mxarray37_ = mclInitializeDouble(-100.0);
    _mxarray38_ = mclInitializeDoubleVector(0, 0, (double *)NULL);
    _mxarray39_ = mclInitializeDouble(2.54);
    _mxarray40_ = mclInitializeDouble(1000.0);
    _mxarray41_ = mclInitializeDouble(100.0);
    _mxarray42_ = mclInitializeDouble(3600.0);
    _mxarray43_ = mclInitializeDouble(.62);
    _mxarray44_ = mclInitializeDouble(9.81);
    _mxarray45_ = mclInitializeDouble(3.14159);
    _mxarray46_ = mclInitializeDouble(6.0);
    _mxarray47_ = mclInitializeDoubleVector(1, 11, _array48_);
    _mxarray49_ = mclInitializeDoubleVector(1, 2, _array50_);
    _mxarray51_ = mclInitializeString(6, _array52_);
    _mxarray53_ = mclInitializeString(7, _array54_);
    _mxarray55_ = mclInitializeDouble(4.0);
    _mxarray56_ = mclInitializeDouble(.5);
    _mxarray57_ = mclInitializeDouble(1.3);
    _mxarray58_ = mclInitializeString(4, _array59_);
    _mxarray60_ = mclInitializeString(1, _array61_);
    _mxarray62_ = mclInitializeString(2, _array63_);
    _mxarray64_ = mclInitializeString(13, _array65_);
    _mxarray66_ = mclInitializeString(12, _array67_);
    _mxarray68_ = mclInitializeString(19, _array69_);
    _mxarray70_ = mclInitializeString(32, _array71_);
    _mxarray72_ = mclInitializeString(35, _array73_);
    _mxarray74_ = mclInitializeString(13, _array75_);
    _mxarray76_ = mclInitializeString(20, _array77_);
    _mxarray78_ = mclInitializeString(19, _array79_);
    _mxarray80_ = mclInitializeString(26, _array81_);
    _mxarray82_ = mclInitializeString(23, _array83_);
    _mxarray84_ = mclInitializeString(19, _array85_);
    _mxarray86_ = mclInitializeString(22, _array87_);
    _mxarray88_ = mclInitializeString(59, _array89_);
    _mxarray90_ = mclInitializeString(60, _array91_);
    _mxarray92_ = mclInitializeString(29, _array93_);
    _mxarray94_ = mclInitializeString(15, _array95_);
    _mxarray96_ = mclInitializeString(14, _array97_);
    _mxarray98_ = mclInitializeString(29, _array99_);
    _mxarray100_ = mclInitializeString(16, _array101_);
    _mxarray102_ = mclInitializeString(44, _array103_);
    _mxarray104_ = mclInitializeString(43, _array105_);
    _mxarray106_ = mclInitializeString(39, _array107_);
    _mxarray108_ = mclInitializeString(31, _array109_);
    _mxarray110_ = mclInitializeString(25, _array111_);
    _mxarray112_ = mclInitializeString(55, _array113_);
    _mxarray114_ = mclInitializeString(23, _array115_);
    _mxarray116_ = mclInitializeString(9, _array117_);
    _mxarray118_ = mclInitializeString(9, _array119_);
    _mxarray120_ = mclInitializeString(23, _array121_);
    _mxarray122_ = mclInitializeString(23, _array123_);
    _mxarray124_ = mclInitializeString(23, _array125_);
    _mxarray126_ = mclInitializeString(23, _array127_);
    _mxarray128_ = mclInitializeString(23, _array129_);
    _mxarray130_ = mclInitializeString(23, _array131_);
    _mxarray132_ = mclInitializeString(35, _array133_);
    _mxarray134_ = mclInitializeString(23, _array135_);
    _mxarray136_ = mclInitializeString(23, _array137_);
    _mxarray138_ = mclInitializeString(30, _array139_);
    _mxarray140_ = mclInitializeString(22, _array141_);
    _mxarray142_ = mclInitializeString(24, _array143_);
    _mxarray144_ = mclInitializeString(20, _array145_);
    _mxarray146_ = mclInitializeString(31, _array147_);
    _mxarray148_ = mclInitializeString(23, _array149_);
    _mxarray150_ = mclInitializeString(22, _array151_);
    _mxarray152_ = mclInitializeString(32, _array153_);
    _mxarray154_ = mclInitializeString(20, _array155_);
    _mxarray156_ = mclInitializeString(14, _array157_);
    _mxarray158_ = mclInitializeString(19, _array159_);
    _mxarray160_ = mclInitializeString(21, _array161_);
    _mxarray162_ = mclInitializeString(15, _array163_);
    _mxarray164_ = mclInitializeString(29, _array165_);
    _mxarray166_ = mclInitializeString(32, _array167_);
    _mxarray168_ = mclInitializeString(30, _array169_);
    _mxarray170_ = mclInitializeString(38, _array171_);
    _mxarray172_ = mclInitializeString(45, _array173_);
    _mxarray174_ = mclInitializeString(38, _array175_);
    _mxarray176_ = mclInitializeString(36, _array177_);
    _mxarray178_ = mclInitializeString(43, _array179_);
    _mxarray180_ = mclInitializeString(36, _array181_);
    _mxarray182_ = mclInitializeString(46, _array183_);
    _mxarray184_ = mclInitializeString(15, _array185_);
    _mxarray186_ = mclInitializeString(45, _array187_);
    _mxarray188_ = mclInitializeString(38, _array189_);
    _mxarray190_ = mclInitializeString(40, _array191_);
    _mxarray192_ = mclInitializeString(7, _array193_);
    _mxarray194_ = mclInitializeString(127, _array195_);
    _mxarray196_ = mclInitializeString(86, _array197_);
}

void TerminateModule_recarga(void) {
    mxDestroyArray(_mxarray196_);
    mxDestroyArray(_mxarray194_);
    mxDestroyArray(_mxarray192_);
    mxDestroyArray(_mxarray190_);
    mxDestroyArray(_mxarray188_);
    mxDestroyArray(_mxarray186_);
    mxDestroyArray(_mxarray184_);
    mxDestroyArray(_mxarray182_);
    mxDestroyArray(_mxarray180_);
    mxDestroyArray(_mxarray178_);
    mxDestroyArray(_mxarray176_);
    mxDestroyArray(_mxarray174_);
    mxDestroyArray(_mxarray172_);
    mxDestroyArray(_mxarray170_);
    mxDestroyArray(_mxarray168_);
    mxDestroyArray(_mxarray166_);
    mxDestroyArray(_mxarray164_);
    mxDestroyArray(_mxarray162_);
    mxDestroyArray(_mxarray160_);
    mxDestroyArray(_mxarray158_);
    mxDestroyArray(_mxarray156_);
    mxDestroyArray(_mxarray154_);
    mxDestroyArray(_mxarray152_);
    mxDestroyArray(_mxarray150_);
    mxDestroyArray(_mxarray148_);
    mxDestroyArray(_mxarray146_);
    mxDestroyArray(_mxarray144_);
    mxDestroyArray(_mxarray142_);
    mxDestroyArray(_mxarray140_);
    mxDestroyArray(_mxarray138_);
    mxDestroyArray(_mxarray136_);
    mxDestroyArray(_mxarray134_);
    mxDestroyArray(_mxarray132_);
    mxDestroyArray(_mxarray130_);
    mxDestroyArray(_mxarray128_);
    mxDestroyArray(_mxarray126_);
    mxDestroyArray(_mxarray124_);
    mxDestroyArray(_mxarray122_);
    mxDestroyArray(_mxarray120_);
    mxDestroyArray(_mxarray118_);
    mxDestroyArray(_mxarray116_);
    mxDestroyArray(_mxarray114_);
    mxDestroyArray(_mxarray112_);
    mxDestroyArray(_mxarray110_);
    mxDestroyArray(_mxarray108_);
    mxDestroyArray(_mxarray106_);
    mxDestroyArray(_mxarray104_);
    mxDestroyArray(_mxarray102_);
    mxDestroyArray(_mxarray100_);
    mxDestroyArray(_mxarray98_);
    mxDestroyArray(_mxarray96_);
    mxDestroyArray(_mxarray94_);
    mxDestroyArray(_mxarray92_);
    mxDestroyArray(_mxarray90_);
    mxDestroyArray(_mxarray88_);
    mxDestroyArray(_mxarray86_);
    mxDestroyArray(_mxarray84_);
    mxDestroyArray(_mxarray82_);
    mxDestroyArray(_mxarray80_);
    mxDestroyArray(_mxarray78_);
    mxDestroyArray(_mxarray76_);
    mxDestroyArray(_mxarray74_);
    mxDestroyArray(_mxarray72_);
    mxDestroyArray(_mxarray70_);
    mxDestroyArray(_mxarray68_);
    mxDestroyArray(_mxarray66_);
    mxDestroyArray(_mxarray64_);
    mxDestroyArray(_mxarray62_);
    mxDestroyArray(_mxarray60_);
    mxDestroyArray(_mxarray58_);
    mxDestroyArray(_mxarray57_);
    mxDestroyArray(_mxarray56_);
    mxDestroyArray(_mxarray55_);
    mxDestroyArray(_mxarray53_);
    mxDestroyArray(_mxarray51_);
    mxDestroyArray(_mxarray49_);
    mxDestroyArray(_mxarray47_);
    mxDestroyArray(_mxarray46_);
    mxDestroyArray(_mxarray45_);
    mxDestroyArray(_mxarray44_);
    mxDestroyArray(_mxarray43_);
    mxDestroyArray(_mxarray42_);
    mxDestroyArray(_mxarray41_);
    mxDestroyArray(_mxarray40_);
    mxDestroyArray(_mxarray39_);
    mxDestroyArray(_mxarray38_);
    mxDestroyArray(_mxarray37_);
    mxDestroyArray(_mxarray36_);
    mxDestroyArray(_mxarray34_);
    mxDestroyArray(_mxarray32_);
    mxDestroyArray(_mxarray31_);
    mxDestroyArray(_mxarray30_);
    mxDestroyArray(_mxarray28_);
    mxDestroyArray(_mxarray26_);
    mxDestroyArray(_mxarray25_);
    mxDestroyArray(_mxarray24_);
    mxDestroyArray(_mxarray23_);
    mxDestroyArray(_mxarray22_);
    mxDestroyArray(_mxarray21_);
    mxDestroyArray(_mxarray20_);
    mxDestroyArray(_mxarray19_);
    mxDestroyArray(_mxarray18_);
    mxDestroyArray(_mxarray16_);
    mxDestroyArray(_mxarray14_);
    mxDestroyArray(_mxarray13_);
    mxDestroyArray(_mxarray12_);
    mxDestroyArray(_mxarray10_);
    mxDestroyArray(_mxarray8_);
    mxDestroyArray(_mxarray6_);
    mxDestroyArray(_mxarray5_);
    mxDestroyArray(_mxarray4_);
    mxDestroyArray(_mxarray3_);
    mxDestroyArray(_mxarray2_);
    mxDestroyArray(_mxarray0_);
}

static mxArray * Mrecarga(mxArray * * Dimp,
                          mxArray * * Dperv,
                          mxArray * * Drunon,
                          mxArray * * Drunoff,
                          mxArray * * Devap,
                          mxArray * * Drecharge,
                          mxArray * * Dudrain,
                          mxArray * * Dstore,
                          mxArray * * tctpmax,
                          mxArray * * tstponded,
                          mxArray * * inot,
                          int nargout_,
                          mxArray * nprecip,
                          mxArray * Sd,
                          mxArray * avet,
                          mxArray * rundays,
                          mxArray * outpre,
                          mxArray * Summary,
                          mxArray * Record,
                          mxArray * dmin,
                          mxArray * Srz,
                          mxArray * St,
                          mxArray * Rsoilnumber,
                          mxArray * Ssoilnumber,
                          mxArray * Nsoilnumber,
                          mxArray * Krz,
                          mxArray * Kst,
                          mxArray * Kcz,
                          mxArray * Urate,
                          mxArray * CN,
                          mxArray * EIA,
                          mxArray * Arg,
                          mxArray * TribA);

_mexLocalFunctionTable _local_function_table_recarga
  = { 0, (mexFunctionTableEntry *)NULL };

/*
 * The function "mlfRecarga" contains the normal interface for the "recarga"
 * M-function from file "u:\final_recarga6\recarga\recarga.m" (lines 1-1229).
 * This function processes any input arguments and passes them to the
 * implementation version of the function, appearing above.
 */
mxArray * mlfRecarga(mxArray * * Dimp,
                     mxArray * * Dperv,
                     mxArray * * Drunon,
                     mxArray * * Drunoff,
                     mxArray * * Devap,
                     mxArray * * Drecharge,
                     mxArray * * Dudrain,
                     mxArray * * Dstore,
                     mxArray * * tctpmax,
                     mxArray * * tstponded,
                     mxArray * * inot,
                     mxArray * nprecip,
                     mxArray * Sd,
                     mxArray * avet,
                     mxArray * rundays,
                     mxArray * outpre,
                     mxArray * Summary,
                     mxArray * Record,
                     mxArray * dmin,
                     mxArray * Srz,
                     mxArray * St,
                     mxArray * Rsoilnumber,
                     mxArray * Ssoilnumber,
                     mxArray * Nsoilnumber,
                     mxArray * Krz,
                     mxArray * Kst,
                     mxArray * Kcz,
                     mxArray * Urate,
                     mxArray * CN,
                     mxArray * EIA,
                     mxArray * Arg,
                     mxArray * TribA) {
    int nargout = 1;
    mxArray * Dprecip = NULL;
    mxArray * Dimp__ = NULL;
    mxArray * Dperv__ = NULL;
    mxArray * Drunon__ = NULL;
    mxArray * Drunoff__ = NULL;
    mxArray * Devap__ = NULL;
    mxArray * Drecharge__ = NULL;
    mxArray * Dudrain__ = NULL;
    mxArray * Dstore__ = NULL;
    mxArray * tctpmax__ = NULL;
    mxArray * tstponded__ = NULL;
    mxArray * inot__ = NULL;
    mlfEnterNewContext(
      11,
      21,
      Dimp,
      Dperv,
      Drunon,
      Drunoff,
      Devap,
      Drecharge,
      Dudrain,
      Dstore,
      tctpmax,
      tstponded,
      inot,
      nprecip,
      Sd,
      avet,
      rundays,
      outpre,
      Summary,
      Record,
      dmin,
      Srz,
      St,
      Rsoilnumber,
      Ssoilnumber,
      Nsoilnumber,
      Krz,
      Kst,
      Kcz,
      Urate,
      CN,
      EIA,
      Arg,
      TribA);
    if (Dimp != NULL) {
        ++nargout;
    }
    if (Dperv != NULL) {
        ++nargout;
    }
    if (Drunon != NULL) {
        ++nargout;
    }
    if (Drunoff != NULL) {
        ++nargout;
    }
    if (Devap != NULL) {
        ++nargout;
    }
    if (Drecharge != NULL) {
        ++nargout;
    }
    if (Dudrain != NULL) {
        ++nargout;
    }
    if (Dstore != NULL) {
        ++nargout;
    }
    if (tctpmax != NULL) {
        ++nargout;
    }
    if (tstponded != NULL) {
        ++nargout;
    }
    if (inot != NULL) {
        ++nargout;
    }
    Dprecip
      = Mrecarga(
          &Dimp__,
          &Dperv__,
          &Drunon__,
          &Drunoff__,
          &Devap__,
          &Drecharge__,
          &Dudrain__,
          &Dstore__,
          &tctpmax__,
          &tstponded__,
          &inot__,
          nargout,
          nprecip,
          Sd,
          avet,
          rundays,
          outpre,
          Summary,
          Record,
          dmin,
          Srz,
          St,
          Rsoilnumber,
          Ssoilnumber,
          Nsoilnumber,
          Krz,
          Kst,
          Kcz,
          Urate,
          CN,
          EIA,
          Arg,
          TribA);
    mlfRestorePreviousContext(
      11,
      21,
      Dimp,
      Dperv,
      Drunon,
      Drunoff,
      Devap,
      Drecharge,
      Dudrain,
      Dstore,
      tctpmax,
      tstponded,
      inot,
      nprecip,
      Sd,
      avet,
      rundays,
      outpre,
      Summary,
      Record,
      dmin,
      Srz,
      St,
      Rsoilnumber,
      Ssoilnumber,
      Nsoilnumber,
      Krz,
      Kst,
      Kcz,
      Urate,
      CN,
      EIA,
      Arg,
      TribA);
    if (Dimp != NULL) {
        mclCopyOutputArg(Dimp, Dimp__);
    } else {
        mxDestroyArray(Dimp__);
    }
    if (Dperv != NULL) {
        mclCopyOutputArg(Dperv, Dperv__);
    } else {
        mxDestroyArray(Dperv__);
    }
    if (Drunon != NULL) {
        mclCopyOutputArg(Drunon, Drunon__);
    } else {
        mxDestroyArray(Drunon__);
    }
    if (Drunoff != NULL) {
        mclCopyOutputArg(Drunoff, Drunoff__);
    } else {
        mxDestroyArray(Drunoff__);
    }
    if (Devap != NULL) {
        mclCopyOutputArg(Devap, Devap__);
    } else {
        mxDestroyArray(Devap__);
    }
    if (Drecharge != NULL) {
        mclCopyOutputArg(Drecharge, Drecharge__);
    } else {
        mxDestroyArray(Drecharge__);
    }
    if (Dudrain != NULL) {
        mclCopyOutputArg(Dudrain, Dudrain__);
    } else {
        mxDestroyArray(Dudrain__);
    }
    if (Dstore != NULL) {
        mclCopyOutputArg(Dstore, Dstore__);
    } else {
        mxDestroyArray(Dstore__);
    }
    if (tctpmax != NULL) {
        mclCopyOutputArg(tctpmax, tctpmax__);
    } else {
        mxDestroyArray(tctpmax__);
    }
    if (tstponded != NULL) {
        mclCopyOutputArg(tstponded, tstponded__);
    } else {
        mxDestroyArray(tstponded__);
    }
    if (inot != NULL) {
        mclCopyOutputArg(inot, inot__);
    } else {
        mxDestroyArray(inot__);
    }
    return mlfReturnValue(Dprecip);
}

/*
 * The function "mlxRecarga" contains the feval interface for the "recarga"
 * M-function from file "u:\final_recarga6\recarga\recarga.m" (lines 1-1229).
 * The feval function calls the implementation version of recarga through this
 * function. This function processes any input arguments and passes them to the
 * implementation version of the function, appearing above.
 */
void mlxRecarga(int nlhs, mxArray * plhs[], int nrhs, mxArray * prhs[]) {
    mxArray * mprhs[21];
    mxArray * mplhs[12];
    int i;
    if (nlhs > 12) {
        mlfError(
          mxCreateString(
            "Run-time Error: File: recarga Line: 1 Column: "
            "1 The function \"recarga\" was called with mor"
            "e than the declared number of outputs (12)."),
          NULL);
    }
    if (nrhs > 21) {
        mlfError(
          mxCreateString(
            "Run-time Error: File: recarga Line: 1 Column: "
            "1 The function \"recarga\" was called with mor"
            "e than the declared number of inputs (21)."),
          NULL);
    }
    for (i = 0; i < 12; ++i) {
        mplhs[i] = NULL;
    }
    for (i = 0; i < 21 && i < nrhs; ++i) {
        mprhs[i] = prhs[i];
    }
    for (; i < 21; ++i) {
        mprhs[i] = NULL;
    }
    mlfEnterNewContext(
      0,
      21,
      mprhs[0],
      mprhs[1],
      mprhs[2],
      mprhs[3],
      mprhs[4],
      mprhs[5],
      mprhs[6],
      mprhs[7],
      mprhs[8],
      mprhs[9],
      mprhs[10],
      mprhs[11],
      mprhs[12],
      mprhs[13],
      mprhs[14],
      mprhs[15],
      mprhs[16],
      mprhs[17],
      mprhs[18],
      mprhs[19],
      mprhs[20]);
    mplhs[0]
      = Mrecarga(
          &mplhs[1],
          &mplhs[2],
          &mplhs[3],
          &mplhs[4],
          &mplhs[5],
          &mplhs[6],
          &mplhs[7],
          &mplhs[8],
          &mplhs[9],
          &mplhs[10],
          &mplhs[11],
          nlhs,
          mprhs[0],
          mprhs[1],
          mprhs[2],
          mprhs[3],
          mprhs[4],
          mprhs[5],
          mprhs[6],
          mprhs[7],
          mprhs[8],
          mprhs[9],
          mprhs[10],
          mprhs[11],
          mprhs[12],
          mprhs[13],
          mprhs[14],
          mprhs[15],
          mprhs[16],
          mprhs[17],
          mprhs[18],
          mprhs[19],
          mprhs[20]);
    mlfRestorePreviousContext(
      0,
      21,
      mprhs[0],
      mprhs[1],
      mprhs[2],
      mprhs[3],
      mprhs[4],
      mprhs[5],
      mprhs[6],
      mprhs[7],
      mprhs[8],
      mprhs[9],
      mprhs[10],
      mprhs[11],
      mprhs[12],
      mprhs[13],
      mprhs[14],
      mprhs[15],
      mprhs[16],
      mprhs[17],
      mprhs[18],
      mprhs[19],
      mprhs[20]);
    plhs[0] = mplhs[0];
    for (i = 1; i < 12 && i < nlhs; ++i) {
        plhs[i] = mplhs[i];
    }
    for (; i < 12; ++i) {
        mxDestroyArray(mplhs[i]);
    }
}

/*
 * The function "Mrecarga" is the implementation version of the "recarga"
 * M-function from file "u:\final_recarga6\recarga\recarga.m" (lines 1-1229).
 * It contains the actual compiled code for that M-function. It is a static
 * function and must only be called from one of the interface functions,
 * appearing below.
 */
/*
 * function [Dprecip,Dimp,Dperv,Drunon,Drunoff,Devap,Drecharge,Dudrain,Dstore,tctpmax,tstponded,inot] = recarga(nprecip,Sd,avet,rundays,outpre,Summary,Record,dmin,Srz,St,Rsoilnumber,Ssoilnumber,Nsoilnumber,Krz,Kst,Kcz,Urate,CN,EIA,Arg,TribA);
 */
static mxArray * Mrecarga(mxArray * * Dimp,
                          mxArray * * Dperv,
                          mxArray * * Drunon,
                          mxArray * * Drunoff,
                          mxArray * * Devap,
                          mxArray * * Drecharge,
                          mxArray * * Dudrain,
                          mxArray * * Dstore,
                          mxArray * * tctpmax,
                          mxArray * * tstponded,
                          mxArray * * inot,
                          int nargout_,
                          mxArray * nprecip,
                          mxArray * Sd,
                          mxArray * avet,
                          mxArray * rundays,
                          mxArray * outpre,
                          mxArray * Summary,
                          mxArray * Record,
                          mxArray * dmin,
                          mxArray * Srz,
                          mxArray * St,
                          mxArray * Rsoilnumber,
                          mxArray * Ssoilnumber,
                          mxArray * Nsoilnumber,
                          mxArray * Krz,
                          mxArray * Kst,
                          mxArray * Kcz,
                          mxArray * Urate,
                          mxArray * CN,
                          mxArray * EIA,
                          mxArray * Arg,
                          mxArray * TribA) {
    mexLocalFunctionTable save_local_function_table_
      = mclSetCurrentLocalFunctionTable(&_local_function_table_recarga);
    mxArray * Dprecip = NULL;
    mxArray * WBfile = NULL;
    mxArray * fid = NULL;
    mxArray * outfile = NULL;
    mxArray * ttcpu = NULL;
    mxArray * dVs = NULL;
    mxArray * unitud = NULL;
    mxArray * unitroff = NULL;
    mxArray * unitron = NULL;
    mxArray * unitet = NULL;
    mxArray * dVsoil = NULL;
    mxArray * unitgwr = NULL;
    mxArray * gwravgdepth = NULL;
    mxArray * gwrdiff = NULL;
    mxArray * Vr = NULL;
    mxArray * infdiff = NULL;
    mxArray * Vinsoil = NULL;
    mxArray * dVsfce = NULL;
    mxArray * unitSM = NULL;
    mxArray * DSSW = NULL;
    mxArray * Pcompl = NULL;
    mxArray * DSSW3 = NULL;
    mxArray * DSSW2 = NULL;
    mxArray * DSSW1 = NULL;
    mxArray * et = NULL;
    mxArray * Pevap = NULL;
    mxArray * aw = NULL;
    mxArray * hpold = NULL;
    mxArray * xswater = NULL;
    mxArray * leftover = NULL;
    mxArray * deficit = NULL;
    mxArray * theta_CZold = NULL;
    mxArray * theta_TZold = NULL;
    mxArray * theta_RZold = NULL;
    mxArray * uflow = NULL;
    mxArray * uflowmax = NULL;
    mxArray * uarea = NULL;
    mxArray * hu = NULL;
    mxArray * F1 = NULL;
    mxArray * A = NULL;
    mxArray * Fp = NULL;
    mxArray * FGA = NULL;
    mxArray * test = NULL;
    mxArray * T = NULL;
    mxArray * thetainiCZ = NULL;
    mxArray * thetainiTZ = NULL;
    mxArray * thetainiRZ = NULL;
    mxArray * te = NULL;
    mxArray * test3 = NULL;
    mxArray * test2 = NULL;
    mxArray * tcpu = NULL;
    mxArray * n = NULL;
    mxArray * ttstorm = NULL;
    mxArray * indstorm = NULL;
    mxArray * inp = NULL;
    mxArray * runon = NULL;
    mxArray * Dbasin = NULL;
    mxArray * Vperv = NULL;
    mxArray * runon2 = NULL;
    mxArray * prain = NULL;
    mxArray * pdatafile = NULL;
    mxArray * Gap_P = NULL;
    mxArray * rain = NULL;
    mxArray * rgrain = NULL;
    mxArray * Vimp = NULL;
    mxArray * runon1 = NULL;
    mxArray * evap = NULL;
    mxArray * irain = NULL;
    mxArray * datafile = NULL;
    mxArray * intfile = NULL;
    mxArray * intensity = NULL;
    mxArray * DARain = NULL;
    mxArray * CummRain = NULL;
    mxArray * tk = NULL;
    mxArray * len = NULL;
    mxArray * pdfile = NULL;
    mxArray * es = NULL;
    mxArray * qs = NULL;
    mxArray * wtimes = NULL;
    mxArray * inwt = NULL;
    mxArray * prevof = NULL;
    mxArray * oftimes = NULL;
    mxArray * ptimes = NULL;
    mxArray * inpt = NULL;
    mxArray * nstimes = NULL;
    mxArray * inst = NULL;
    mxArray * tctrzwiltmax = NULL;
    mxArray * tctrzwilt = NULL;
    mxArray * tstrzwilt = NULL;
    mxArray * tcttznsat = NULL;
    mxArray * tctrznsat = NULL;
    mxArray * tcttznsatmax = NULL;
    mxArray * tctrznsatmax = NULL;
    mxArray * tsttznsat = NULL;
    mxArray * tstrznsat = NULL;
    mxArray * tctoverflow = NULL;
    mxArray * tctponded = NULL;
    mxArray * V_pond = NULL;
    mxArray * V_et = NULL;
    mxArray * V_ud = NULL;
    mxArray * V_rech = NULL;
    mxArray * V_roff = NULL;
    mxArray * V_inf = NULL;
    mxArray * V_inp = NULL;
    mxArray * perc3lim = NULL;
    mxArray * perc2lim = NULL;
    mxArray * perc1lim = NULL;
    mxArray * s = NULL;
    mxArray * x = NULL;
    mxArray * Recharge = NULL;
    mxArray * EvapoT = NULL;
    mxArray * Runoff = NULL;
    mxArray * infilt = NULL;
    mxArray * udrain = NULL;
    mxArray * Inflow = NULL;
    mxArray * WBSum = NULL;
    mxArray * perc3 = NULL;
    mxArray * perc2 = NULL;
    mxArray * perc1 = NULL;
    mxArray * roff = NULL;
    mxArray * hp = NULL;
    mxArray * Fcum = NULL;
    mxArray * flim = NULL;
    mxArray * f = NULL;
    mxArray * tp = NULL;
    mxArray * ttp = NULL;
    mxArray * Fprev = NULL;
    mxArray * F = NULL;
    mxArray * ponding = NULL;
    mxArray * curve_area = NULL;
    mxArray * theta_CZ = NULL;
    mxArray * theta_TZ = NULL;
    mxArray * theta_RZ = NULL;
    mxArray * t = NULL;
    mxArray * Storm_gap = NULL;
    mxArray * A_perv = NULL;
    mxArray * hp_0 = NULL;
    mxArray * pi = NULL;
    mxArray * Co = NULL;
    mxArray * hmax = NULL;
    mxArray * udiam = NULL;
    mxArray * Umax = NULL;
    mxArray * Rp = NULL;
    mxArray * Ri = NULL;
    mxArray * Rb = NULL;
    mxArray * Sd_perv = NULL;
    mxArray * theta_init = NULL;
    mxArray * i = NULL;
    mxArray * hin = NULL;
    mxArray * thick = NULL;
    mxArray * KSAT = NULL;
    mxArray * NRES = NULL;
    mxArray * NSAT = NULL;
    mxArray * hwf = NULL;
    mxArray * NREScz = NULL;
    mxArray * NSATcz = NULL;
    mxArray * h0 = NULL;
    mxArray * N = NULL;
    mxArray * Scz = NULL;
    mxArray * stime = NULL;
    mxArray * dt = NULL;
    mxArray * time = NULL;
    mxArray * nsfactor = NULL;
    mxArray * nwiltfactor = NULL;
    mxArray * rzwilt = NULL;
    mxArray * b = NULL;
    mxArray * he = NULL;
    mxArray * Ss = NULL;
    mxArray * AA = NULL;
    mxArray * NRESrz = NULL;
    mxArray * NSATrz = NULL;
    mxArray * NRESsz = NULL;
    mxArray * NSATsz = NULL;
    mxArray * soilsdata = NULL;
    mxArray * Cpan = NULL;
    mxArray * Evap_perv = NULL;
    mxArray * Evap_imp = NULL;
    mxArray * Si_rg = NULL;
    mxArray * Si_perv = NULL;
    mxArray * Sd_imp = NULL;
    mxArray * g = NULL;
    mxArray * ans = NULL;
    mclCopyArray(&nprecip);
    mclCopyArray(&Sd);
    mclCopyArray(&avet);
    mclCopyArray(&rundays);
    mclCopyArray(&outpre);
    mclCopyArray(&Summary);
    mclCopyArray(&Record);
    mclCopyArray(&dmin);
    mclCopyArray(&Srz);
    mclCopyArray(&St);
    mclCopyArray(&Rsoilnumber);
    mclCopyArray(&Ssoilnumber);
    mclCopyArray(&Nsoilnumber);
    mclCopyArray(&Krz);
    mclCopyArray(&Kst);
    mclCopyArray(&Kcz);
    mclCopyArray(&Urate);
    mclCopyArray(&CN);
    mclCopyArray(&EIA);
    mclCopyArray(&Arg);
    mclCopyArray(&TribA);
    /*
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %  BASIC PROGRAM INFORMATION                                              
     * %
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %
     * %  FILE NAME:
     * %  recarga25.m
     * %
     * %  GENERAL DESCRIPTION:
     * %  This program computes the values of the water balance terms (inflow,
     * %  infiltration, runoff, evapotranspiration, andrecharge) for a raingarden 
     * %  or area of lawn receiving runoff from animpervious surface, such as a 
     * %  roof.
     * %  *MULTILAYER GREEN&AMPT!!!
     * %
     * %  ASSUMPTIONS:
     * %  1. Conditions corresponding to a summer rainy season (no snow).  
     * %  2. Vertical, one-dimensional infiltration in a uniform, partially 
     * %     saturated soil.
     * %  3. hwf (wetting front suction) is constant, and equal to root zone soil value
     * %     [this could be modified: follow wf and update hwf in time, but
     * %      it has been shown GA is not very sensitive to hwf parameter]
     * %
     * %  CALCULATION PROCEDURES:
     * %  1. Infiltration is calculated according to Green-Ampt.
     * %  2. Evaporation/evapotranspiration calculated using pan evaporation
     * %  records
     * %
     * %  REFERENCES:
     * %  1. van Genuchten, 1980  (soil characteristic functions)
     * %
     * %  REQUIRED FUNCTION M-FILES:
     * %  1. storageout.m
     * %  2. stormcheck.m, eventinfo.m
     * %  3. iga.m
     * %  4. relcond.m
     * %
     * %  HOW TO USE THE PROGRAM:
     * %  For each model run, the user will want to adjust the parameters and
     * %  input file in order to obtain the information he desires.  Adjust the
     * %  code at the appropriate locations prior to each model run.  Make sure
     * %  all units are consistent (see comments within code).  The input file
     * %  must be a text file containing two tab-delimited columns, time [h] and
     * %  precipitation depth [cm].  The program assumes the first row of each 
     * %  column is the column name, and eliminates it from calculations.  The
     * %  second must be 0 hours, O precipitation depth.
     * %
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %  PARAMETERS                                                             
     * %
     * format long
     */
    mlfFormat(_mxarray0_, NULL);
    /*
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * % TRIBUTARY BASIN PARAMETERS
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * g=9.8;		% acceleration of gravity [m/s2]
     */
    mlfAssign(&g, _mxarray2_);
    /*
     * Sd_imp = Sd;           % Depression storage depth of impervious surface [cm]
     */
    mlfAssign(&Sd_imp, mclVa(Sd, "Sd"));
    /*
     * Si_perv = Sd;          % Interception storage depth of pervious surface [cm]
     */
    mlfAssign(&Si_perv, mclVa(Sd, "Sd"));
    /*
     * Si_rg = Sd;            % Interception storage depth of pervious surface [cm]
     */
    mlfAssign(&Si_rg, mclVa(Sd, "Sd"));
    /*
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * % EVAPORATION PARAMETERS
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * 
     * Evap_imp = 10/10/24;    % Average rate of evaporation from impervious surface [cm/h]
     */
    mlfAssign(&Evap_imp, _mxarray3_);
    /*
     * Evap_perv = avet/24;  % Average rate of evaporation from pervious surface [cm/h]
     */
    mlfAssign(&Evap_perv, mclMrdivide(mclVa(avet, "avet"), _mxarray4_));
    /*
     * Cpan = .75;            % {DBA} Pan Evaporation correction coefficient
     */
    mlfAssign(&Cpan, _mxarray5_);
    /*
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * % SOIL SENSITIVITY PARAMETERS:
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * soilsdata=dlmread('soilsDNR.txt','\t',[0 1 10 4]);
     */
    mlfAssign(
      &soilsdata, mlfDlmread(_mxarray6_, _mxarray8_, _mxarray10_, NULL, NULL));
    /*
     * %soilnames=textread('soilslist.txt','%s','delimiter','\n','whitespace','');
     * %STexture = soilnames(Ssoilnumber,1);
     * %RTexture = soilnames(Rsoilnumber,1);
     * NSATsz = soilsdata(Ssoilnumber,1);  %Porosity of control zone
     */
    mlfAssign(
      &NSATsz,
      mclArrayRef2(
        mclVv(soilsdata, "soilsdata"),
        mclVa(Ssoilnumber, "Ssoilnumber"),
        _mxarray12_));
    /*
     * NRESsz = soilsdata(Ssoilnumber,1)-soilsdata(Ssoilnumber,2); %Residual water content, NSAT - NEFF
     */
    mlfAssign(
      &NRESsz,
      mclMinus(
        mclArrayRef2(
          mclVv(soilsdata, "soilsdata"),
          mclVa(Ssoilnumber, "Ssoilnumber"),
          _mxarray12_),
        mclArrayRef2(
          mclVv(soilsdata, "soilsdata"),
          mclVa(Ssoilnumber, "Ssoilnumber"),
          _mxarray13_)));
    /*
     * NSATrz = soilsdata(Rsoilnumber,1);  %Porosity of control zone
     */
    mlfAssign(
      &NSATrz,
      mclArrayRef2(
        mclVv(soilsdata, "soilsdata"),
        mclVa(Rsoilnumber, "Rsoilnumber"),
        _mxarray12_));
    /*
     * NRESrz = soilsdata(Rsoilnumber,1)-soilsdata(Rsoilnumber,2); %Residual water content, NSAT - NEFF
     */
    mlfAssign(
      &NRESrz,
      mclMinus(
        mclArrayRef2(
          mclVv(soilsdata, "soilsdata"),
          mclVa(Rsoilnumber, "Rsoilnumber"),
          _mxarray12_),
        mclArrayRef2(
          mclVv(soilsdata, "soilsdata"),
          mclVa(Rsoilnumber, "Rsoilnumber"),
          _mxarray13_)));
    /*
     * AA = [0.01 0.01 0.01];				% alpha [cm-1]
     */
    mlfAssign(&AA, _mxarray14_);
    /*
     * Ss=[1e-6 1e-6 1e-6];   % specific storage [cm-1]
     */
    mlfAssign(&Ss, _mxarray16_);
    /*
     * he=-1./AA;   % air entry suction [cm]
     */
    mlfAssign(&he, mclRdivide(_mxarray18_, mclVv(AA, "AA")));
    /*
     * b = 3.1;                         %{DBA} Pore size distribution index, from C&N
     */
    mlfAssign(&b, _mxarray19_);
    /*
     * rzwilt = 0.1;                       %{DBA} Root zone wilting point, sandy loam from C&N
     */
    mlfAssign(&rzwilt, _mxarray20_);
    /*
     * nwiltfactor=1.21; %{DBA} for estimation of near wilting point
     */
    mlfAssign(&nwiltfactor, _mxarray21_);
    /*
     * nsfactor=0.95;  % for tracking near saturated conditions in the root zone
     */
    mlfAssign(&nsfactor, _mxarray22_);
    /*
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * % MODELING PARAMETERS
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * % TEMPORAL:
     * 
     * time=rundays*24*60;			% min in data file
     */
    mlfAssign(
      &time,
      mclMtimes(mclMtimes(mclVa(rundays, "rundays"), _mxarray4_), _mxarray23_));
    /*
     * dt=15;						% TIME STEP [min]
     */
    mlfAssign(&dt, _mxarray24_);
    /*
     * stime=time/60;				% total time of simulation [h]
     */
    mlfAssign(&stime, mclMrdivide(mclVv(time, "time"), _mxarray23_));
    /*
     * dt=dt/60;					% TIME STEP [h]
     */
    mlfAssign(&dt, mclMrdivide(mclVv(dt, "dt"), _mxarray23_));
    /*
     * 
     * Scz = 10;             %Control Zone Thickness
     */
    mlfAssign(&Scz, _mxarray25_);
    /*
     * N  = [1.6 1.378 1.242];
     */
    mlfAssign(&N, _mxarray26_);
    /*
     * h0=[-2.4246e-004 -4.7992e-007 -1.2314e-010];
     */
    mlfAssign(&h0, _mxarray28_);
    /*
     * %NTexture = soilnames(Nsoilnumber,1);
     * NSATcz = soilsdata(Nsoilnumber,1);  %Porosity of control zone
     */
    mlfAssign(
      &NSATcz,
      mclArrayRef2(
        mclVv(soilsdata, "soilsdata"),
        mclVa(Nsoilnumber, "Nsoilnumber"),
        _mxarray12_));
    /*
     * NREScz = soilsdata(Nsoilnumber,1)-soilsdata(Nsoilnumber,2); %Residual water content, NSAT - NEFF
     */
    mlfAssign(
      &NREScz,
      mclMinus(
        mclArrayRef2(
          mclVv(soilsdata, "soilsdata"),
          mclVa(Nsoilnumber, "Nsoilnumber"),
          _mxarray12_),
        mclArrayRef2(
          mclVv(soilsdata, "soilsdata"),
          mclVa(Nsoilnumber, "Nsoilnumber"),
          _mxarray13_)));
    /*
     * if St==0    %Allow for no storage by dividing control zone into two equal layers
     */
    if (mclEqBool(mclVa(St, "St"), _mxarray30_)) {
        /*
         * St=5;
         */
        mlfAssign(&St, _mxarray31_);
        /*
         * Scz=5;
         */
        mlfAssign(&Scz, _mxarray31_);
        /*
         * Kst = Kcz;
         */
        mlfAssign(&Kst, mclVa(Kcz, "Kcz"));
        /*
         * NSATsz = NSATcz;
         */
        mlfAssign(&NSATsz, mclVv(NSATcz, "NSATcz"));
        /*
         * NRESsz = NREScz;
         */
        mlfAssign(&NRESsz, mclVv(NREScz, "NREScz"));
        /*
         * N  = [1.6 1.242 1.242];
         */
        mlfAssign(&N, _mxarray32_);
        /*
         * h0=[-2.4246e-004 -1.2314e-010 -1.2314e-010];
         */
        mlfAssign(&h0, _mxarray34_);
    /*
     * end
     */
    }
    /*
     * hwf=3*(2+3*(N(1)-1))/(1+3*(N(1)-1))*he(1)/2;%11;% wetting front suction of upper layer, Rawls&Brakensiek
     */
    mlfAssign(
      &hwf,
      mclMrdivide(
        mclMtimes(
          mclMrdivide(
            mclMtimes(
              _mxarray36_,
              mclPlus(
                _mxarray13_,
                mclMtimes(
                  _mxarray36_,
                  mclMinus(mclIntArrayRef1(mclVv(N, "N"), 1), _mxarray12_)))),
            mclPlus(
              _mxarray12_,
              mclMtimes(
                _mxarray36_,
                mclMinus(mclIntArrayRef1(mclVv(N, "N"), 1), _mxarray12_)))),
          mclIntArrayRef1(mclVv(he, "he"), 1)),
        _mxarray13_));
    /*
     * NSAT = [NSATrz NSATsz NSATcz];				% water content at saturation
     */
    mlfAssign(
      &NSAT,
      mlfHorzcat(
        mclVv(NSATrz, "NSATrz"),
        mclVv(NSATsz, "NSATsz"),
        mclVv(NSATcz, "NSATcz"),
        NULL));
    /*
     * NRES = [NRESrz NRESsz NREScz];				% residual water content
     */
    mlfAssign(
      &NRES,
      mlfHorzcat(
        mclVv(NRESrz, "NRESrz"),
        mclVv(NRESsz, "NRESsz"),
        mclVv(NREScz, "NREScz"),
        NULL));
    /*
     * KSAT = [Krz Kst Kcz];	% sat. hydraulic conductivity [cm/h]
     */
    mlfAssign(
      &KSAT,
      mlfHorzcat(
        mclVa(Krz, "Krz"), mclVa(Kst, "Kst"), mclVa(Kcz, "Kcz"), NULL));
    /*
     * thick = [Srz St Scz];		% thickness of RZ, TZ, CZ layers [cm]
     */
    mlfAssign(
      &thick,
      mlfHorzcat(mclVa(Srz, "Srz"), mclVa(St, "St"), mclVv(Scz, "Scz"), NULL));
    /*
     * 
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * % INITIAL CONDITIONS:
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * hin=-100;					% [cm]
     */
    mlfAssign(&hin, _mxarray37_);
    /*
     * for i=1:3
     */
    {
        int v_ = mclForIntStart(1);
        int e_ = 3;
        if (v_ > e_) {
            mlfAssign(&i, _mxarray38_);
        } else {
            /*
             * theta_init(i)=THETAH(hin,AA(i),N(i),NSAT(i),NRES(i),Ss(i),h0(i)); 
             * %set init. moisture for initial pressure head
             * %theta_init(i)=NRES(i);  %For use to begin simulation at field capacity
             * end
             */
            for (; ; ) {
                mclIntArrayAssign1(
                  &theta_init,
                  mlfThetah(
                    mclVv(hin, "hin"),
                    mclIntArrayRef1(mclVv(AA, "AA"), v_),
                    mclIntArrayRef1(mclVv(N, "N"), v_),
                    mclIntArrayRef1(mclVv(NSAT, "NSAT"), v_),
                    mclIntArrayRef1(mclVv(NRES, "NRES"), v_),
                    mclIntArrayRef1(mclVv(Ss, "Ss"), v_),
                    mclIntArrayRef1(mclVv(h0, "h0"), v_)),
                  v_);
                if (v_ == e_) {
                    break;
                }
                ++v_;
            }
            mlfAssign(&i, mlfScalar(v_));
        }
    }
    /*
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * Sd_perv = 2.54*(1000/CN-10);     % Depression storage depth of pervious surface CN=80 [cm]
     */
    mlfAssign(
      &Sd_perv,
      mclMtimes(
        _mxarray39_,
        mclMinus(mclMrdivide(_mxarray40_, mclVa(CN, "CN")), _mxarray25_)));
    /*
     * Rb = Arg/TribA;  %Ratio of Bioretention to Total Area
     */
    mlfAssign(&Rb, mclMrdivide(mclVa(Arg, "Arg"), mclVa(TribA, "TribA")));
    /*
     * Ri = EIA/100/Rb;  %Ratio of impervious area to facility area
     */
    mlfAssign(
      &Ri,
      mclMrdivide(
        mclMrdivide(mclVa(EIA, "EIA"), _mxarray41_), mclVv(Rb, "Rb")));
    /*
     * Rp = (1-EIA/100)/Rb; %Ratio of pervious area to facility area
     */
    mlfAssign(
      &Rp,
      mclMrdivide(
        mclMinus(_mxarray12_, mclMrdivide(mclVa(EIA, "EIA"), _mxarray41_)),
        mclVv(Rb, "Rb")));
    /*
     * 
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %  UNDERDRAIN SIZING
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * Umax = (Urate/100)*(Arg/3600); %Max. Underdrain flow [m3/sec]
     */
    mlfAssign(
      &Umax,
      mclMtimes(
        mclMrdivide(mclVa(Urate, "Urate"), _mxarray41_),
        mclMrdivide(mclVa(Arg, "Arg"), _mxarray42_)));
    /*
     * if Umax == 0
     */
    if (mclEqBool(mclVv(Umax, "Umax"), _mxarray30_)) {
        /*
         * udiam=0;
         */
        mlfAssign(&udiam, _mxarray30_);
    /*
     * end
     */
    }
    /*
     * hmax = (dmin + Srz)/100;  %Max. Head on underdrain [m]
     */
    mlfAssign(
      &hmax,
      mclMrdivide(
        mclPlus(mclVa(dmin, "dmin"), mclVa(Srz, "Srz")), _mxarray41_));
    /*
     * Co = 0.62; g=9.81; pi=3.14159; %constants
     */
    mlfAssign(&Co, _mxarray43_);
    mlfAssign(&g, _mxarray44_);
    mlfAssign(&pi, _mxarray45_);
    /*
     * udiam = 2*sqrt(Umax/(Co*pi*sqrt(2*g*hmax)));  %underdrain diameter [m]
     */
    mlfAssign(
      &udiam,
      mclMtimes(
        _mxarray13_,
        mlfSqrt(
          mclMrdivide(
            mclVv(Umax, "Umax"),
            mclMtimes(
              mclMtimes(mclVv(Co, "Co"), mclVv(pi, "pi")),
              mlfSqrt(
                mclMtimes(
                  mclMtimes(_mxarray13_, mclVv(g, "g")),
                  mclVv(hmax, "hmax"))))))));
    /*
     * 
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %  PERVIOUS SURFACE              %
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %
     * hp_0 = 0;               % Ponded depth at impervious surface initialized [cm]
     */
    mlfAssign(&hp_0, _mxarray30_);
    /*
     * A_perv = Arg;           % Area of pervious surface (raingarden or lawn area) [m2]
     */
    mlfAssign(&A_perv, mclVa(Arg, "Arg"));
    /*
     * %
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %  STORM SIMULATION              %
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %
     * Storm_gap = 6;  % Time with no rain defining gap between storms [h]
     */
    mlfAssign(&Storm_gap, _mxarray46_);
    /*
     * t=0;
     */
    mlfAssign(&t, _mxarray30_);
    /*
     * % dt = 5/60;      % Time step used in model (use divisor of data interval) [h]
     * % stime=4;        % [h], end of simulation
     * %
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %  SOIL                           %
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %
     * % theta_init=NSAT;  % Initial soil water moisture for Green-Ampt
     * %
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %  INTIIALIZE VARIABLES AND COMPUTE VALUES                                
     * %
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %  INFILTRATION                  %
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %
     * 
     * theta_RZ = theta_init(1);         % Soil moisture of root zone initialized
     */
    mlfAssign(&theta_RZ, mclIntArrayRef1(mclVv(theta_init, "theta_init"), 1));
    /*
     * %theta_RZ = .2;         % Use to match event model
     * theta_TZ=theta_init(2); theta_CZ=theta_init(3);
     */
    mlfAssign(&theta_TZ, mclIntArrayRef1(mclVv(theta_init, "theta_init"), 2));
    mlfAssign(&theta_CZ, mclIntArrayRef1(mclVv(theta_init, "theta_init"), 3));
    /*
     * curve_area = 0;             % Area under storm event intensity curve initialized [cm]
     */
    mlfAssign(&curve_area, _mxarray30_);
    /*
     * ponding = 0;                % Ponding indicator initialized to false [0 or 1]
     */
    mlfAssign(&ponding, _mxarray30_);
    /*
     * F = 0;                      % Cumulative infiltration depth initialized [cm]
     */
    mlfAssign(&F, _mxarray30_);
    /*
     * Fprev=0;
     */
    mlfAssign(&Fprev, _mxarray30_);
    /*
     * ttp = 0;                    % Total time ponded initialized [h]
     */
    mlfAssign(&ttp, _mxarray30_);
    /*
     * tp = 0;                     % time of ponding for storm [h]
     */
    mlfAssign(&tp, _mxarray30_);
    /*
     * %
     * %
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %  WATER BALANCE                 %
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %
     * f = 0;          % Infiltration rate initialized [cm/h]
     */
    mlfAssign(&f, _mxarray30_);
    /*
     * flim = 0;       % Used when saturated conditions impose limits on f.  [cm/h]
     */
    mlfAssign(&flim, _mxarray30_);
    /*
     * Fcum=0;         % Cumulative infiltration depth [cm]
     */
    mlfAssign(&Fcum, _mxarray30_);
    /*
     * hp=hp_0;
     */
    mlfAssign(&hp, mclVv(hp_0, "hp_0"));
    /*
     * roff = 0;       % Runoff rate initialized [cm/h]
     */
    mlfAssign(&roff, _mxarray30_);
    /*
     * perc1=0; perc2=0;perc3(:,1)=0;       % Percolation rate initialized [cm/h]
     */
    mlfAssign(&perc1, _mxarray30_);
    mlfAssign(&perc2, _mxarray30_);
    mclArrayAssign2(&perc3, _mxarray30_, mlfCreateColonIndex(), _mxarray12_);
    /*
     * WBSum(1,:)=[0 0 0 0 0 0 0 0 0 0 0];                       %Water Budget Summary initialized
     */
    mclArrayAssign2(&WBSum, _mxarray47_, _mxarray12_, mlfCreateColonIndex());
    /*
     * Inflow(:,1)=0; udrain(:,1)= 0; infilt(:,1)= 0; Runoff(:,1)= 0; 
     */
    mclArrayAssign2(&Inflow, _mxarray30_, mlfCreateColonIndex(), _mxarray12_);
    mclArrayAssign2(&udrain, _mxarray30_, mlfCreateColonIndex(), _mxarray12_);
    mclArrayAssign2(&infilt, _mxarray30_, mlfCreateColonIndex(), _mxarray12_);
    mclArrayAssign2(&Runoff, _mxarray30_, mlfCreateColonIndex(), _mxarray12_);
    /*
     * EvapoT(:,1)= 0; Recharge(:,1)=0;
     */
    mclArrayAssign2(&EvapoT, _mxarray30_, mlfCreateColonIndex(), _mxarray12_);
    mclArrayAssign2(&Recharge, _mxarray30_, mlfCreateColonIndex(), _mxarray12_);
    /*
     * x= 1; s = 1;
     */
    mlfAssign(&x, _mxarray12_);
    mlfAssign(&s, _mxarray12_);
    /*
     * perc1lim=0; perc2lim=0; perc3lim=0; % Used when saturated condtns impose limits on perc. [cm/h]
     */
    mlfAssign(&perc1lim, _mxarray30_);
    mlfAssign(&perc2lim, _mxarray30_);
    mlfAssign(&perc3lim, _mxarray30_);
    /*
     * V_inp = 0;      % Inflow volume initialized [m^3]
     */
    mlfAssign(&V_inp, _mxarray30_);
    /*
     * V_inf = 0;      % Infiltrated volume initialized [m^3]
     */
    mlfAssign(&V_inf, _mxarray30_);
    /*
     * V_roff = 0;     % Runoff volume initialized [m^3]
     */
    mlfAssign(&V_roff, _mxarray30_);
    /*
     * V_rech = 0;     % Recharge volume initialized [m^3]
     */
    mlfAssign(&V_rech, _mxarray30_);
    /*
     * V_ud = 0;       % Underdrain volume initialized [m^3]
     */
    mlfAssign(&V_ud, _mxarray30_);
    /*
     * V_et = 0;       % Evapotranspiration volume initialized [m^3]
     */
    mlfAssign(&V_et, _mxarray30_);
    /*
     * V_pond = 0;     % Ponding volume initialized [m^3]  %TAKE OUT
     */
    mlfAssign(&V_pond, _mxarray30_);
    /*
     * tctpmax=0; tstponded=0; tctponded=0;    % initial counters for ponded time
     */
    mlfAssign(tctpmax, _mxarray30_);
    mlfAssign(tstponded, _mxarray30_);
    mlfAssign(&tctponded, _mxarray30_);
    /*
     * tctoverflow=0;  %initial counters for overflow time
     */
    mlfAssign(&tctoverflow, _mxarray30_);
    /*
     * % track near saturated total time (>nsfactor):
     * tstrznsat=0; tsttznsat=0; tctrznsatmax=0; tcttznsatmax=0; tctrznsat=0; 
     */
    mlfAssign(&tstrznsat, _mxarray30_);
    mlfAssign(&tsttznsat, _mxarray30_);
    mlfAssign(&tctrznsatmax, _mxarray30_);
    mlfAssign(&tcttznsatmax, _mxarray30_);
    mlfAssign(&tctrznsat, _mxarray30_);
    /*
     * tcttznsat=0;
     */
    mlfAssign(&tcttznsat, _mxarray30_);
    /*
     * tstrzwilt=0; tctrzwilt=0; tctrzwiltmax=0; %{DBA}
     */
    mlfAssign(&tstrzwilt, _mxarray30_);
    mlfAssign(&tctrzwilt, _mxarray30_);
    mlfAssign(&tctrzwiltmax, _mxarray30_);
    /*
     * inst=0; nstimes(1,:)=[0 0]; % counter for vector 'nstimes' that will store nsat times for simulation
     */
    mlfAssign(&inst, _mxarray30_);
    mclArrayAssign2(&nstimes, _mxarray49_, _mxarray12_, mlfCreateColonIndex());
    /*
     * inpt=0; ptimes(1,:)=[0 0]; % counter for vector 'ptimes' that will store ponded times for simulation
     */
    mlfAssign(&inpt, _mxarray30_);
    mclArrayAssign2(&ptimes, _mxarray49_, _mxarray12_, mlfCreateColonIndex());
    /*
     * inot=0; oftimes(1,:)=[0 0]; prevof=0; % counter for vector 'oftimes' that will store overflow times for simulation
     */
    mlfAssign(inot, _mxarray30_);
    mclArrayAssign2(&oftimes, _mxarray49_, _mxarray12_, mlfCreateColonIndex());
    mlfAssign(&prevof, _mxarray30_);
    /*
     * inwt=0; wtimes(1,:)=[0 0]; % {DBA} counter for vector 'wtimes' that will store wilted times for simulation
     */
    mlfAssign(&inwt, _mxarray30_);
    mclArrayAssign2(&wtimes, _mxarray49_, _mxarray12_, mlfCreateColonIndex());
    /*
     * %
     * %
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %  DATA RETRIEVAL AND PREPROCESSING                                       
     * %
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %
     * % Read data file from working directory.  The file must be a tab-delimited
     * % text file with two columns.  The first column must be time [h] and the 
     * % second column must be precipitation depth [cm].  The first row of each
     * % column should contain the title of the column, as the model will not look
     * % at that row.
     * %
     * qs=nprecip(1,2);    %first value for surface flux [cm h-1]
     */
    mlfAssign(&qs, mclIntArrayRef2(mclVa(nprecip, "nprecip"), 1, 2));
    /*
     * es=nprecip(1,3);    %first value for surface flux
     */
    mlfAssign(&es, mclIntArrayRef2(mclVa(nprecip, "nprecip"), 1, 3));
    /*
     * pdfile = nprecip; %dlmread('lawntest.txt','\t',1,0); % Read data file [h cm]
     */
    mlfAssign(&pdfile, mclVa(nprecip, "nprecip"));
    /*
     * %
     * %
     * len = size(pdfile,1); % Calculate number of data points %PUT BACK IN
     */
    mlfAssign(
      &len, mlfSize(mclValueVarargout(), mclVv(pdfile, "pdfile"), _mxarray12_));
    /*
     * if stime>pdfile(len,1);stime=pdfile(len,1);end; % Set max simulation time as final time in data set [h] % PUT BACK IN
     */
    if (mclGtBool(
          mclVv(stime, "stime"),
          mclArrayRef2(
            mclVv(pdfile, "pdfile"), mclVv(len, "len"), _mxarray12_))) {
        mlfAssign(
          &stime,
          mclArrayRef2(
            mclVv(pdfile, "pdfile"), mclVv(len, "len"), _mxarray12_));
    }
    /*
     * %
     * %
     * % Linearly interpolate the data set using the model time step.  This will
     * % produce a new two-column array with values of time and precipitation 
     * % depth at each multiple of the time step.
     * %
     * tk = 0:dt:stime; tk = tk'; % Create a column vector with multiples of dt
     */
    mlfAssign(
      &tk, mlfColon(_mxarray30_, mclVv(dt, "dt"), mclVv(stime, "stime")));
    mlfAssign(&tk, mlfCtranspose(mclVv(tk, "tk")));
    /*
     * CummRain = cumsum(pdfile(:,2));  % Cummulative Rainfall
     */
    mlfAssign(
      &CummRain,
      mlfCumsum(
        mclArrayRef2(
          mclVv(pdfile, "pdfile"), mlfCreateColonIndex(), _mxarray13_),
        NULL));
    /*
     * DARain = interp1(pdfile(:,1),CummRain,tk,'*cubic'); % Disaggregated Rainfall
     */
    mlfAssign(
      &DARain,
      mlfInterp1(
        mclArrayRef2(
          mclVv(pdfile, "pdfile"), mlfCreateColonIndex(), _mxarray12_),
        mclVv(CummRain, "CummRain"),
        mclVv(tk, "tk"),
        _mxarray51_,
        NULL));
    /*
     * %intensity = diff(DARain)/dt;  % Return to intensity column vector
     * %intensity = [DARain(1);intensity];  % Offset data by 1
     * intensity = interp1(pdfile(:,1),pdfile(:,2),tk,'*linear'); % Intensity vector
     */
    mlfAssign(
      &intensity,
      mlfInterp1(
        mclArrayRef2(
          mclVv(pdfile, "pdfile"), mlfCreateColonIndex(), _mxarray12_),
        mclArrayRef2(
          mclVv(pdfile, "pdfile"), mlfCreateColonIndex(), _mxarray13_),
        mclVv(tk, "tk"),
        _mxarray53_,
        NULL));
    /*
     * %Pevap = interp1(pdfile(:,1),pdfile(:,3),tk,'*linear'); %{DBA} Pot. Evaporation
     * intfile = [tk intensity]; % Interpolated time/intensity array [h cm/h]
     */
    mlfAssign(
      &intfile,
      mlfHorzcat(mclVv(tk, "tk"), mclVv(intensity, "intensity"), NULL));
    /*
     * Storm_gap = Storm_gap/dt; % Convert storm gap from [h] to # time steps
     */
    mlfAssign(
      &Storm_gap, mclMrdivide(mclVv(Storm_gap, "Storm_gap"), mclVv(dt, "dt")));
    /*
     * %
     * % Take the interpolated time/intensity array and determine which time steps
     * % correspond to storm events.  This is done using the stormcheck function.
     * % The function accepts a two-column time/intensity array and the storm gap.
     * % It adds a third column indicating whether the time step is part of a storm
     * % event [0 or 1] and a fourth column indicating the length of the event,
     * % shown at the first time step in an event.
     * %
     * datafile = stormcheck(intfile,Storm_gap); % Obtain storm information
     */
    mlfAssign(
      &datafile,
      mlfStormcheck(mclVv(intfile, "intfile"), mclVv(Storm_gap, "Storm_gap")));
    /*
     * Dprecip = sum(nprecip(:,2));
     */
    mlfAssign(
      &Dprecip,
      mlfSum(
        mclArrayRef2(
          mclVa(nprecip, "nprecip"), mlfCreateColonIndex(), _mxarray13_),
        NULL));
    /*
     * %
     * % Remove the depression storage in impervious surface from runon:
     * %
     * irain = [tk dt*datafile(:,2)*Ri];    % Compute total ppt on impervious area
     */
    mlfAssign(
      &irain,
      mlfHorzcat(
        mclVv(tk, "tk"),
        mclMtimes(
          mclMtimes(
            mclVv(dt, "dt"),
            mclArrayRef2(
              mclVv(datafile, "datafile"), mlfCreateColonIndex(), _mxarray13_)),
          mclVv(Ri, "Ri")),
        NULL));
    /*
     * evap = Evap_imp*dt;                     % Compute the evaporated depth per time step [cm]
     */
    mlfAssign(&evap, mclMtimes(mclVv(Evap_imp, "Evap_imp"), mclVv(dt, "dt")));
    /*
     * %runon = storageout(irain,Sd_imp,evap);  % Remove water stored in impervious
     * runon1 = storageout(irain,Sd_imp*Ri,evap*Ri,datafile);  % Remove water stored in impervious
     */
    mlfAssign(
      &runon1,
      mlfStorageout(
        mclVv(irain, "irain"),
        mclMtimes(mclVv(Sd_imp, "Sd_imp"), mclVv(Ri, "Ri")),
        mclMtimes(mclVv(evap, "evap"), mclVv(Ri, "Ri")),
        mclVv(datafile, "datafile")));
    /*
     * runon1(:,2) = runon1(:,2) / dt; 
     */
    mclArrayAssign2(
      &runon1,
      mclMrdivide(
        mclArrayRef2(
          mclVv(runon1, "runon1"), mlfCreateColonIndex(), _mxarray13_),
        mclVv(dt, "dt")),
      mlfCreateColonIndex(),
      _mxarray13_);
    /*
     * Vimp = sum(runon1(:,2))* dt;
     */
    mlfAssign(
      &Vimp,
      mclMtimes(
        mlfSum(
          mclArrayRef2(
            mclVv(runon1, "runon1"), mlfCreateColonIndex(), _mxarray13_),
          NULL),
        mclVv(dt, "dt")));
    /*
     * if Ri == 0;
     */
    if (mclEqBool(mclVv(Ri, "Ri"), _mxarray30_)) {
        /*
         * Dimp = 0;
         */
        mlfAssign(Dimp, _mxarray30_);
    /*
     * else 
     */
    } else {
        /*
         * Dimp = Vimp/Ri;
         */
        mlfAssign(Dimp, mclMrdivide(mclVv(Vimp, "Vimp"), mclVv(Ri, "Ri")));
    /*
     * end
     */
    }
    /*
     * %
     * % Remove the interception storage in pervious surface from rain:
     * % 
     * rgrain = [tk dt*datafile(:,2)];          % Compute the total ppt on pervious area
     */
    mlfAssign(
      &rgrain,
      mlfHorzcat(
        mclVv(tk, "tk"),
        mclMtimes(
          mclVv(dt, "dt"),
          mclArrayRef2(
            mclVv(datafile, "datafile"), mlfCreateColonIndex(), _mxarray13_)),
        NULL));
    /*
     * rain = storageout(rgrain,Si_rg,evap,datafile);  % Remove water stored in vegetation
     */
    mlfAssign(
      &rain,
      mlfStorageout(
        mclVv(rgrain, "rgrain"),
        mclVv(Si_rg, "Si_rg"),
        mclVv(evap, "evap"),
        mclVv(datafile, "datafile")));
    /*
     * rain (:,2) = rain(:,2) / dt;
     */
    mclArrayAssign2(
      &rain,
      mclMrdivide(
        mclArrayRef2(mclVv(rain, "rain"), mlfCreateColonIndex(), _mxarray13_),
        mclVv(dt, "dt")),
      mlfCreateColonIndex(),
      _mxarray13_);
    /*
     * % Calc runoff from pervious areas using NRCS Method
     * Gap_P = 24;  % 24 hour gap for pervious areas
     */
    mlfAssign(&Gap_P, _mxarray4_);
    /*
     * Gap_P = Gap_P/dt; % Convert storm gap from [h] to # time steps
     */
    mlfAssign(&Gap_P, mclMrdivide(mclVv(Gap_P, "Gap_P"), mclVv(dt, "dt")));
    /*
     * pdatafile = stormcheck(intfile,Gap_P); % Obtain storm information
     */
    mlfAssign(
      &pdatafile,
      mlfStormcheck(mclVv(intfile, "intfile"), mclVv(Gap_P, "Gap_P")));
    /*
     * prain=[tk dt*pdatafile(:,2)];
     */
    mlfAssign(
      &prain,
      mlfHorzcat(
        mclVv(tk, "tk"),
        mclMtimes(
          mclVv(dt, "dt"),
          mclArrayRef2(
            mclVv(pdatafile, "pdatafile"), mlfCreateColonIndex(), _mxarray13_)),
        NULL));
    /*
     * if Rp == 0;
     */
    if (mclEqBool(mclVv(Rp, "Rp"), _mxarray30_)) {
        /*
         * runon2 = storageout(prain,0,evap,pdatafile);
         */
        mlfAssign(
          &runon2,
          mlfStorageout(
            mclVv(prain, "prain"),
            _mxarray30_,
            mclVv(evap, "evap"),
            mclVv(pdatafile, "pdatafile")));
    /*
     * else
     */
    } else {
        /*
         * runon2 = scsII(prain,CN,evap,pdatafile);
         */
        mlfAssign(
          &runon2,
          mlfScsii(
            mclVv(prain, "prain"),
            mclVa(CN, "CN"),
            mclVv(evap, "evap"),
            mclVv(pdatafile, "pdatafile")));
    /*
     * end
     */
    }
    /*
     * runon2(:,2) = Rp*runon2(:,2) / dt;
     */
    mclArrayAssign2(
      &runon2,
      mclMrdivide(
        mclMtimes(
          mclVv(Rp, "Rp"),
          mclArrayRef2(
            mclVv(runon2, "runon2"), mlfCreateColonIndex(), _mxarray13_)),
        mclVv(dt, "dt")),
      mlfCreateColonIndex(),
      _mxarray13_);
    /*
     * Vperv = sum(runon2(:,2)) *  dt;
     */
    mlfAssign(
      &Vperv,
      mclMtimes(
        mlfSum(
          mclArrayRef2(
            mclVv(runon2, "runon2"), mlfCreateColonIndex(), _mxarray13_),
          NULL),
        mclVv(dt, "dt")));
    /*
     * if Rp == 0;
     */
    if (mclEqBool(mclVv(Rp, "Rp"), _mxarray30_)) {
        /*
         * Dperv=0;
         */
        mlfAssign(Dperv, _mxarray30_);
    /*
     * else
     */
    } else {
        /*
         * Dperv = Vperv/Rp;
         */
        mlfAssign(Dperv, mclMrdivide(mclVv(Vperv, "Vperv"), mclVv(Rp, "Rp")));
    /*
     * end
     */
    }
    /*
     * Dbasin=(Dprecip-Dimp)*(EIA)+(Dprecip-Dperv)*(1-EIA);
     */
    mlfAssign(
      &Dbasin,
      mclPlus(
        mclMtimes(
          mclMinus(mclVv(Dprecip, "Dprecip"), mclVv(*Dimp, "Dimp")),
          mclVa(EIA, "EIA")),
        mclMtimes(
          mclMinus(mclVv(Dprecip, "Dprecip"), mclVv(*Dperv, "Dperv")),
          mclMinus(_mxarray12_, mclVa(EIA, "EIA")))));
    /*
     * %
     * % Correct input file
     * %
     * runon(:,2) = runon1(:,2) + runon2(:,2);
     */
    mclArrayAssign2(
      &runon,
      mclPlus(
        mclArrayRef2(
          mclVv(runon1, "runon1"), mlfCreateColonIndex(), _mxarray13_),
        mclArrayRef2(
          mclVv(runon2, "runon2"), mlfCreateColonIndex(), _mxarray13_)),
      mlfCreateColonIndex(),
      _mxarray13_);
    /*
     * datafile(:,2)=rain(:,2)+runon(:,2); % Sum effects of impervious & pervious
     */
    mclArrayAssign2(
      &datafile,
      mclPlus(
        mclArrayRef2(mclVv(rain, "rain"), mlfCreateColonIndex(), _mxarray13_),
        mclArrayRef2(
          mclVv(runon, "runon"), mlfCreateColonIndex(), _mxarray13_)),
      mlfCreateColonIndex(),
      _mxarray13_);
    /*
     * %
     * % Create vectors for convenient use of key values
     * %
     * inp = datafile(:,2);            % [cm/h]
     */
    mlfAssign(
      &inp,
      mclArrayRef2(
        mclVv(datafile, "datafile"), mlfCreateColonIndex(), _mxarray13_));
    /*
     * 
     * indstorm = datafile(:,3);       % 0/1
     */
    mlfAssign(
      &indstorm,
      mclArrayRef2(
        mclVv(datafile, "datafile"), mlfCreateColonIndex(), _mxarray36_));
    /*
     * ttstorm = datafile(:,4) * dt;
     */
    mlfAssign(
      &ttstorm,
      mclMtimes(
        mclArrayRef2(
          mclVv(datafile, "datafile"), mlfCreateColonIndex(), _mxarray55_),
        mclVv(dt, "dt")));
    /*
     * %
     * % Create output file for storing results of simulation
     * % delete and open file for output, print simulation parameters:
     * % outfile='lawntest_out.txt';
     * 
     * %
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %  MAIN SIMULATION LOOP                                                   
     * %
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %
     * t=0; n=1;
     */
    mlfAssign(&t, _mxarray30_);
    mlfAssign(&n, _mxarray12_);
    /*
     * % This indicates that the following loop begins at the BEGINNING of the first
     * % time step.
     * tcpu=cputime;    % starts counting CPU time
     */
    mlfAssign(&tcpu, mlfCputime());
    /*
     * %
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * while t < stime % Check that simulation time not passed
     */
    while (mclLtBool(mclVv(t, "t"), mclVv(stime, "stime"))) {
        /*
         * 
         * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         * % CHECK IF STORM STARTS TO GET INTO LOOP:
         * if indstorm(n)==1
         */
        if (mclEqBool(
              mclArrayRef1(mclVv(indstorm, "indstorm"), mclVv(n, "n")),
              _mxarray12_)) {
            /*
             * %(n==1 & indstorm(n)==1) | (indstorm(n)==1 & indstorm(n-1)==0)
             * test2 = indstorm(n); %TAKE OUT
             */
            mlfAssign(
              &test2, mclArrayRef1(mclVv(indstorm, "indstorm"), mclVv(n, "n")));
            /*
             * test3 = indstorm(n-1); %TAKE OUT
             */
            mlfAssign(
              &test3,
              mclArrayRef1(
                mclVv(indstorm, "indstorm"),
                mclMinus(mclVv(n, "n"), _mxarray12_)));
            /*
             * % reset counters before beginning storm loop:
             * te = 0; ponding = 0;
             */
            mlfAssign(&te, _mxarray30_);
            mlfAssign(&ponding, _mxarray30_);
            /*
             * F = 0; curve_area = 0;
             */
            mlfAssign(&F, _mxarray30_);
            mlfAssign(&curve_area, _mxarray30_);
            /*
             * thetainiRZ = theta_RZ; thetainiTZ=theta_TZ; 
             */
            mlfAssign(&thetainiRZ, mclVv(theta_RZ, "theta_RZ"));
            mlfAssign(&thetainiTZ, mclVv(theta_TZ, "theta_TZ"));
            /*
             * thetainiCZ=theta_CZ; % thetaini is initial soil moisture for event
             */
            mlfAssign(&thetainiCZ, mclVv(theta_CZ, "theta_CZ"));
            /*
             * T = ttstorm(n); % Check length of storm, store as T
             */
            mlfAssign(
              &T, mclArrayRef1(mclVv(ttstorm, "ttstorm"), mclVv(n, "n")));
            /*
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * while (indstorm(n)==1 | ponding==1) & t <= stime   % INFILTRATION/SUPPLY LOOP%%%%%
             */
            for (;;) {
                mxArray * a_
                  = mclInitialize(
                      mclEq(
                        mclArrayRef1(
                          mclVv(indstorm, "indstorm"), mclVv(n, "n")),
                        _mxarray12_));
                if (mlfTobool(a_)) {
                    mlfAssign(&a_, mlfScalar(1));
                } else {
                    mlfAssign(
                      &a_,
                      mclOr(a_, mclEq(mclVv(ponding, "ponding"), _mxarray12_)));
                }
                if (mlfTobool(a_)
                    && mlfTobool(
                         mclAnd(
                           a_, mclLe(mclVv(t, "t"), mclVv(stime, "stime"))))) {
                    mxDestroyArray(a_);
                } else {
                    mxDestroyArray(a_);
                    break;
                }
                /*
                 * %(te < T | ponding==1) & t <= stime   % INFILTRATION/SUPPLY LOOP%%%%%
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * %DETERMINE f (INFILTRATION RATE)
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * test = inp(n); %TAKE OUT
                 */
                mlfAssign(
                  &test, mclArrayRef1(mclVv(inp, "inp"), mclVv(n, "n")));
                /*
                 * Fprev=F;    % store previous value of storm cumulative infiltration for accounting
                 */
                mlfAssign(&Fprev, mclVv(F, "F"));
                /*
                 * if ponding == 0 % Rainfall infiltration case
                 */
                if (mclEqBool(mclVv(ponding, "ponding"), _mxarray30_)) {
                    /*
                     * % Keep track of area under water input (infiltration) curve
                     * if n>1; curve_area = curve_area + (inp(n) + inp(n-1)) / 2 * dt; end; % trapezoidal rule
                     */
                    if (mclGtBool(mclVv(n, "n"), _mxarray12_)) {
                        mlfAssign(
                          &curve_area,
                          mclPlus(
                            mclVv(curve_area, "curve_area"),
                            mclMtimes(
                              mclMrdivide(
                                mclPlus(
                                  mclArrayRef1(
                                    mclVv(inp, "inp"), mclVv(n, "n")),
                                  mclArrayRef1(
                                    mclVv(inp, "inp"),
                                    mclMinus(mclVv(n, "n"), _mxarray12_))),
                                _mxarray13_),
                              mclVv(dt, "dt"))));
                    }
                    /*
                     * f = inp(n); F = F + f * dt;     % update storm infiltration rate & depth
                     */
                    mlfAssign(
                      &f, mclArrayRef1(mclVv(inp, "inp"), mclVv(n, "n")));
                    mlfAssign(
                      &F,
                      mclPlus(
                        mclVv(F, "F"),
                        mclMtimes(mclVv(f, "f"), mclVv(dt, "dt"))));
                    /*
                     * 
                     * 
                     * inp(n)
                     */
                    mclPrintAns(
                      &ans, mclArrayRef1(mclVv(inp, "inp"), mclVv(n, "n")));
                    /*
                     * 
                     * 
                     * 
                     * 
                     * if inp(n) > KSAT(1);              % only ponding if i>ksat
                     */
                    if (mclGtBool(
                          mclArrayRef1(mclVv(inp, "inp"), mclVv(n, "n")),
                          mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1))) {
                        /*
                         * % Compute F assuming Green-Ampt > 0
                         * FGA = KSAT(1)*abs(hwf)*(NSAT(1) - thetainiRZ)/(inp(n) - KSAT(1));
                         */
                        mlfAssign(
                          &FGA,
                          mclMrdivide(
                            mclMtimes(
                              mclMtimes(
                                mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1),
                                mlfAbs(mclVv(hwf, "hwf"))),
                              mclMinus(
                                mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1),
                                mclVv(thetainiRZ, "thetainiRZ"))),
                            mclMinus(
                              mclArrayRef1(mclVv(inp, "inp"), mclVv(n, "n")),
                              mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1))));
                        /*
                         * if curve_area > FGA
                         */
                        if (mclGtBool(
                              mclVv(curve_area, "curve_area"),
                              mclVv(FGA, "FGA"))) {
                            /*
                             * ponding = 1;
                             */
                            mlfAssign(&ponding, _mxarray12_);
                            /*
                             * tp = te;            % relative time (t is absolute)
                             */
                            mlfAssign(&tp, mclVv(te, "te"));
                            /*
                             * Fp = curve_area;    % will overestimate tp and Fp
                             */
                            mlfAssign(&Fp, mclVv(curve_area, "curve_area"));
                        /*
                         * end
                         */
                        }
                    /*
                     * end
                     */
                    }
                /*
                 * end
                 */
                }
                /*
                 * if ponding == 1 % Ponded infiltration case
                 */
                if (mclEqBool(mclVv(ponding, "ponding"), _mxarray12_)) {
                    /*
                     * A = (abs(hwf) + hp) * (NSAT(1) - thetainiRZ);        % Infilt. impetus
                     */
                    mlfAssign(
                      &A,
                      mclMtimes(
                        mclPlus(mlfAbs(mclVv(hwf, "hwf")), mclVv(hp, "hp")),
                        mclMinus(
                          mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1),
                          mclVv(thetainiRZ, "thetainiRZ"))));
                    /*
                     * [time,F1]=ode45(@iga,[t,t+dt],F,[],KSAT(1),A);    % Solve for F diff eqn
                     */
                    mlfNOde45(
                      0,
                      mlfVarargout(&time, &F1, NULL),
                      mclCreateSimpleFunctionHandle(mlxIga, "iga"),
                      mlfHorzcat(
                        mclVv(t, "t"),
                        mclPlus(mclVv(t, "t"), mclVv(dt, "dt")),
                        NULL),
                      mclVv(F, "F"),
                      _mxarray38_,
                      mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1),
                      mclVv(A, "A"),
                      NULL);
                    /*
                     * F = F1(size(F1,1));                               % Retrieve final value
                     */
                    mlfAssign(
                      &F,
                      mclArrayRef1(
                        mclVv(F1, "F1"),
                        mlfSize(
                          mclValueVarargout(), mclVv(F1, "F1"), _mxarray12_)));
                    /*
                     * f = iga(t,F,KSAT(1),A);                           % Solve for f
                     */
                    mlfAssign(
                      &f,
                      mlfIga(
                        mclVv(t, "t"),
                        mclVv(F, "F"),
                        mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1),
                        mclVv(A, "A")));
                /*
                 * end
                 */
                }
                /*
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * %Check for limits on infiltration
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * %The infil. rate cannot be > input rate plus ponding.  Therefore impose limit below.
                 * f = min(inp(n)+hp/dt,f);
                 */
                mlfAssign(
                  &f,
                  mlfMin(
                    NULL,
                    mclPlus(
                      mclArrayRef1(mclVv(inp, "inp"), mclVv(n, "n")),
                      mclMrdivide(mclVv(hp, "hp"), mclVv(dt, "dt"))),
                    mclVv(f, "f"),
                    NULL));
                /*
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * %When any combination of the soil layers is saturated, infiltration is limited by
                 * %the smallest KSAT of the saturated layers.  The limit is assigned in the section 
                 * %that addresses soil layer saturation and taken into account here.
                 * if flim >0
                 */
                if (mclGtBool(mclVv(flim, "flim"), _mxarray30_)) {
                    /*
                     * f = min(flim,f);
                     */
                    mlfAssign(
                      &f,
                      mlfMin(NULL, mclVv(flim, "flim"), mclVv(f, "f"), NULL));
                /*
                 * end
                 */
                }
                /*
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * %DETERMINE PERCOLATION RATES
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * perc1 = KSAT(1) * relcond(theta_RZ,AA(1),N(1),NRES(1),NSAT(1)); 
                 */
                mlfAssign(
                  &perc1,
                  mclMtimes(
                    mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1),
                    mlfRelcond(
                      mclVv(theta_RZ, "theta_RZ"),
                      mclIntArrayRef1(mclVv(AA, "AA"), 1),
                      mclIntArrayRef1(mclVv(N, "N"), 1),
                      mclIntArrayRef1(mclVv(NRES, "NRES"), 1),
                      mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1))));
                /*
                 * perc2 = KSAT(2) * relcond(theta_TZ,AA(2),N(2),NRES(2),NSAT(2));
                 */
                mlfAssign(
                  &perc2,
                  mclMtimes(
                    mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                    mlfRelcond(
                      mclVv(theta_TZ, "theta_TZ"),
                      mclIntArrayRef1(mclVv(AA, "AA"), 2),
                      mclIntArrayRef1(mclVv(N, "N"), 2),
                      mclIntArrayRef1(mclVv(NRES, "NRES"), 2),
                      mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2))));
                /*
                 * perc3(n,1) = KSAT(3) * relcond(theta_CZ,AA(3),N(3),NRES(3),NSAT(3));
                 */
                mclArrayAssign2(
                  &perc3,
                  mclMtimes(
                    mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                    mlfRelcond(
                      mclVv(theta_CZ, "theta_CZ"),
                      mclIntArrayRef1(mclVv(AA, "AA"), 3),
                      mclIntArrayRef1(mclVv(N, "N"), 3),
                      mclIntArrayRef1(mclVv(NRES, "NRES"), 3),
                      mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3))),
                  mclVv(n, "n"),
                  _mxarray12_);
                /*
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * %Check for limits on percolation.
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * %The percolation rate out of a soil layer cannot be greater than the KSAT of that layer.
                 * %It does not matter whether the layer is saturated or unsaturated.
                 * if perc1 > KSAT(1)
                 */
                if (mclGtBool(
                      mclVv(perc1, "perc1"),
                      mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1))) {
                    /*
                     * perc1 = KSAT(1);
                     */
                    mlfAssign(&perc1, mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1));
                /*
                 * end
                 */
                }
                /*
                 * if perc2 > KSAT(2)
                 */
                if (mclGtBool(
                      mclVv(perc2, "perc2"),
                      mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2))) {
                    /*
                     * perc2 = KSAT(2);
                     */
                    mlfAssign(&perc2, mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2));
                /*
                 * end
                 */
                }
                /*
                 * if perc3(n,1) > KSAT(3)
                 */
                if (mclGtBool(
                      mclArrayRef2(
                        mclVv(perc3, "perc3"), mclVv(n, "n"), _mxarray12_),
                      mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3))) {
                    /*
                     * perc3(n,1) = KSAT(3);
                     */
                    mclArrayAssign2(
                      &perc3,
                      mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                      mclVv(n, "n"),
                      _mxarray12_);
                /*
                 * end
                 */
                }
                /*
                 * 
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * %When any combination of the soil layers is saturated, percolation rates are
                 * %limited by the smallest KSAT of the saturated layers.  Limits are assigned
                 * %in the section that addresses soil layer saturation and taken into account here.
                 * if perc1lim> 0
                 */
                if (mclGtBool(mclVv(perc1lim, "perc1lim"), _mxarray30_)) {
                    /*
                     * perc1 = min(perc1lim,perc1);
                     */
                    mlfAssign(
                      &perc1,
                      mlfMin(
                        NULL,
                        mclVv(perc1lim, "perc1lim"),
                        mclVv(perc1, "perc1"),
                        NULL));
                /*
                 * end
                 */
                }
                /*
                 * if perc2lim > 0
                 */
                if (mclGtBool(mclVv(perc2lim, "perc2lim"), _mxarray30_)) {
                    /*
                     * perc2 = min(perc2lim,perc2);
                     */
                    mlfAssign(
                      &perc2,
                      mlfMin(
                        NULL,
                        mclVv(perc2lim, "perc2lim"),
                        mclVv(perc2, "perc2"),
                        NULL));
                /*
                 * end
                 */
                }
                /*
                 * if perc3lim > 0
                 */
                if (mclGtBool(mclVv(perc3lim, "perc3lim"), _mxarray30_)) {
                    /*
                     * perc3(n,1) = min(perc3lim,perc3(n,1));
                     */
                    mclArrayAssign2(
                      &perc3,
                      mlfMin(
                        NULL,
                        mclVv(perc3lim, "perc3lim"),
                        mclArrayRef2(
                          mclVv(perc3, "perc3"), mclVv(n, "n"), _mxarray12_),
                        NULL),
                      mclVv(n, "n"),
                      _mxarray12_);
                /*
                 * end
                 */
                }
                /*
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * % Determine Underdrain discharge rate.  Assume underdrain
                 * % discharge occurs only when Storage zone is saturated.  Assume
                 * % underdrain limited only by outgoing pipe size (orifice
                 * % equation) not by perforation in underdrain itself.
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * if udiam > 0
                 */
                if (mclGtBool(mclVv(udiam, "udiam"), _mxarray30_)) {
                    /*
                     * if theta_TZ==NSAT(2)
                     */
                    if (mclEqBool(
                          mclVv(theta_TZ, "theta_TZ"),
                          mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2))) {
                        /*
                         * hu = (Srz*(theta_RZ-NRES(1))/(NSAT(1)-NRES(1)))/100;
                         */
                        mlfAssign(
                          &hu,
                          mclMrdivide(
                            mclMrdivide(
                              mclMtimes(
                                mclVa(Srz, "Srz"),
                                mclMinus(
                                  mclVv(theta_RZ, "theta_RZ"),
                                  mclIntArrayRef1(mclVv(NRES, "NRES"), 1))),
                              mclMinus(
                                mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1),
                                mclIntArrayRef1(mclVv(NRES, "NRES"), 1))),
                            _mxarray41_));
                        /*
                         * if hp > 0
                         */
                        if (mclGtBool(mclVv(hp, "hp"), _mxarray30_)) {
                            /*
                             * hu = (hp + Srz)/100;
                             */
                            mlfAssign(
                              &hu,
                              mclMrdivide(
                                mclPlus(mclVv(hp, "hp"), mclVa(Srz, "Srz")),
                                _mxarray41_));
                        /*
                         * end
                         */
                        }
                    /*
                     * else
                     */
                    } else {
                        /*
                         * hu = 0;
                         */
                        mlfAssign(&hu, _mxarray30_);
                    /*
                     * end
                     */
                    }
                    /*
                     * uarea = pi*(udiam/2)^2;
                     */
                    mlfAssign(
                      &uarea,
                      mclMtimes(
                        mclVv(pi, "pi"),
                        mclMpower(
                          mclMrdivide(mclVv(udiam, "udiam"), _mxarray13_),
                          _mxarray13_)));
                    /*
                     * uflowmax = (Co*uarea*(2*g*hu)^.5)/Arg*3600*100;  
                     */
                    mlfAssign(
                      &uflowmax,
                      mclMtimes(
                        mclMtimes(
                          mclMrdivide(
                            mclMtimes(
                              mclMtimes(mclVv(Co, "Co"), mclVv(uarea, "uarea")),
                              mclMpower(
                                mclMtimes(
                                  mclMtimes(_mxarray13_, mclVv(g, "g")),
                                  mclVv(hu, "hu")),
                                _mxarray56_)),
                            mclVa(Arg, "Arg")),
                          _mxarray42_),
                        _mxarray41_));
                    /*
                     * %underdrain flow [cm/hr]
                     * uflow = min(uflowmax,perc1-perc2);   
                     */
                    mlfAssign(
                      &uflow,
                      mlfMin(
                        NULL,
                        mclVv(uflowmax, "uflowmax"),
                        mclMinus(mclVv(perc1, "perc1"), mclVv(perc2, "perc2")),
                        NULL));
                    /*
                     * uflow = max(uflow,0);
                     */
                    mlfAssign(
                      &uflow,
                      mlfMax(NULL, mclVv(uflow, "uflow"), _mxarray30_, NULL));
                /*
                 * else
                 */
                } else {
                    /*
                     * uflow = 0; uflowmax = 0;
                     */
                    mlfAssign(&uflow, _mxarray30_);
                    mlfAssign(&uflowmax, _mxarray30_);
                /*
                 * end
                 */
                }
                /*
                 * 
                 * 
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * %When soil moisture decreases to residual levels, percolation equals zero.  But percolation
                 * %function makes percolation essentially equal zero.  So don't need to include special
                 * %limits.
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * % Update SOIL MOISTURE values and keep old values:
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * theta_RZold=theta_RZ; theta_TZold=theta_TZ; theta_CZold=theta_CZ;
                 */
                mlfAssign(&theta_RZold, mclVv(theta_RZ, "theta_RZ"));
                mlfAssign(&theta_TZold, mclVv(theta_TZ, "theta_TZ"));
                mlfAssign(&theta_CZold, mclVv(theta_CZ, "theta_CZ"));
                /*
                 * theta_RZ = theta_RZ + (f - perc1) * dt / thick(1); % no et when ponded
                 */
                mlfAssign(
                  &theta_RZ,
                  mclPlus(
                    mclVv(theta_RZ, "theta_RZ"),
                    mclMrdivide(
                      mclMtimes(
                        mclMinus(mclVv(f, "f"), mclVv(perc1, "perc1")),
                        mclVv(dt, "dt")),
                      mclIntArrayRef1(mclVv(thick, "thick"), 1))));
                /*
                 * theta_TZ = theta_TZ + (perc1- uflow - perc2) * dt / thick(2);
                 */
                mlfAssign(
                  &theta_TZ,
                  mclPlus(
                    mclVv(theta_TZ, "theta_TZ"),
                    mclMrdivide(
                      mclMtimes(
                        mclMinus(
                          mclMinus(
                            mclVv(perc1, "perc1"), mclVv(uflow, "uflow")),
                          mclVv(perc2, "perc2")),
                        mclVv(dt, "dt")),
                      mclIntArrayRef1(mclVv(thick, "thick"), 2))));
                /*
                 * theta_CZ = theta_CZ + (perc2 - perc3(n,1)) * dt / thick(3);
                 */
                mlfAssign(
                  &theta_CZ,
                  mclPlus(
                    mclVv(theta_CZ, "theta_CZ"),
                    mclMrdivide(
                      mclMtimes(
                        mclMinus(
                          mclVv(perc2, "perc2"),
                          mclArrayRef2(
                            mclVv(perc3, "perc3"), mclVv(n, "n"), _mxarray12_)),
                        mclVv(dt, "dt")),
                      mclIntArrayRef1(mclVv(thick, "thick"), 3))));
                /*
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * % Checks for limits on soil moisture:
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * % 1. Lower: 
                 * %The soil moisture content of the layers cannot be less than the residual soil moisture.
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * %if theta_RZ <= NRES(1)
                 * %deficit = NRES(1) - theta_RZ;
                 * %theta_RZ = NRES(1);
                 * %perc1 = perc1 - deficit*thick(1)/dt; 
                 * %theta_TZ = theta_TZold + (perc1-perc2)*dt/thick(2);
                 * %end
                 * if theta_TZ <= NRES(2)
                 */
                if (mclLeBool(
                      mclVv(theta_TZ, "theta_TZ"),
                      mclIntArrayRef1(mclVv(NRES, "NRES"), 2))) {
                    /*
                     * deficit = NRES(2) - theta_TZ;
                     */
                    mlfAssign(
                      &deficit,
                      mclMinus(
                        mclIntArrayRef1(mclVv(NRES, "NRES"), 2),
                        mclVv(theta_TZ, "theta_TZ")));
                    /*
                     * theta_TZ = NRES(2);
                     */
                    mlfAssign(
                      &theta_TZ, mclIntArrayRef1(mclVv(NRES, "NRES"), 2));
                    /*
                     * perc2 = perc2 - deficit*thick(2)/dt; 
                     */
                    mlfAssign(
                      &perc2,
                      mclMinus(
                        mclVv(perc2, "perc2"),
                        mclMrdivide(
                          mclMtimes(
                            mclVv(deficit, "deficit"),
                            mclIntArrayRef1(mclVv(thick, "thick"), 2)),
                          mclVv(dt, "dt"))));
                    /*
                     * theta_CZ = theta_CZold + (perc2-perc3(n,1))*dt/thick(3);
                     */
                    mlfAssign(
                      &theta_CZ,
                      mclPlus(
                        mclVv(theta_CZold, "theta_CZold"),
                        mclMrdivide(
                          mclMtimes(
                            mclMinus(
                              mclVv(perc2, "perc2"),
                              mclArrayRef2(
                                mclVv(perc3, "perc3"),
                                mclVv(n, "n"),
                                _mxarray12_)),
                            mclVv(dt, "dt")),
                          mclIntArrayRef1(mclVv(thick, "thick"), 3))));
                /*
                 * end
                 */
                }
                /*
                 * if theta_CZ <= NRES(3)
                 */
                if (mclLeBool(
                      mclVv(theta_CZ, "theta_CZ"),
                      mclIntArrayRef1(mclVv(NRES, "NRES"), 3))) {
                    /*
                     * deficit = NRES(3) - theta_CZ;
                     */
                    mlfAssign(
                      &deficit,
                      mclMinus(
                        mclIntArrayRef1(mclVv(NRES, "NRES"), 3),
                        mclVv(theta_CZ, "theta_CZ")));
                    /*
                     * theta_CZ = NRES(3);
                     */
                    mlfAssign(
                      &theta_CZ, mclIntArrayRef1(mclVv(NRES, "NRES"), 3));
                    /*
                     * perc3(n,1) = perc3(n,1) - deficit*thick(3)/dt;
                     */
                    mclArrayAssign2(
                      &perc3,
                      mclMinus(
                        mclArrayRef2(
                          mclVv(perc3, "perc3"), mclVv(n, "n"), _mxarray12_),
                        mclMrdivide(
                          mclMtimes(
                            mclVv(deficit, "deficit"),
                            mclIntArrayRef1(mclVv(thick, "thick"), 3)),
                          mclVv(dt, "dt"))),
                      mclVv(n, "n"),
                      _mxarray12_);
                /*
                 * end
                 */
                }
                /*
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * % 2. Upper:
                 * %When any combination of the soil layers is saturated, infiltration and percolation
                 * %are limited by the smallest KSAT of the saturated layers. If the layers become
                 * %saturated in the current time step, then this section assigns appropriate limits 
                 * %for the next time step.  
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * flim=0; 
                 */
                mlfAssign(&flim, _mxarray30_);
                /*
                 * perc1lim=0; perc2lim=0; perc3lim=0;
                 */
                mlfAssign(&perc1lim, _mxarray30_);
                mlfAssign(&perc2lim, _mxarray30_);
                mlfAssign(&perc3lim, _mxarray30_);
                /*
                 * if theta_RZ >= NSAT(1)  %only RZ saturated
                 */
                if (mclGeBool(
                      mclVv(theta_RZ, "theta_RZ"),
                      mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1))) {
                    /*
                     * leftover = theta_RZ-NSAT(1);
                     */
                    mlfAssign(
                      &leftover,
                      mclMinus(
                        mclVv(theta_RZ, "theta_RZ"),
                        mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1)));
                    /*
                     * xswater = leftover*thick(1); %[cm]
                     */
                    mlfAssign(
                      &xswater,
                      mclMtimes(
                        mclVv(leftover, "leftover"),
                        mclIntArrayRef1(mclVv(thick, "thick"), 1)));
                    /*
                     * theta_RZ = NSAT(1);
                     */
                    mlfAssign(
                      &theta_RZ, mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1));
                    /*
                     * f = f - xswater/dt; %[cm/h]
                     */
                    mlfAssign(
                      &f,
                      mclMinus(
                        mclVv(f, "f"),
                        mclMrdivide(
                          mclVv(xswater, "xswater"), mclVv(dt, "dt"))));
                    /*
                     * flim = KSAT(1); %[cm/h]
                     */
                    mlfAssign(&flim, mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1));
                    /*
                     * perc1lim = KSAT(1); %[cm/h]
                     */
                    mlfAssign(
                      &perc1lim, mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1));
                    /*
                     * if theta_TZ >= NSAT(2)  %RZ and TZ saturated
                     */
                    if (mclGeBool(
                          mclVv(theta_TZ, "theta_TZ"),
                          mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2))) {
                        /*
                         * leftover = theta_TZ-NSAT(2);
                         */
                        mlfAssign(
                          &leftover,
                          mclMinus(
                            mclVv(theta_TZ, "theta_TZ"),
                            mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2)));
                        /*
                         * xswater = leftover*thick(2); %[cm]
                         */
                        mlfAssign(
                          &xswater,
                          mclMtimes(
                            mclVv(leftover, "leftover"),
                            mclIntArrayRef1(mclVv(thick, "thick"), 2)));
                        /*
                         * theta_TZ = NSAT(2);
                         */
                        mlfAssign(
                          &theta_TZ, mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2));
                        /*
                         * %perc1 = perc1 - xswater/dt;
                         * %perc1 used to calc. theta_RZ and theta_TZ but already know those values so
                         * %not necessary to keep track of perc1.
                         * f = f-xswater/dt; %[cm/h]
                         */
                        mlfAssign(
                          &f,
                          mclMinus(
                            mclVv(f, "f"),
                            mclMrdivide(
                              mclVv(xswater, "xswater"), mclVv(dt, "dt"))));
                        /*
                         * flim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]
                         */
                        mlfAssign(
                          &flim,
                          mlfMin(
                            NULL,
                            mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1),
                            mclPlus(
                              mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                              mclVv(uflowmax, "uflowmax")),
                            NULL));
                        /*
                         * perc1lim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]
                         */
                        mlfAssign(
                          &perc1lim,
                          mlfMin(
                            NULL,
                            mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1),
                            mclPlus(
                              mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                              mclVv(uflowmax, "uflowmax")),
                            NULL));
                        /*
                         * perc2lim = KSAT(2); %[cm/h]
                         */
                        mlfAssign(
                          &perc2lim, mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2));
                        /*
                         * if theta_CZ >= NSAT(3) %RZ, TZ, and CZ saturated
                         */
                        if (mclGeBool(
                              mclVv(theta_CZ, "theta_CZ"),
                              mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3))) {
                            /*
                             * leftover = theta_CZ - NSAT(3);
                             */
                            mlfAssign(
                              &leftover,
                              mclMinus(
                                mclVv(theta_CZ, "theta_CZ"),
                                mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3)));
                            /*
                             * xswater = leftover*thick(3); %[cm]
                             */
                            mlfAssign(
                              &xswater,
                              mclMtimes(
                                mclVv(leftover, "leftover"),
                                mclIntArrayRef1(mclVv(thick, "thick"), 3)));
                            /*
                             * theta_CZ = NSAT(3);
                             */
                            mlfAssign(
                              &theta_CZ,
                              mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3));
                            /*
                             * %perc2 = perc2 - xswater/dt; [cm/h]
                             * %perc1 = perc1 - xswater/dt; [cm/h]
                             * f = f-xswater/dt; %[cm/h]
                             */
                            mlfAssign(
                              &f,
                              mclMinus(
                                mclVv(f, "f"),
                                mclMrdivide(
                                  mclVv(xswater, "xswater"), mclVv(dt, "dt"))));
                            /*
                             * flim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]
                             */
                            mlfAssign(
                              &flim,
                              mlfMin(
                                NULL,
                                mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1),
                                mclPlus(
                                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                                  mclVv(uflowmax, "uflowmax")),
                                NULL));
                            /*
                             * flim = min(KSAT(3)+uflowmax,flim); %[cm/h]
                             */
                            mlfAssign(
                              &flim,
                              mlfMin(
                                NULL,
                                mclPlus(
                                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                                  mclVv(uflowmax, "uflowmax")),
                                mclVv(flim, "flim"),
                                NULL));
                            /*
                             * perc1lim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]
                             */
                            mlfAssign(
                              &perc1lim,
                              mlfMin(
                                NULL,
                                mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1),
                                mclPlus(
                                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                                  mclVv(uflowmax, "uflowmax")),
                                NULL));
                            /*
                             * perc1lim = min(KSAT(3)+uflowmax,perc1lim); %[cm/h]
                             */
                            mlfAssign(
                              &perc1lim,
                              mlfMin(
                                NULL,
                                mclPlus(
                                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                                  mclVv(uflowmax, "uflowmax")),
                                mclVv(perc1lim, "perc1lim"),
                                NULL));
                            /*
                             * perc2lim = min(KSAT(2),KSAT(3)); %[cm/h]
                             */
                            mlfAssign(
                              &perc2lim,
                              mlfMin(
                                NULL,
                                mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                                mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                                NULL));
                            /*
                             * perc3lim = KSAT(3); %[cm/h]
                             */
                            mlfAssign(
                              &perc3lim,
                              mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3));
                        /*
                         * end
                         */
                        }
                    /*
                     * elseif theta_CZ >= NSAT(3)  %RZ and CZ saturated, TZ being filled
                     */
                    } else if (mclGeBool(
                                 mclVv(theta_CZ, "theta_CZ"),
                                 mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3))) {
                        /*
                         * leftover = theta_CZ-NSAT(3);
                         */
                        mlfAssign(
                          &leftover,
                          mclMinus(
                            mclVv(theta_CZ, "theta_CZ"),
                            mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3)));
                        /*
                         * xswater = leftover*thick(3); %[cm]
                         */
                        mlfAssign(
                          &xswater,
                          mclMtimes(
                            mclVv(leftover, "leftover"),
                            mclIntArrayRef1(mclVv(thick, "thick"), 3)));
                        /*
                         * theta_CZ = NSAT(3);
                         */
                        mlfAssign(
                          &theta_CZ, mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3));
                        /*
                         * theta_TZ = theta_TZ + xswater/thick(2);                  
                         */
                        mlfAssign(
                          &theta_TZ,
                          mclPlus(
                            mclVv(theta_TZ, "theta_TZ"),
                            mclMrdivide(
                              mclVv(xswater, "xswater"),
                              mclIntArrayRef1(mclVv(thick, "thick"), 2))));
                        /*
                         * %perc2 = perc2 - xswater/dt;  %[cm/h] no need to track perc2
                         * %flim = KSAT(1); [cm/h]  value assigned above
                         * %perc1lim = KSAT(1); [cm/h]  value assigned above
                         * perc2lim = KSAT(3); %[cm/h] %
                         */
                        mlfAssign(
                          &perc2lim, mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3));
                        /*
                         * perc3lim = KSAT(3); %[cm/h]
                         */
                        mlfAssign(
                          &perc3lim, mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3));
                        /*
                         * if theta_TZ >= NSAT(2);  %RZ, TZ, and CZ saturated
                         */
                        if (mclGeBool(
                              mclVv(theta_TZ, "theta_TZ"),
                              mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2))) {
                            /*
                             * leftover = theta_TZ-NSAT(2);
                             */
                            mlfAssign(
                              &leftover,
                              mclMinus(
                                mclVv(theta_TZ, "theta_TZ"),
                                mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2)));
                            /*
                             * xswater = leftover*thick(2); %[cm]
                             */
                            mlfAssign(
                              &xswater,
                              mclMtimes(
                                mclVv(leftover, "leftover"),
                                mclIntArrayRef1(mclVv(thick, "thick"), 2)));
                            /*
                             * theta_TZ = NSAT(2);
                             */
                            mlfAssign(
                              &theta_TZ,
                              mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2));
                            /*
                             * %perc1 = perc1 - xswater/dt; [cm/h] no need to track perc1
                             * f = f - xswater/dt; %[cm/h]
                             */
                            mlfAssign(
                              &f,
                              mclMinus(
                                mclVv(f, "f"),
                                mclMrdivide(
                                  mclVv(xswater, "xswater"), mclVv(dt, "dt"))));
                            /*
                             * flim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]
                             */
                            mlfAssign(
                              &flim,
                              mlfMin(
                                NULL,
                                mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1),
                                mclPlus(
                                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                                  mclVv(uflowmax, "uflowmax")),
                                NULL));
                            /*
                             * flim = min(KSAT(3)+uflowmax,flim); %[cm/h]
                             */
                            mlfAssign(
                              &flim,
                              mlfMin(
                                NULL,
                                mclPlus(
                                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                                  mclVv(uflowmax, "uflowmax")),
                                mclVv(flim, "flim"),
                                NULL));
                            /*
                             * perc1lim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]
                             */
                            mlfAssign(
                              &perc1lim,
                              mlfMin(
                                NULL,
                                mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1),
                                mclPlus(
                                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                                  mclVv(uflowmax, "uflowmax")),
                                NULL));
                            /*
                             * perc1lim = min(KSAT(3)+uflowmax,perc1lim); %[cm/h]
                             */
                            mlfAssign(
                              &perc1lim,
                              mlfMin(
                                NULL,
                                mclPlus(
                                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                                  mclVv(uflowmax, "uflowmax")),
                                mclVv(perc1lim, "perc1lim"),
                                NULL));
                            /*
                             * perc2lim = min(KSAT(2),KSAT(3)); %[cm/h] 
                             */
                            mlfAssign(
                              &perc2lim,
                              mlfMin(
                                NULL,
                                mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                                mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                                NULL));
                        /*
                         * end
                         */
                        }
                    /*
                     * end
                     */
                    }
                /*
                 * elseif theta_TZ >= NSAT(2)  %only TZ saturated, RZ being filled
                 */
                } else if (mclGeBool(
                             mclVv(theta_TZ, "theta_TZ"),
                             mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2))) {
                    /*
                     * leftover = theta_TZ - NSAT(2);
                     */
                    mlfAssign(
                      &leftover,
                      mclMinus(
                        mclVv(theta_TZ, "theta_TZ"),
                        mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2)));
                    /*
                     * xswater = leftover*thick(2);  %[cm]
                     */
                    mlfAssign(
                      &xswater,
                      mclMtimes(
                        mclVv(leftover, "leftover"),
                        mclIntArrayRef1(mclVv(thick, "thick"), 2)));
                    /*
                     * theta_TZ = NSAT(2);
                     */
                    mlfAssign(
                      &theta_TZ, mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2));
                    /*
                     * theta_RZ = theta_RZ + xswater/thick(1);
                     */
                    mlfAssign(
                      &theta_RZ,
                      mclPlus(
                        mclVv(theta_RZ, "theta_RZ"),
                        mclMrdivide(
                          mclVv(xswater, "xswater"),
                          mclIntArrayRef1(mclVv(thick, "thick"), 1))));
                    /*
                     * %perc1 = perc1 - xswater/dt; %[cm/h] no need to track perc1
                     * perc1lim = KSAT(2)+uflowmax; %[cm/h] 
                     */
                    mlfAssign(
                      &perc1lim,
                      mclPlus(
                        mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                        mclVv(uflowmax, "uflowmax")));
                    /*
                     * perc2lim = KSAT(2); %[cm/h]
                     */
                    mlfAssign(
                      &perc2lim, mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2));
                    /*
                     * if theta_RZ >= NSAT(1)%RZ and TZ saturated
                     */
                    if (mclGeBool(
                          mclVv(theta_RZ, "theta_RZ"),
                          mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1))) {
                        /*
                         * leftover = theta_RZ - NSAT(1);
                         */
                        mlfAssign(
                          &leftover,
                          mclMinus(
                            mclVv(theta_RZ, "theta_RZ"),
                            mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1)));
                        /*
                         * xswater = leftover*thick(1); %[cm]
                         */
                        mlfAssign(
                          &xswater,
                          mclMtimes(
                            mclVv(leftover, "leftover"),
                            mclIntArrayRef1(mclVv(thick, "thick"), 1)));
                        /*
                         * theta_RZ = NSAT(1);  
                         */
                        mlfAssign(
                          &theta_RZ, mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1));
                        /*
                         * f = f - xswater/dt; %[cm/h]
                         */
                        mlfAssign(
                          &f,
                          mclMinus(
                            mclVv(f, "f"),
                            mclMrdivide(
                              mclVv(xswater, "xswater"), mclVv(dt, "dt"))));
                        /*
                         * flim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]
                         */
                        mlfAssign(
                          &flim,
                          mlfMin(
                            NULL,
                            mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1),
                            mclPlus(
                              mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                              mclVv(uflowmax, "uflowmax")),
                            NULL));
                        /*
                         * perc1lim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h] 
                         */
                        mlfAssign(
                          &perc1lim,
                          mlfMin(
                            NULL,
                            mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1),
                            mclPlus(
                              mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                              mclVv(uflowmax, "uflowmax")),
                            NULL));
                    /*
                     * %perc2lim = KSAT(2); value assigned above
                     * elseif theta_CZ >= NSAT(3) %TZ and CZ saturated, RZ being filled
                     */
                    } else if (mclGeBool(
                                 mclVv(theta_CZ, "theta_CZ"),
                                 mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3))) {
                        /*
                         * leftover = theta_CZ - NSAT(3);
                         */
                        mlfAssign(
                          &leftover,
                          mclMinus(
                            mclVv(theta_CZ, "theta_CZ"),
                            mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3)));
                        /*
                         * xswater = leftover*thick(3); %[cm]
                         */
                        mlfAssign(
                          &xswater,
                          mclMtimes(
                            mclVv(leftover, "leftover"),
                            mclIntArrayRef1(mclVv(thick, "thick"), 3)));
                        /*
                         * theta_CZ = NSAT(3);
                         */
                        mlfAssign(
                          &theta_CZ, mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3));
                        /*
                         * theta_RZ = theta_RZ + xswater/thick(1);
                         */
                        mlfAssign(
                          &theta_RZ,
                          mclPlus(
                            mclVv(theta_RZ, "theta_RZ"),
                            mclMrdivide(
                              mclVv(xswater, "xswater"),
                              mclIntArrayRef1(mclVv(thick, "thick"), 1))));
                        /*
                         * %perc2 = perc2 - xswater/dt; [cm/h] no need to track perc2
                         * %perc1 = perc1 - xswater/dt; %[cm/h] no need to track perc1
                         * perc1lim = min(KSAT(2)+uflowmax,KSAT(3)+uflowmax); %[cm/h] 
                         */
                        mlfAssign(
                          &perc1lim,
                          mlfMin(
                            NULL,
                            mclPlus(
                              mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                              mclVv(uflowmax, "uflowmax")),
                            mclPlus(
                              mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                              mclVv(uflowmax, "uflowmax")),
                            NULL));
                        /*
                         * perc2lim = min(KSAT(2),KSAT(3)); %[cm/h]
                         */
                        mlfAssign(
                          &perc2lim,
                          mlfMin(
                            NULL,
                            mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                            mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                            NULL));
                        /*
                         * perc3lim = KSAT(3); %[cm/h]
                         */
                        mlfAssign(
                          &perc3lim, mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3));
                        /*
                         * if theta_RZ >=NSAT(1) %RZ,TZ, and CZ saturated
                         */
                        if (mclGeBool(
                              mclVv(theta_RZ, "theta_RZ"),
                              mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1))) {
                            /*
                             * leftover = theta_RZ - NSAT(1);
                             */
                            mlfAssign(
                              &leftover,
                              mclMinus(
                                mclVv(theta_RZ, "theta_RZ"),
                                mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1)));
                            /*
                             * xswater = leftover*thick(1); %[cm]
                             */
                            mlfAssign(
                              &xswater,
                              mclMtimes(
                                mclVv(leftover, "leftover"),
                                mclIntArrayRef1(mclVv(thick, "thick"), 1)));
                            /*
                             * theta_RZ=NSAT(1);  
                             */
                            mlfAssign(
                              &theta_RZ,
                              mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1));
                            /*
                             * f = f - xswater/dt; %[cm/h] 
                             */
                            mlfAssign(
                              &f,
                              mclMinus(
                                mclVv(f, "f"),
                                mclMrdivide(
                                  mclVv(xswater, "xswater"), mclVv(dt, "dt"))));
                            /*
                             * flim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]
                             */
                            mlfAssign(
                              &flim,
                              mlfMin(
                                NULL,
                                mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1),
                                mclPlus(
                                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                                  mclVv(uflowmax, "uflowmax")),
                                NULL));
                            /*
                             * flim = min(KSAT(3)+uflowmax,flim); %[cm/h]
                             */
                            mlfAssign(
                              &flim,
                              mlfMin(
                                NULL,
                                mclPlus(
                                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                                  mclVv(uflowmax, "uflowmax")),
                                mclVv(flim, "flim"),
                                NULL));
                            /*
                             * perc1lim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]  
                             */
                            mlfAssign(
                              &perc1lim,
                              mlfMin(
                                NULL,
                                mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1),
                                mclPlus(
                                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                                  mclVv(uflowmax, "uflowmax")),
                                NULL));
                            /*
                             * perc1lim = min(KSAT(3)+uflowmax,perc1lim); %[cm/h] 
                             */
                            mlfAssign(
                              &perc1lim,
                              mlfMin(
                                NULL,
                                mclPlus(
                                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                                  mclVv(uflowmax, "uflowmax")),
                                mclVv(perc1lim, "perc1lim"),
                                NULL));
                        /*
                         * %perc2lim = min(KSAT(2),KSAT(3)); % [cm/h]  value assigned above
                         * end    
                         */
                        }
                    /*
                     * end    
                     */
                    }
                /*
                 * elseif theta_CZ >= NSAT(3)  %only CZ saturated, TZ being filled
                 */
                } else if (mclGeBool(
                             mclVv(theta_CZ, "theta_CZ"),
                             mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3))) {
                    /*
                     * leftover = theta_CZ - NSAT(3);
                     */
                    mlfAssign(
                      &leftover,
                      mclMinus(
                        mclVv(theta_CZ, "theta_CZ"),
                        mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3)));
                    /*
                     * xswater = leftover*thick(3);  %[cm]
                     */
                    mlfAssign(
                      &xswater,
                      mclMtimes(
                        mclVv(leftover, "leftover"),
                        mclIntArrayRef1(mclVv(thick, "thick"), 3)));
                    /*
                     * theta_CZ = NSAT(3);
                     */
                    mlfAssign(
                      &theta_CZ, mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3));
                    /*
                     * theta_TZ = theta_TZ + xswater/thick(2);               
                     */
                    mlfAssign(
                      &theta_TZ,
                      mclPlus(
                        mclVv(theta_TZ, "theta_TZ"),
                        mclMrdivide(
                          mclVv(xswater, "xswater"),
                          mclIntArrayRef1(mclVv(thick, "thick"), 2))));
                    /*
                     * %perc2 = perc2 - xswater/dt; %[cm/h] no need to track perc2
                     * perc2lim = KSAT(3); %[cm/h]  
                     */
                    mlfAssign(
                      &perc2lim, mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3));
                    /*
                     * perc3lim = KSAT(3); %[cm/h]
                     */
                    mlfAssign(
                      &perc3lim, mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3));
                    /*
                     * if theta_TZ >=NSAT(2) %TZ and CZ saturated, RZ being filled
                     */
                    if (mclGeBool(
                          mclVv(theta_TZ, "theta_TZ"),
                          mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2))) {
                        /*
                         * leftover = theta_TZ - NSAT(2); 
                         */
                        mlfAssign(
                          &leftover,
                          mclMinus(
                            mclVv(theta_TZ, "theta_TZ"),
                            mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2)));
                        /*
                         * xswater = leftover*thick(2); %[cm]
                         */
                        mlfAssign(
                          &xswater,
                          mclMtimes(
                            mclVv(leftover, "leftover"),
                            mclIntArrayRef1(mclVv(thick, "thick"), 2)));
                        /*
                         * theta_TZ = NSAT(2);  
                         */
                        mlfAssign(
                          &theta_TZ, mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2));
                        /*
                         * theta_RZ = theta_RZ + xswater/thick(1);                  
                         */
                        mlfAssign(
                          &theta_RZ,
                          mclPlus(
                            mclVv(theta_RZ, "theta_RZ"),
                            mclMrdivide(
                              mclVv(xswater, "xswater"),
                              mclIntArrayRef1(mclVv(thick, "thick"), 1))));
                        /*
                         * %perc1 = perc1 - xswater/dt; %[cm/h] no need to track perc1
                         * perc1lim = min(KSAT(2)+uflowmax,KSAT(3)+uflowmax); %[cm/h] 
                         */
                        mlfAssign(
                          &perc1lim,
                          mlfMin(
                            NULL,
                            mclPlus(
                              mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                              mclVv(uflowmax, "uflowmax")),
                            mclPlus(
                              mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                              mclVv(uflowmax, "uflowmax")),
                            NULL));
                        /*
                         * perc2lim = min(KSAT(2),KSAT(3)); %[cm/h] 
                         */
                        mlfAssign(
                          &perc2lim,
                          mlfMin(
                            NULL,
                            mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                            mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                            NULL));
                        /*
                         * if theta_RZ >= NSAT(1) %RZ, TZ, and CZ saturated
                         */
                        if (mclGeBool(
                              mclVv(theta_RZ, "theta_RZ"),
                              mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1))) {
                            /*
                             * leftover = theta_RZ - NSAT(1);
                             */
                            mlfAssign(
                              &leftover,
                              mclMinus(
                                mclVv(theta_RZ, "theta_RZ"),
                                mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1)));
                            /*
                             * xswater = leftover*thick(1); %[cm]
                             */
                            mlfAssign(
                              &xswater,
                              mclMtimes(
                                mclVv(leftover, "leftover"),
                                mclIntArrayRef1(mclVv(thick, "thick"), 1)));
                            /*
                             * theta_RZ = NSAT(1);  
                             */
                            mlfAssign(
                              &theta_RZ,
                              mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1));
                            /*
                             * f = f - xswater/dt; %[cm/h]
                             */
                            mlfAssign(
                              &f,
                              mclMinus(
                                mclVv(f, "f"),
                                mclMrdivide(
                                  mclVv(xswater, "xswater"), mclVv(dt, "dt"))));
                            /*
                             * flim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h]
                             */
                            mlfAssign(
                              &flim,
                              mlfMin(
                                NULL,
                                mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1),
                                mclPlus(
                                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                                  mclVv(uflowmax, "uflowmax")),
                                NULL));
                            /*
                             * flim = min(KSAT(3)+uflowmax,flim); %[cm/h]
                             */
                            mlfAssign(
                              &flim,
                              mlfMin(
                                NULL,
                                mclPlus(
                                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                                  mclVv(uflowmax, "uflowmax")),
                                mclVv(flim, "flim"),
                                NULL));
                            /*
                             * perc1lim = min(KSAT(1),KSAT(2)+uflowmax); %[cm/h] 
                             */
                            mlfAssign(
                              &perc1lim,
                              mlfMin(
                                NULL,
                                mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1),
                                mclPlus(
                                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                                  mclVv(uflowmax, "uflowmax")),
                                NULL));
                            /*
                             * perc1lim = min(KSAT(3)+uflowmax,perc1lim); %[cm/h] 
                             */
                            mlfAssign(
                              &perc1lim,
                              mlfMin(
                                NULL,
                                mclPlus(
                                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                                  mclVv(uflowmax, "uflowmax")),
                                mclVv(perc1lim, "perc1lim"),
                                NULL));
                        /*
                         * %perc2lim = min(KSAT(2),KSAT(3)); [cm/h]  value assigned above
                         * end
                         */
                        }
                    /*
                     * end
                     */
                    }
                /*
                 * end
                 */
                }
                /*
                 * 
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * %SURFACE WATER MASS BALANCE CALCULATIONS:
                 * hpold=hp;   % keep previous value of surface ponding depth
                 */
                mlfAssign(&hpold, mclVv(hp, "hp"));
                /*
                 * hp = hp + dt * (inp(n) - f); % Potential (Note: no evaporation)
                 */
                mlfAssign(
                  &hp,
                  mclPlus(
                    mclVv(hp, "hp"),
                    mclMtimes(
                      mclVv(dt, "dt"),
                      mclMinus(
                        mclArrayRef1(mclVv(inp, "inp"), mclVv(n, "n")),
                        mclVv(f, "f")))));
                /*
                 * 
                 * if hp <= 0   % check for no ponding, drainage in the time step
                 */
                if (mclLeBool(mclVv(hp, "hp"), _mxarray30_)) {
                    /*
                     * hp = 0;
                     */
                    mlfAssign(&hp, _mxarray30_);
                    /*
                     * ponding=0;  % end of ponding as hp goes to zero
                     */
                    mlfAssign(&ponding, _mxarray30_);
                /*
                 * else
                 */
                } else {
                    /*
                     * ponding=1;  % because hp>0
                     */
                    mlfAssign(&ponding, _mxarray12_);
                /*
                 * end
                 */
                }
                /*
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * if hp>dmin   % means water can spill from depression storage
                 */
                if (mclGtBool(mclVv(hp, "hp"), mclVa(dmin, "dmin"))) {
                    /*
                     * roff = (hp-dmin)/dt;
                     */
                    mlfAssign(
                      &roff,
                      mclMrdivide(
                        mclMinus(mclVv(hp, "hp"), mclVa(dmin, "dmin")),
                        mclVv(dt, "dt")));
                    /*
                     * hp = hp - dt * roff;
                     */
                    mlfAssign(
                      &hp,
                      mclMinus(
                        mclVv(hp, "hp"),
                        mclMtimes(mclVv(dt, "dt"), mclVv(roff, "roff"))));
                /*
                 * else
                 */
                } else {
                    /*
                     * roff = 0;
                     */
                    mlfAssign(&roff, _mxarray30_);
                /*
                 * end
                 */
                }
                /*
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * %if hp <= 0%Always 0!!!
                 * %    ttp = ttp + (te-tp);
                 * %    %roff = 0;
                 * %end;
                 * %             fprintf('V_inp= %f\n',V_inp); %TAKE OUT
                 * %             fprintf('hp= %f\n',hp); %TAKE OUT
                 * %                           fprintf('time= %f\n',t); %TAKE OUT
                 * %            fprintf('n= %f\n',n); % TAKE OUT
                 * %                         fprintf('Theta_RZ= %f\n',theta_RZ); %[TAKE THIS OUT]
                 * %             fprintf('Theta_TZ= %f\n',theta_TZ); %[TAKE THIS OUT]
                 * %fprintf('Theta_CZ= %f\n',theta_CZ); %[TAKE THIS OUT]
                 * %             fprintf('f= %f\n',f); %TAKE OUT
                 * %             fprintf('inp(n)= %f\n',inp(n)); %TAKE OUT
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * % UPDATE CUMULATIVE WATER BUDGET TERMS:
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * V_inp = V_inp + inp(n) / 100 * A_perv * dt;
                 */
                mlfAssign(
                  &V_inp,
                  mclPlus(
                    mclVv(V_inp, "V_inp"),
                    mclMtimes(
                      mclMtimes(
                        mclMrdivide(
                          mclArrayRef1(mclVv(inp, "inp"), mclVv(n, "n")),
                          _mxarray41_),
                        mclVv(A_perv, "A_perv")),
                      mclVv(dt, "dt"))));
                /*
                 * V_roff = V_roff + roff / 100 * A_perv * dt;     % Runoff volume [m^3]
                 */
                mlfAssign(
                  &V_roff,
                  mclPlus(
                    mclVv(V_roff, "V_roff"),
                    mclMtimes(
                      mclMtimes(
                        mclMrdivide(mclVv(roff, "roff"), _mxarray41_),
                        mclVv(A_perv, "A_perv")),
                      mclVv(dt, "dt"))));
                /*
                 * V_inf = V_inf + f / 100 * A_perv * dt;          % Infiltrated volume [m^3]
                 */
                mlfAssign(
                  &V_inf,
                  mclPlus(
                    mclVv(V_inf, "V_inf"),
                    mclMtimes(
                      mclMtimes(
                        mclMrdivide(mclVv(f, "f"), _mxarray41_),
                        mclVv(A_perv, "A_perv")),
                      mclVv(dt, "dt"))));
                /*
                 * V_ud = V_ud + uflow / 100 * A_perv * dt;              % Underdrain volume [m^3]
                 */
                mlfAssign(
                  &V_ud,
                  mclPlus(
                    mclVv(V_ud, "V_ud"),
                    mclMtimes(
                      mclMtimes(
                        mclMrdivide(mclVv(uflow, "uflow"), _mxarray41_),
                        mclVv(A_perv, "A_perv")),
                      mclVv(dt, "dt"))));
                /*
                 * %V_rech = V_rech + perc3(n,1) / 100 * A_perv * dt;
                 * %             V_pond = V_pond + hp/ 100 * A_perv * dt;  %TAKE OUT
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * % Keep track of saturated soil and ponding indicators:
                 * %
                 * % Update values for average water content in root zone and transmission zone:
                 * if theta_RZold > nsfactor*NSAT(1)    % near saturated in past time step, either the beginning or intermediate
                 */
                if (mclGtBool(
                      mclVv(theta_RZold, "theta_RZold"),
                      mclMtimes(
                        mclVv(nsfactor, "nsfactor"),
                        mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1)))) {
                    /*
                     * if theta_RZ > nsfactor*NSAT(1)   % near sat. now too, so keep adding:
                     */
                    if (mclGtBool(
                          mclVv(theta_RZ, "theta_RZ"),
                          mclMtimes(
                            mclVv(nsfactor, "nsfactor"),
                            mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1)))) {
                        /*
                         * tstrznsat=tstrznsat+dt; % track total simulation near sat. time
                         */
                        mlfAssign(
                          &tstrznsat,
                          mclPlus(
                            mclVv(tstrznsat, "tstrznsat"), mclVv(dt, "dt")));
                        /*
                         * tctrznsat=tctrznsat+dt; % continuous near sat. time
                         */
                        mlfAssign(
                          &tctrznsat,
                          mclPlus(
                            mclVv(tctrznsat, "tctrznsat"), mclVv(dt, "dt")));
                    /*
                     * else      % means check if maximum value in simulation so far:
                     */
                    } else {
                        /*
                         * if tctrznsat>tctrznsatmax; tctrznsatmax=tctrznsat; 
                         */
                        if (mclGtBool(
                              mclVv(tctrznsat, "tctrznsat"),
                              mclVv(tctrznsatmax, "tctrznsatmax"))) {
                            mlfAssign(
                              &tctrznsatmax, mclVv(tctrznsat, "tctrznsat"));
                        /*
                         * end; % track max cont. near sat. time
                         */
                        }
                        /*
                         * inst=inst+1; nstimes(inst,:)=[t tctrznsat];      % keep all times stored in vector [t_end tct]
                         */
                        mlfAssign(
                          &inst, mclPlus(mclVv(inst, "inst"), _mxarray12_));
                        mclArrayAssign2(
                          &nstimes,
                          mlfHorzcat(
                            mclVv(t, "t"), mclVv(tctrznsat, "tctrznsat"), NULL),
                          mclVv(inst, "inst"),
                          mlfCreateColonIndex());
                        /*
                         * tctrznsat=0;  % reset event counter
                         */
                        mlfAssign(&tctrznsat, _mxarray30_);
                    /*
                     * end
                     */
                    }
                /*
                 * end
                 */
                }
                /*
                 * %  {DBA}
                 * if theta_RZold < nwiltfactor*rzwilt
                 */
                if (mclLtBool(
                      mclVv(theta_RZold, "theta_RZold"),
                      mclMtimes(
                        mclVv(nwiltfactor, "nwiltfactor"),
                        mclVv(rzwilt, "rzwilt")))) {
                    /*
                     * if theta_RZ < nwiltfactor*rzwilt
                     */
                    if (mclLtBool(
                          mclVv(theta_RZ, "theta_RZ"),
                          mclMtimes(
                            mclVv(nwiltfactor, "nwiltfactor"),
                            mclVv(rzwilt, "rzwilt")))) {
                        /*
                         * tstrzwilt=tstrzwilt+dt; % track total simulation at wilting point
                         */
                        mlfAssign(
                          &tstrzwilt,
                          mclPlus(
                            mclVv(tstrzwilt, "tstrzwilt"), mclVv(dt, "dt")));
                        /*
                         * tctrzwilt=tctrzwilt+dt; % continuous wilting time
                         */
                        mlfAssign(
                          &tctrzwilt,
                          mclPlus(
                            mclVv(tctrzwilt, "tctrzwilt"), mclVv(dt, "dt")));
                    /*
                     * else % check if max value in simulation
                     */
                    } else {
                        /*
                         * if tctrzwilt>tctrzwiltmax; tctrzwiltmax=tctrzwilt; 
                         */
                        if (mclGtBool(
                              mclVv(tctrzwilt, "tctrzwilt"),
                              mclVv(tctrzwiltmax, "tctrzwiltmax"))) {
                            mlfAssign(
                              &tctrzwiltmax, mclVv(tctrzwilt, "tctrzwilt"));
                        /*
                         * end; % track max wilting time
                         */
                        }
                        /*
                         * inwt=inwt+1; wtimes(inwt,:)=[t tctrzwilt];    % keep track of all wilting times
                         */
                        mlfAssign(
                          &inwt, mclPlus(mclVv(inwt, "inwt"), _mxarray12_));
                        mclArrayAssign2(
                          &wtimes,
                          mlfHorzcat(
                            mclVv(t, "t"), mclVv(tctrzwilt, "tctrzwilt"), NULL),
                          mclVv(inwt, "inwt"),
                          mlfCreateColonIndex());
                        /*
                         * tctrzwilt=0; %  reset counter
                         */
                        mlfAssign(&tctrzwilt, _mxarray30_);
                    /*
                     * end
                     */
                    }
                /*
                 * end
                 */
                }
                /*
                 * % end {DBA} edits
                 * if theta_TZold > nsfactor*NSAT(2)    % near saturated in past time step, either the beginning or intermediate
                 */
                if (mclGtBool(
                      mclVv(theta_TZold, "theta_TZold"),
                      mclMtimes(
                        mclVv(nsfactor, "nsfactor"),
                        mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2)))) {
                    /*
                     * if theta_TZ > nsfactor*NSAT(2)   % near sat. now too, so keep adding:
                     */
                    if (mclGtBool(
                          mclVv(theta_TZ, "theta_TZ"),
                          mclMtimes(
                            mclVv(nsfactor, "nsfactor"),
                            mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2)))) {
                        /*
                         * tsttznsat=tsttznsat+dt; % track total simulation near sat. time
                         */
                        mlfAssign(
                          &tsttznsat,
                          mclPlus(
                            mclVv(tsttznsat, "tsttznsat"), mclVv(dt, "dt")));
                        /*
                         * tcttznsat=tcttznsat+dt; % continuous near sat. time
                         */
                        mlfAssign(
                          &tcttznsat,
                          mclPlus(
                            mclVv(tcttznsat, "tcttznsat"), mclVv(dt, "dt")));
                    /*
                     * else      % means check if maximum value in simulation so far:
                     */
                    } else {
                        /*
                         * if tcttznsat>tcttznsatmax; tcttznsatmax=tcttznsat; 
                         */
                        if (mclGtBool(
                              mclVv(tcttznsat, "tcttznsat"),
                              mclVv(tcttznsatmax, "tcttznsatmax"))) {
                            mlfAssign(
                              &tcttznsatmax, mclVv(tcttznsat, "tcttznsat"));
                        /*
                         * end; % track max cont. near sat. time
                         */
                        }
                        /*
                         * tcttznsat=0;  % reset event counter
                         */
                        mlfAssign(&tcttznsat, _mxarray30_);
                    /*
                     * end
                     */
                    }
                /*
                 * end
                 */
                }
                /*
                 * % Update values for ponding times:
                 * if hpold >= dmin  % Track overflow events
                 */
                if (mclGeBool(mclVv(hpold, "hpold"), mclVa(dmin, "dmin"))) {
                    /*
                     * if hp >= dmin
                     */
                    if (mclGeBool(mclVv(hp, "hp"), mclVa(dmin, "dmin"))) {
                        /*
                         * tctoverflow=tctoverflow+dt;
                         */
                        mlfAssign(
                          &tctoverflow,
                          mclPlus(
                            mclVv(tctoverflow, "tctoverflow"),
                            mclVv(dt, "dt")));
                    /*
                     * else
                     */
                    } else {
                        /*
                         * if t > prevof+24; inot=inot+1; oftimes(inot,:)=[t tctoverflow]; tctoverflow=0; prevof=t;end; % keep all times stored in vector
                         */
                        if (mclGtBool(
                              mclVv(t, "t"),
                              mclPlus(mclVv(prevof, "prevof"), _mxarray4_))) {
                            mlfAssign(
                              inot, mclPlus(mclVv(*inot, "inot"), _mxarray12_));
                            mclArrayAssign2(
                              &oftimes,
                              mlfHorzcat(
                                mclVv(t, "t"),
                                mclVv(tctoverflow, "tctoverflow"),
                                NULL),
                              mclVv(*inot, "inot"),
                              mlfCreateColonIndex());
                            mlfAssign(&tctoverflow, _mxarray30_);
                            mlfAssign(&prevof, mclVv(t, "t"));
                        }
                    /*
                     * end
                     */
                    }
                /*
                 * end  % Track overflow events
                 */
                }
                /*
                 * if hpold > 0    % ponded in past time step, either the beginning or intermediate
                 */
                if (mclGtBool(mclVv(hpold, "hpold"), _mxarray30_)) {
                    /*
                     * if hp > 0   % ponded now too, so keep adding:
                     */
                    if (mclGtBool(mclVv(hp, "hp"), _mxarray30_)) {
                        /*
                         * tstponded=tstponded+dt; % track total ponded time
                         */
                        mlfAssign(
                          tstponded,
                          mclPlus(
                            mclVv(*tstponded, "tstponded"), mclVv(dt, "dt")));
                        /*
                         * tctponded=tctponded+dt; % continuous ponded time
                         */
                        mlfAssign(
                          &tctponded,
                          mclPlus(
                            mclVv(tctponded, "tctponded"), mclVv(dt, "dt")));
                    /*
                     * else      % means ponding event ended, so check if maximum value in simulation so far:
                     */
                    } else {
                        /*
                         * if tctponded>tctpmax; tctpmax=tctponded; end; % track max cont. ponded t
                         */
                        if (mclGtBool(
                              mclVv(tctponded, "tctponded"),
                              mclVv(*tctpmax, "tctpmax"))) {
                            mlfAssign(tctpmax, mclVv(tctponded, "tctponded"));
                        }
                        /*
                         * inpt=inpt+1; ptimes(inpt,:)=[t tctponded];      % keep all times stored in vector [t_end tct]
                         */
                        mlfAssign(
                          &inpt, mclPlus(mclVv(inpt, "inpt"), _mxarray12_));
                        mclArrayAssign2(
                          &ptimes,
                          mlfHorzcat(
                            mclVv(t, "t"), mclVv(tctponded, "tctponded"), NULL),
                          mclVv(inpt, "inpt"),
                          mlfCreateColonIndex());
                        /*
                         * tctponded=0;  % reset event counter
                         */
                        mlfAssign(&tctponded, _mxarray30_);
                    /*
                     * end
                     */
                    }
                /*
                 * end
                 */
                }
                /*
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * %RECORD WATER BALANCE TERMS
                 * 
                 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 * Inflow(s,1)=inp(n)/4;
                 */
                mclArrayAssign2(
                  &Inflow,
                  mclMrdivide(
                    mclArrayRef1(mclVv(inp, "inp"), mclVv(n, "n")),
                    _mxarray55_),
                  mclVv(s, "s"),
                  _mxarray12_);
                /*
                 * udrain(s,1)= uflow/4;
                 */
                mclArrayAssign2(
                  &udrain,
                  mclMrdivide(mclVv(uflow, "uflow"), _mxarray55_),
                  mclVv(s, "s"),
                  _mxarray12_);
                /*
                 * infilt(s,1)= f/4;
                 */
                mclArrayAssign2(
                  &infilt,
                  mclMrdivide(mclVv(f, "f"), _mxarray55_),
                  mclVv(s, "s"),
                  _mxarray12_);
                /*
                 * Runoff(s,1)= roff/4;
                 */
                mclArrayAssign2(
                  &Runoff,
                  mclMrdivide(mclVv(roff, "roff"), _mxarray55_),
                  mclVv(s, "s"),
                  _mxarray12_);
                /*
                 * EvapoT(s,1)= 0;
                 */
                mclArrayAssign2(
                  &EvapoT, _mxarray30_, mclVv(s, "s"), _mxarray12_);
                /*
                 * Recharge(s,1)=perc3(n,1)/4;
                 */
                mclArrayAssign2(
                  &Recharge,
                  mclMrdivide(
                    mclArrayRef2(
                      mclVv(perc3, "perc3"), mclVv(n, "n"), _mxarray12_),
                    _mxarray55_),
                  mclVv(s, "s"),
                  _mxarray12_);
                /*
                 * s=s+1;
                 */
                mlfAssign(&s, mclPlus(mclVv(s, "s"), _mxarray12_));
                /*
                 * if mod(t,1)==0
                 */
                if (mclEqBool(
                      mlfMod(mclVv(t, "t"), _mxarray12_), _mxarray30_)) {
                    /*
                     * WBSum(x,:)=[t sum(Inflow(:,1)) hp sum(infilt(:,1)) sum(Runoff(:,1)) sum(udrain(:,1)) sum(Recharge(:,1)) sum(EvapoT(:,1)) theta_RZ theta_TZ theta_CZ];
                     */
                    mclArrayAssign2(
                      &WBSum,
                      mlfHorzcat(
                        mclVv(t, "t"),
                        mlfSum(
                          mclArrayRef2(
                            mclVv(Inflow, "Inflow"),
                            mlfCreateColonIndex(),
                            _mxarray12_),
                          NULL),
                        mclVv(hp, "hp"),
                        mlfSum(
                          mclArrayRef2(
                            mclVv(infilt, "infilt"),
                            mlfCreateColonIndex(),
                            _mxarray12_),
                          NULL),
                        mlfSum(
                          mclArrayRef2(
                            mclVv(Runoff, "Runoff"),
                            mlfCreateColonIndex(),
                            _mxarray12_),
                          NULL),
                        mlfSum(
                          mclArrayRef2(
                            mclVv(udrain, "udrain"),
                            mlfCreateColonIndex(),
                            _mxarray12_),
                          NULL),
                        mlfSum(
                          mclArrayRef2(
                            mclVv(Recharge, "Recharge"),
                            mlfCreateColonIndex(),
                            _mxarray12_),
                          NULL),
                        mlfSum(
                          mclArrayRef2(
                            mclVv(EvapoT, "EvapoT"),
                            mlfCreateColonIndex(),
                            _mxarray12_),
                          NULL),
                        mclVv(theta_RZ, "theta_RZ"),
                        mclVv(theta_TZ, "theta_TZ"),
                        mclVv(theta_CZ, "theta_CZ"),
                        NULL),
                      mclVv(x, "x"),
                      mlfCreateColonIndex());
                    /*
                     * Inflow(:,1)=0;
                     */
                    mclArrayAssign2(
                      &Inflow, _mxarray30_, mlfCreateColonIndex(), _mxarray12_);
                    /*
                     * udrain(:,1)= 0;
                     */
                    mclArrayAssign2(
                      &udrain, _mxarray30_, mlfCreateColonIndex(), _mxarray12_);
                    /*
                     * infilt(:,1)= 0;
                     */
                    mclArrayAssign2(
                      &infilt, _mxarray30_, mlfCreateColonIndex(), _mxarray12_);
                    /*
                     * Runoff(:,1)= 0;
                     */
                    mclArrayAssign2(
                      &Runoff, _mxarray30_, mlfCreateColonIndex(), _mxarray12_);
                    /*
                     * EvapoT(:,1)= 0;
                     */
                    mclArrayAssign2(
                      &EvapoT, _mxarray30_, mlfCreateColonIndex(), _mxarray12_);
                    /*
                     * Recharge(:,1)=0;
                     */
                    mclArrayAssign2(
                      &Recharge,
                      _mxarray30_,
                      mlfCreateColonIndex(),
                      _mxarray12_);
                    /*
                     * x=x+1;
                     */
                    mlfAssign(&x, mclPlus(mclVv(x, "x"), _mxarray12_));
                    /*
                     * s=1;
                     */
                    mlfAssign(&s, _mxarray12_);
                /*
                 * end
                 */
                }
                /*
                 * 
                 * te = te + dt;
                 */
                mlfAssign(&te, mclPlus(mclVv(te, "te"), mclVv(dt, "dt")));
                /*
                 * Fcum=Fcum+(F-Fprev);
                 */
                mlfAssign(
                  &Fcum,
                  mclPlus(
                    mclVv(Fcum, "Fcum"),
                    mclMinus(mclVv(F, "F"), mclVv(Fprev, "Fprev"))));
                /*
                 * %qst(n,:)=[t Fcum f hp theta_RZ theta_TZ theta_CZ perc1 perc2 perc3];   % value to record
                 * t = t + dt;
                 */
                mlfAssign(&t, mclPlus(mclVv(t, "t"), mclVv(dt, "dt")));
                /*
                 * %set(handles.Progress,'String',num2str(Pcompl));
                 * n = n + 1;
                 */
                mlfAssign(&n, mclPlus(mclVv(n, "n"), _mxarray12_));
                /*
                 * if t>stime, break, end
                 */
                if (mclGtBool(mclVv(t, "t"), mclVv(stime, "stime"))) {
                    break;
                }
            /*
             * 
             * %fprintf('n= %-3.0f| ',n);fprintf('In Storm\n',t);  % {DBA} screen print, omit if want to speed up program execution
             * end % OF STORM/SUPPLY LOOP
             */
            }
        /*
         * 
         * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         * else    % FREE DRAINAGE CASE (BETWEEN STORMS):
         */
        } else {
            /*
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * f = 0;      % no infiltration during hiatus periods
             */
            mlfAssign(&f, _mxarray30_);
            /*
             * roff = 0;   % no runoff either
             */
            mlfAssign(&roff, _mxarray30_);
            /*
             * test = inp(n); %TAKE OUT
             */
            mlfAssign(&test, mclArrayRef1(mclVv(inp, "inp"), mclVv(n, "n")));
            /*
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * %DETERMINE PERCOLATION RATES
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * perc1 = KSAT(1) * relcond(theta_RZ,AA(1),N(1),NRES(1),NSAT(1)); 
             */
            mlfAssign(
              &perc1,
              mclMtimes(
                mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1),
                mlfRelcond(
                  mclVv(theta_RZ, "theta_RZ"),
                  mclIntArrayRef1(mclVv(AA, "AA"), 1),
                  mclIntArrayRef1(mclVv(N, "N"), 1),
                  mclIntArrayRef1(mclVv(NRES, "NRES"), 1),
                  mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1))));
            /*
             * perc2 = KSAT(2) * relcond(theta_TZ,AA(2),N(2),NRES(2),NSAT(2));
             */
            mlfAssign(
              &perc2,
              mclMtimes(
                mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                mlfRelcond(
                  mclVv(theta_TZ, "theta_TZ"),
                  mclIntArrayRef1(mclVv(AA, "AA"), 2),
                  mclIntArrayRef1(mclVv(N, "N"), 2),
                  mclIntArrayRef1(mclVv(NRES, "NRES"), 2),
                  mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2))));
            /*
             * perc3(n,1) = KSAT(3) * relcond(theta_CZ,AA(3),N(3),NRES(3),NSAT(3));
             */
            mclArrayAssign2(
              &perc3,
              mclMtimes(
                mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                mlfRelcond(
                  mclVv(theta_CZ, "theta_CZ"),
                  mclIntArrayRef1(mclVv(AA, "AA"), 3),
                  mclIntArrayRef1(mclVv(N, "N"), 3),
                  mclIntArrayRef1(mclVv(NRES, "NRES"), 3),
                  mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3))),
              mclVv(n, "n"),
              _mxarray12_);
            /*
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * %Check for limits on percolation
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
             * %The percolation rate out of a soil layer cannot be greater than the KSAT of that layer.
             * %It does not matter whether the layer is saturated or unsaturated.
             * if perc1 > KSAT(1)
             */
            if (mclGtBool(
                  mclVv(perc1, "perc1"),
                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1))) {
                /*
                 * perc1 = KSAT(1);
                 */
                mlfAssign(&perc1, mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1));
            /*
             * end
             */
            }
            /*
             * if perc2 > KSAT(2)
             */
            if (mclGtBool(
                  mclVv(perc2, "perc2"),
                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2))) {
                /*
                 * perc2 = KSAT(2);
                 */
                mlfAssign(&perc2, mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2));
            /*
             * end
             */
            }
            /*
             * if perc3(n,1) > KSAT(3)
             */
            if (mclGtBool(
                  mclArrayRef2(
                    mclVv(perc3, "perc3"), mclVv(n, "n"), _mxarray12_),
                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3))) {
                /*
                 * perc3(n,1) = KSAT(3);
                 */
                mclArrayAssign2(
                  &perc3,
                  mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                  mclVv(n, "n"),
                  _mxarray12_);
            /*
             * end
             */
            }
            /*
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * % Determine Underdrain discharge rate.  Assume underdrain
             * % discharge occurs only when Storage zone is saturated.  Assume
             * % underdrain limited only by outgoing pipe size (orifice
             * % equation) not by perforation in underdrain itself.
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * if udiam > 0
             */
            if (mclGtBool(mclVv(udiam, "udiam"), _mxarray30_)) {
                /*
                 * if theta_TZ==NSAT(2)
                 */
                if (mclEqBool(
                      mclVv(theta_TZ, "theta_TZ"),
                      mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2))) {
                    /*
                     * hu = (Srz*(theta_RZ-NRES(1))/(NSAT(1)-NRES(1)))/100;
                     */
                    mlfAssign(
                      &hu,
                      mclMrdivide(
                        mclMrdivide(
                          mclMtimes(
                            mclVa(Srz, "Srz"),
                            mclMinus(
                              mclVv(theta_RZ, "theta_RZ"),
                              mclIntArrayRef1(mclVv(NRES, "NRES"), 1))),
                          mclMinus(
                            mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1),
                            mclIntArrayRef1(mclVv(NRES, "NRES"), 1))),
                        _mxarray41_));
                    /*
                     * if hp > 0
                     */
                    if (mclGtBool(mclVv(hp, "hp"), _mxarray30_)) {
                        /*
                         * hu = (hp + Srz)/100;
                         */
                        mlfAssign(
                          &hu,
                          mclMrdivide(
                            mclPlus(mclVv(hp, "hp"), mclVa(Srz, "Srz")),
                            _mxarray41_));
                    /*
                     * end
                     */
                    }
                /*
                 * else
                 */
                } else {
                    /*
                     * hu = 0;
                     */
                    mlfAssign(&hu, _mxarray30_);
                /*
                 * end
                 */
                }
                /*
                 * uarea = pi*(udiam/2)^2;
                 */
                mlfAssign(
                  &uarea,
                  mclMtimes(
                    mclVv(pi, "pi"),
                    mclMpower(
                      mclMrdivide(mclVv(udiam, "udiam"), _mxarray13_),
                      _mxarray13_)));
                /*
                 * uflowmax = (Co*uarea*(2*g*hu)^.5/Arg)*3600*100;  
                 */
                mlfAssign(
                  &uflowmax,
                  mclMtimes(
                    mclMtimes(
                      mclMrdivide(
                        mclMtimes(
                          mclMtimes(mclVv(Co, "Co"), mclVv(uarea, "uarea")),
                          mclMpower(
                            mclMtimes(
                              mclMtimes(_mxarray13_, mclVv(g, "g")),
                              mclVv(hu, "hu")),
                            _mxarray56_)),
                        mclVa(Arg, "Arg")),
                      _mxarray42_),
                    _mxarray41_));
            /*
             * %underdrain flow [cm/hr]
             * else
             */
            } else {
                /*
                 * uflowmax = 0;
                 */
                mlfAssign(&uflowmax, _mxarray30_);
            /*
             * end
             */
            }
            /*
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * %When any combination of the soil layers is saturated, percolation rates are
             * %limited by the smallest KSAT of the saturated layers.  The limits are assigned
             * %in the section that addresses soil layer saturation and taken into account here.
             * if perc1lim> 0
             */
            if (mclGtBool(mclVv(perc1lim, "perc1lim"), _mxarray30_)) {
                /*
                 * perc1 = min(perc1lim,perc1);
                 */
                mlfAssign(
                  &perc1,
                  mlfMin(
                    NULL,
                    mclVv(perc1lim, "perc1lim"),
                    mclVv(perc1, "perc1"),
                    NULL));
                /*
                 * perc1 = max(perc1,0);
                 */
                mlfAssign(
                  &perc1,
                  mlfMax(NULL, mclVv(perc1, "perc1"), _mxarray30_, NULL));
            /*
             * end
             */
            }
            /*
             * if perc2lim > 0
             */
            if (mclGtBool(mclVv(perc2lim, "perc2lim"), _mxarray30_)) {
                /*
                 * perc2 = min(perc2lim,perc2);
                 */
                mlfAssign(
                  &perc2,
                  mlfMin(
                    NULL,
                    mclVv(perc2lim, "perc2lim"),
                    mclVv(perc2, "perc2"),
                    NULL));
                /*
                 * perc2 = max(perc2,0);
                 */
                mlfAssign(
                  &perc2,
                  mlfMax(NULL, mclVv(perc2, "perc2"), _mxarray30_, NULL));
            /*
             * end
             */
            }
            /*
             * if perc3lim > 0
             */
            if (mclGtBool(mclVv(perc3lim, "perc3lim"), _mxarray30_)) {
                /*
                 * perc3(n,1) = min(perc3lim,perc3(n,1));
                 */
                mclArrayAssign2(
                  &perc3,
                  mlfMin(
                    NULL,
                    mclVv(perc3lim, "perc3lim"),
                    mclArrayRef2(
                      mclVv(perc3, "perc3"), mclVv(n, "n"), _mxarray12_),
                    NULL),
                  mclVv(n, "n"),
                  _mxarray12_);
                /*
                 * perc3(n,1) = max(perc3(n,1),0);
                 */
                mclArrayAssign2(
                  &perc3,
                  mlfMax(
                    NULL,
                    mclArrayRef2(
                      mclVv(perc3, "perc3"), mclVv(n, "n"), _mxarray12_),
                    _mxarray30_,
                    NULL),
                  mclVv(n, "n"),
                  _mxarray12_);
            /*
             * end
             */
            }
            /*
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * %Limit uflowmax to min. of actual flow calculated per orifice equation
             * %and max. available from percolation
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * uflow = min(uflowmax,perc1-perc2);
             */
            mlfAssign(
              &uflow,
              mlfMin(
                NULL,
                mclVv(uflowmax, "uflowmax"),
                mclMinus(mclVv(perc1, "perc1"), mclVv(perc2, "perc2")),
                NULL));
            /*
             * uflow = max(uflow,0);
             */
            mlfAssign(
              &uflow, mlfMax(NULL, mclVv(uflow, "uflow"), _mxarray30_, NULL));
            /*
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * %When soil moisture decreases to residual levels, percolation equals zero.  But percolation
             * %function makes percolation essentially equal zero.  So don't need to include special
             * %limits.
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * %Set EVAPORATION RATE
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * if theta_RZ <= rzwilt,
             */
            if (mclLeBool(
                  mclVv(theta_RZ, "theta_RZ"), mclVv(rzwilt, "rzwilt"))) {
                /*
                 * aw = 0;
                 */
                mlfAssign(&aw, _mxarray30_);
            /*
             * elseif theta_RZ >= NRES(1),
             */
            } else if (mclGeBool(
                         mclVv(theta_RZ, "theta_RZ"),
                         mclIntArrayRef1(mclVv(NRES, "NRES"), 1))) {
                /*
                 * aw = 1;
                 */
                mlfAssign(&aw, _mxarray12_);
            /*
             * else
             */
            } else {
                /*
                 * aw = (theta_RZ - rzwilt)/(NRES(1)-rzwilt);
                 */
                mlfAssign(
                  &aw,
                  mclMrdivide(
                    mclMinus(
                      mclVv(theta_RZ, "theta_RZ"), mclVv(rzwilt, "rzwilt")),
                    mclMinus(
                      mclIntArrayRef1(mclVv(NRES, "NRES"), 1),
                      mclVv(rzwilt, "rzwilt"))));
            /*
             * end
             */
            }
            /*
             * Pevap = interp1(pdfile(:,1),pdfile(:,3),t,'*linear'); %{DBA} Pot. Evaporation
             */
            mlfAssign(
              &Pevap,
              mlfInterp1(
                mclArrayRef2(
                  mclVv(pdfile, "pdfile"), mlfCreateColonIndex(), _mxarray12_),
                mclArrayRef2(
                  mclVv(pdfile, "pdfile"), mlfCreateColonIndex(), _mxarray36_),
                mclVv(t, "t"),
                _mxarray53_,
                NULL));
            /*
             * et = (1-(1+1.3*aw).^-b)*Pevap;  
             */
            mlfAssign(
              &et,
              mclMtimes(
                mclMinus(
                  _mxarray12_,
                  mlfPower(
                    mclPlus(
                      _mxarray12_, mclMtimes(_mxarray57_, mclVv(aw, "aw"))),
                    mclUminus(mclVv(b, "b")))),
                mclVv(Pevap, "Pevap")));
            /*
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * %Update SOIL MOISTURE values
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * theta_RZold=theta_RZ; theta_TZold=theta_TZ; 
             */
            mlfAssign(&theta_RZold, mclVv(theta_RZ, "theta_RZ"));
            mlfAssign(&theta_TZold, mclVv(theta_TZ, "theta_TZ"));
            /*
             * theta_CZold=theta_CZ;
             */
            mlfAssign(&theta_CZold, mclVv(theta_CZ, "theta_CZ"));
            /*
             * theta_RZ = theta_RZ + (-perc1 - et)*dt/thick(1);
             */
            mlfAssign(
              &theta_RZ,
              mclPlus(
                mclVv(theta_RZ, "theta_RZ"),
                mclMrdivide(
                  mclMtimes(
                    mclMinus(mclUminus(mclVv(perc1, "perc1")), mclVv(et, "et")),
                    mclVv(dt, "dt")),
                  mclIntArrayRef1(mclVv(thick, "thick"), 1))));
            /*
             * theta_TZ = theta_TZ + (perc1 - uflow - perc2) * dt / thick(2);
             */
            mlfAssign(
              &theta_TZ,
              mclPlus(
                mclVv(theta_TZ, "theta_TZ"),
                mclMrdivide(
                  mclMtimes(
                    mclMinus(
                      mclMinus(mclVv(perc1, "perc1"), mclVv(uflow, "uflow")),
                      mclVv(perc2, "perc2")),
                    mclVv(dt, "dt")),
                  mclIntArrayRef1(mclVv(thick, "thick"), 2))));
            /*
             * theta_CZ = theta_CZ + (perc2 - perc3(n,1)) * dt / thick(3);
             */
            mlfAssign(
              &theta_CZ,
              mclPlus(
                mclVv(theta_CZ, "theta_CZ"),
                mclMrdivide(
                  mclMtimes(
                    mclMinus(
                      mclVv(perc2, "perc2"),
                      mclArrayRef2(
                        mclVv(perc3, "perc3"), mclVv(n, "n"), _mxarray12_)),
                    mclVv(dt, "dt")),
                  mclIntArrayRef1(mclVv(thick, "thick"), 3))));
            /*
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * % Checks for limits on soil moisture:
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * % 2 Lower limits:
             * %The soil moisture content of the layers cannot be less than the residual soil moisture.
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * %if theta_RZ <= NRES(1)
             * %deficit = NRES(1) - theta_RZ;
             * %theta_RZ = NRES(1);
             * %perc1 = perc1 - deficit*thick(1)/dt; 
             * %perc1 = max(perc1,0);
             * %theta_TZ = theta_TZold + (perc1-perc2)*dt/thick(2);
             * %end
             * if theta_TZ <= NRES(2)
             */
            if (mclLeBool(
                  mclVv(theta_TZ, "theta_TZ"),
                  mclIntArrayRef1(mclVv(NRES, "NRES"), 2))) {
                /*
                 * deficit = NRES(2) - theta_TZ;
                 */
                mlfAssign(
                  &deficit,
                  mclMinus(
                    mclIntArrayRef1(mclVv(NRES, "NRES"), 2),
                    mclVv(theta_TZ, "theta_TZ")));
                /*
                 * theta_TZ = NRES(2);
                 */
                mlfAssign(&theta_TZ, mclIntArrayRef1(mclVv(NRES, "NRES"), 2));
                /*
                 * perc2 = perc2 - deficit*thick(2)/dt;
                 */
                mlfAssign(
                  &perc2,
                  mclMinus(
                    mclVv(perc2, "perc2"),
                    mclMrdivide(
                      mclMtimes(
                        mclVv(deficit, "deficit"),
                        mclIntArrayRef1(mclVv(thick, "thick"), 2)),
                      mclVv(dt, "dt"))));
                /*
                 * perc2 = max(perc2,0);
                 */
                mlfAssign(
                  &perc2,
                  mlfMax(NULL, mclVv(perc2, "perc2"), _mxarray30_, NULL));
                /*
                 * theta_CZ = theta_CZold + (perc2-perc3(n,1))*dt/thick(3);
                 */
                mlfAssign(
                  &theta_CZ,
                  mclPlus(
                    mclVv(theta_CZold, "theta_CZold"),
                    mclMrdivide(
                      mclMtimes(
                        mclMinus(
                          mclVv(perc2, "perc2"),
                          mclArrayRef2(
                            mclVv(perc3, "perc3"), mclVv(n, "n"), _mxarray12_)),
                        mclVv(dt, "dt")),
                      mclIntArrayRef1(mclVv(thick, "thick"), 3))));
            /*
             * end
             */
            }
            /*
             * if theta_CZ <= NRES(3)
             */
            if (mclLeBool(
                  mclVv(theta_CZ, "theta_CZ"),
                  mclIntArrayRef1(mclVv(NRES, "NRES"), 3))) {
                /*
                 * deficit = NRES(3) - theta_CZ;
                 */
                mlfAssign(
                  &deficit,
                  mclMinus(
                    mclIntArrayRef1(mclVv(NRES, "NRES"), 3),
                    mclVv(theta_CZ, "theta_CZ")));
                /*
                 * theta_CZ = NRES(3);
                 */
                mlfAssign(&theta_CZ, mclIntArrayRef1(mclVv(NRES, "NRES"), 3));
                /*
                 * perc3(n,1) = perc3(n,1) - deficit*thick(3)/dt;
                 */
                mclArrayAssign2(
                  &perc3,
                  mclMinus(
                    mclArrayRef2(
                      mclVv(perc3, "perc3"), mclVv(n, "n"), _mxarray12_),
                    mclMrdivide(
                      mclMtimes(
                        mclVv(deficit, "deficit"),
                        mclIntArrayRef1(mclVv(thick, "thick"), 3)),
                      mclVv(dt, "dt"))),
                  mclVv(n, "n"),
                  _mxarray12_);
                /*
                 * perc3(n,1) = max(perc3(n,1),0);
                 */
                mclArrayAssign2(
                  &perc3,
                  mlfMax(
                    NULL,
                    mclArrayRef2(
                      mclVv(perc3, "perc3"), mclVv(n, "n"), _mxarray12_),
                    _mxarray30_,
                    NULL),
                  mclVv(n, "n"),
                  _mxarray12_);
            /*
             * end
             */
            }
            /*
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * % 2. Upper limits:
             * %When any combination of the soil layers is saturated, infiltration and percolation
             * %are limited by the smallest KSAT of the saturated layers. If the layers become
             * %saturated in the current time step, then this section assigns appropriate limits 
             * %for the next time step.  Since free drainage, and infiltration equals zero, do not
             * %consider the case of the root zone being saturated.
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * perc1lim=0; perc2lim=0; perc3lim=0;
             */
            mlfAssign(&perc1lim, _mxarray30_);
            mlfAssign(&perc2lim, _mxarray30_);
            mlfAssign(&perc3lim, _mxarray30_);
            /*
             * if theta_TZ >= NSAT(2)  % TZ saturated, RZ being filled
             */
            if (mclGeBool(
                  mclVv(theta_TZ, "theta_TZ"),
                  mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2))) {
                /*
                 * leftover = theta_TZ-NSAT(2);
                 */
                mlfAssign(
                  &leftover,
                  mclMinus(
                    mclVv(theta_TZ, "theta_TZ"),
                    mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2)));
                /*
                 * xswater = leftover*thick(2); %[cm]
                 */
                mlfAssign(
                  &xswater,
                  mclMtimes(
                    mclVv(leftover, "leftover"),
                    mclIntArrayRef1(mclVv(thick, "thick"), 2)));
                /*
                 * theta_TZ = NSAT(2);
                 */
                mlfAssign(&theta_TZ, mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2));
                /*
                 * theta_RZ = theta_RZ + xswater/thick(1);        
                 */
                mlfAssign(
                  &theta_RZ,
                  mclPlus(
                    mclVv(theta_RZ, "theta_RZ"),
                    mclMrdivide(
                      mclVv(xswater, "xswater"),
                      mclIntArrayRef1(mclVv(thick, "thick"), 1))));
                /*
                 * %perc1 = perc1 - xswater/dt; %[cm/h] no need to track perc1
                 * perc1lim = KSAT(2)+uflowmax; %[cm/h] 
                 */
                mlfAssign(
                  &perc1lim,
                  mclPlus(
                    mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                    mclVv(uflowmax, "uflowmax")));
                /*
                 * perc2lim = KSAT(2); %[cm/h]
                 */
                mlfAssign(&perc2lim, mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2));
                /*
                 * %             if theta_RZ >=NSAT(1) %RZ and TZ saturated
                 * %                 theta_RZ = NSAT(1);
                 * %                 perc1lim = min(KSAT(1),KSAT(2)); %[cm/h] 
                 * %                 %perc2lim = KSAT(2); [cm/h]  value assigned above
                 * %             elseif theta_CZ >= NSAT(3) %TZ and CZ saturated, RZ being filled
                 * if theta_CZ >= NSAT(3) %TZ and CZ saturated, RZ being filled
                 */
                if (mclGeBool(
                      mclVv(theta_CZ, "theta_CZ"),
                      mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3))) {
                    /*
                     * leftover = theta_CZ - NSAT(3);
                     */
                    mlfAssign(
                      &leftover,
                      mclMinus(
                        mclVv(theta_CZ, "theta_CZ"),
                        mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3)));
                    /*
                     * xswater = leftover*thick(3); %[cm]
                     */
                    mlfAssign(
                      &xswater,
                      mclMtimes(
                        mclVv(leftover, "leftover"),
                        mclIntArrayRef1(mclVv(thick, "thick"), 3)));
                    /*
                     * theta_CZ = NSAT(3);
                     */
                    mlfAssign(
                      &theta_CZ, mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3));
                    /*
                     * theta_RZ = theta_RZ + xswater/thick(1);
                     */
                    mlfAssign(
                      &theta_RZ,
                      mclPlus(
                        mclVv(theta_RZ, "theta_RZ"),
                        mclMrdivide(
                          mclVv(xswater, "xswater"),
                          mclIntArrayRef1(mclVv(thick, "thick"), 1))));
                    /*
                     * %perc2 = perc2 - xswater/dt; [cm/h] no need to track perc2
                     * %perc1 = perc1 - xswater/dt; %[cm/h] no need to track perc2
                     * perc1lim = min(KSAT(2)+uflowmax,KSAT(3)+uflowmax); %[cm/h] 
                     */
                    mlfAssign(
                      &perc1lim,
                      mlfMin(
                        NULL,
                        mclPlus(
                          mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                          mclVv(uflowmax, "uflowmax")),
                        mclPlus(
                          mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                          mclVv(uflowmax, "uflowmax")),
                        NULL));
                    /*
                     * perc2lim = min(KSAT(2),KSAT(3)); %[cm/h]
                     */
                    mlfAssign(
                      &perc2lim,
                      mlfMin(
                        NULL,
                        mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                        mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                        NULL));
                    /*
                     * perc3lim = KSAT(3); %[cm/h]%                 if theta_RZ >=NSAT(1)%                     theta_RZ = NSAT(1);  %RZ, TZ, CZ saturated%                     perc1lim = min(KSAT(1),KSAT(2)); %[cm/h] 
                     */
                    mlfAssign(
                      &perc3lim, mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3));
                /*
                 * %                     perc1lim = min(KSAT(3),perc1lim); %[cm/h] 
                 * %                     %perc2lim = min(KSAT(2),KSAT(3)); [cm/h]  value assigned above
                 * 
                 * end    
                 */
                }
            /*
             * else if theta_CZ >= NSAT(3)  %only CZ saturated, TZ being filled
             */
            } else {
                if (mclGeBool(
                      mclVv(theta_CZ, "theta_CZ"),
                      mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3))) {
                    /*
                     * leftover = theta_CZ - NSAT(3);
                     */
                    mlfAssign(
                      &leftover,
                      mclMinus(
                        mclVv(theta_CZ, "theta_CZ"),
                        mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3)));
                    /*
                     * xswater = leftover*thick(3); %[cm]
                     */
                    mlfAssign(
                      &xswater,
                      mclMtimes(
                        mclVv(leftover, "leftover"),
                        mclIntArrayRef1(mclVv(thick, "thick"), 3)));
                    /*
                     * theta_CZ = NSAT(3);
                     */
                    mlfAssign(
                      &theta_CZ, mclIntArrayRef1(mclVv(NSAT, "NSAT"), 3));
                    /*
                     * theta_TZ = theta_TZ + xswater/thick(2); %WHAT IF THICK(2)=0??????????????????????
                     */
                    mlfAssign(
                      &theta_TZ,
                      mclPlus(
                        mclVv(theta_TZ, "theta_TZ"),
                        mclMrdivide(
                          mclVv(xswater, "xswater"),
                          mclIntArrayRef1(mclVv(thick, "thick"), 2))));
                    /*
                     * %perc2 = perc2 - leftover*thick(3)/dt; %[cm/h] no need to track perc2
                     * perc2lim = KSAT(3); %[cm/h] 
                     */
                    mlfAssign(
                      &perc2lim, mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3));
                    /*
                     * perc3lim = KSAT(3); %[cm/h]
                     */
                    mlfAssign(
                      &perc3lim, mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3));
                    /*
                     * if theta_TZ >= NSAT(2) %TZ and CZ saturated, RZ being filled
                     */
                    if (mclGeBool(
                          mclVv(theta_TZ, "theta_TZ"),
                          mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2))) {
                        /*
                         * leftover = theta_TZ-NSAT(2);
                         */
                        mlfAssign(
                          &leftover,
                          mclMinus(
                            mclVv(theta_TZ, "theta_TZ"),
                            mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2)));
                        /*
                         * xswater = leftover*thick(2); %[cm]
                         */
                        mlfAssign(
                          &xswater,
                          mclMtimes(
                            mclVv(leftover, "leftover"),
                            mclIntArrayRef1(mclVv(thick, "thick"), 2)));
                        /*
                         * theta_TZ = NSAT(2);  
                         */
                        mlfAssign(
                          &theta_TZ, mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2));
                        /*
                         * theta_RZ = theta_RZ + xswater/thick(1);
                         */
                        mlfAssign(
                          &theta_RZ,
                          mclPlus(
                            mclVv(theta_RZ, "theta_RZ"),
                            mclMrdivide(
                              mclVv(xswater, "xswater"),
                              mclIntArrayRef1(mclVv(thick, "thick"), 1))));
                        /*
                         * %perc1 = perc1 - xswater/dt; %[cm/h] no need to track perc1
                         * perc1lim = min(KSAT(2)+uflowmax,KSAT(3)+uflowmax); %[cm/h] 
                         */
                        mlfAssign(
                          &perc1lim,
                          mlfMin(
                            NULL,
                            mclPlus(
                              mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                              mclVv(uflowmax, "uflowmax")),
                            mclPlus(
                              mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                              mclVv(uflowmax, "uflowmax")),
                            NULL));
                        /*
                         * perc2lim = min(KSAT(2),KSAT(3)); %[cm/h] 
                         */
                        mlfAssign(
                          &perc2lim,
                          mlfMin(
                            NULL,
                            mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
                            mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
                            NULL));
                    /*
                     * %perc3lim = KSAT(3); [cm/h]  value assigned above
                     * %                 if theta_RZ >= NSAT(1) %RZ, TZ, CZ saturated
                     * %                     theta_RZ = NSAT(1);
                     * %                     perc1lim = min(KSAT(1),KSAT(2)); %[cm/h] 
                     * %                     perc1lim = min(KSAT(3),perc1lim); %[cm/h] 
                     * %                     % perc2lim = min(KSAT(2),KSAT(3)); [cm/h]  value assigned above
                     * %                     %perc3lim = KSAT(3); [cm/h]  value assigned above
                     * end
                     */
                    }
                /*
                 * end
                 */
                }
            /*
             * end
             */
            }
            /*
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * % UPDATE CUMULATIVE WATER BUDGET TERMS:
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * V_inp = V_inp + inp(n) / 100 * A_perv * dt;
             */
            mlfAssign(
              &V_inp,
              mclPlus(
                mclVv(V_inp, "V_inp"),
                mclMtimes(
                  mclMtimes(
                    mclMrdivide(
                      mclArrayRef1(mclVv(inp, "inp"), mclVv(n, "n")),
                      _mxarray41_),
                    mclVv(A_perv, "A_perv")),
                  mclVv(dt, "dt"))));
            /*
             * V_roff = V_roff + roff / 100 * A_perv * dt; % Runoff volume [m^3]
             */
            mlfAssign(
              &V_roff,
              mclPlus(
                mclVv(V_roff, "V_roff"),
                mclMtimes(
                  mclMtimes(
                    mclMrdivide(mclVv(roff, "roff"), _mxarray41_),
                    mclVv(A_perv, "A_perv")),
                  mclVv(dt, "dt"))));
            /*
             * V_inf = V_inf + f / 100 * A_perv * dt; % Only exfiltration in this case [m^3]
             */
            mlfAssign(
              &V_inf,
              mclPlus(
                mclVv(V_inf, "V_inf"),
                mclMtimes(
                  mclMtimes(
                    mclMrdivide(mclVv(f, "f"), _mxarray41_),
                    mclVv(A_perv, "A_perv")),
                  mclVv(dt, "dt"))));
            /*
             * %V_rech = V_rech + perc3(n,1) / 100 * A_perv * dt;
             * V_et = V_et + et / 100 * A_perv * dt; %TAKE OUT
             */
            mlfAssign(
              &V_et,
              mclPlus(
                mclVv(V_et, "V_et"),
                mclMtimes(
                  mclMtimes(
                    mclMrdivide(mclVv(et, "et"), _mxarray41_),
                    mclVv(A_perv, "A_perv")),
                  mclVv(dt, "dt"))));
            /*
             * V_ud = V_ud + uflow / 100 * A_perv * dt;
             */
            mlfAssign(
              &V_ud,
              mclPlus(
                mclVv(V_ud, "V_ud"),
                mclMtimes(
                  mclMtimes(
                    mclMrdivide(mclVv(uflow, "uflow"), _mxarray41_),
                    mclVv(A_perv, "A_perv")),
                  mclVv(dt, "dt"))));
            /*
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * % Keep track of saturated soil and ponding indicators:
             * %
             * % Update values for average water content in root zone and transmission zone:
             * if theta_RZold > nsfactor*NSAT(1)    % near saturated in past time step, either the beginning or intermediate
             */
            if (mclGtBool(
                  mclVv(theta_RZold, "theta_RZold"),
                  mclMtimes(
                    mclVv(nsfactor, "nsfactor"),
                    mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1)))) {
                /*
                 * if theta_RZ > nsfactor*NSAT(1)   % near sat. now too, so keep adding:
                 */
                if (mclGtBool(
                      mclVv(theta_RZ, "theta_RZ"),
                      mclMtimes(
                        mclVv(nsfactor, "nsfactor"),
                        mclIntArrayRef1(mclVv(NSAT, "NSAT"), 1)))) {
                    /*
                     * tstrznsat=tstrznsat+dt; % track total simulation near sat. time
                     */
                    mlfAssign(
                      &tstrznsat,
                      mclPlus(mclVv(tstrznsat, "tstrznsat"), mclVv(dt, "dt")));
                    /*
                     * tctrznsat=tctrznsat+dt; % continuous near sat. time
                     */
                    mlfAssign(
                      &tctrznsat,
                      mclPlus(mclVv(tctrznsat, "tctrznsat"), mclVv(dt, "dt")));
                /*
                 * else      % means check if maximum value in simulation so far:
                 */
                } else {
                    /*
                     * if tctrznsat>tctrznsatmax; tctrznsatmax=tctrznsat; 
                     */
                    if (mclGtBool(
                          mclVv(tctrznsat, "tctrznsat"),
                          mclVv(tctrznsatmax, "tctrznsatmax"))) {
                        mlfAssign(&tctrznsatmax, mclVv(tctrznsat, "tctrznsat"));
                    /*
                     * end; 
                     */
                    }
                    /*
                     * % track max cont. near sat. time
                     * inst=inst+1; nstimes(inst,:)=[t tctrznsat];      % keep all times stored in vector [t_end tct]
                     */
                    mlfAssign(&inst, mclPlus(mclVv(inst, "inst"), _mxarray12_));
                    mclArrayAssign2(
                      &nstimes,
                      mlfHorzcat(
                        mclVv(t, "t"), mclVv(tctrznsat, "tctrznsat"), NULL),
                      mclVv(inst, "inst"),
                      mlfCreateColonIndex());
                    /*
                     * tctrznsat=0;  % reset event counter
                     */
                    mlfAssign(&tctrznsat, _mxarray30_);
                /*
                 * end
                 */
                }
            /*
             * end
             */
            }
            /*
             * %  {DBA}
             * if theta_RZold < nwiltfactor*rzwilt
             */
            if (mclLtBool(
                  mclVv(theta_RZold, "theta_RZold"),
                  mclMtimes(
                    mclVv(nwiltfactor, "nwiltfactor"),
                    mclVv(rzwilt, "rzwilt")))) {
                /*
                 * if theta_RZ < nwiltfactor*rzwilt
                 */
                if (mclLtBool(
                      mclVv(theta_RZ, "theta_RZ"),
                      mclMtimes(
                        mclVv(nwiltfactor, "nwiltfactor"),
                        mclVv(rzwilt, "rzwilt")))) {
                    /*
                     * tstrzwilt=tstrzwilt+dt; % track total simulation at wilting point
                     */
                    mlfAssign(
                      &tstrzwilt,
                      mclPlus(mclVv(tstrzwilt, "tstrzwilt"), mclVv(dt, "dt")));
                    /*
                     * tctrzwilt=tctrzwilt+dt; % continuous wilting time
                     */
                    mlfAssign(
                      &tctrzwilt,
                      mclPlus(mclVv(tctrzwilt, "tctrzwilt"), mclVv(dt, "dt")));
                /*
                 * else % check if max value in simulation
                 */
                } else {
                    /*
                     * if tctrzwilt>tctrzwiltmax; tctrzwiltmax=tctrzwilt; end; % track max wilting time
                     */
                    if (mclGtBool(
                          mclVv(tctrzwilt, "tctrzwilt"),
                          mclVv(tctrzwiltmax, "tctrzwiltmax"))) {
                        mlfAssign(&tctrzwiltmax, mclVv(tctrzwilt, "tctrzwilt"));
                    }
                    /*
                     * inwt=inwt+1; wtimes(inwt,:)=[t tctrzwilt];    % keep track of all wilting times
                     */
                    mlfAssign(&inwt, mclPlus(mclVv(inwt, "inwt"), _mxarray12_));
                    mclArrayAssign2(
                      &wtimes,
                      mlfHorzcat(
                        mclVv(t, "t"), mclVv(tctrzwilt, "tctrzwilt"), NULL),
                      mclVv(inwt, "inwt"),
                      mlfCreateColonIndex());
                    /*
                     * tctrzwilt=0; %  reset counter
                     */
                    mlfAssign(&tctrzwilt, _mxarray30_);
                /*
                 * end
                 */
                }
            /*
             * end
             */
            }
            /*
             * % end {DBA} edits
             * if theta_TZold > nsfactor*NSAT(2)    % near saturated in past time step, either the beginning or intermediate
             */
            if (mclGtBool(
                  mclVv(theta_TZold, "theta_TZold"),
                  mclMtimes(
                    mclVv(nsfactor, "nsfactor"),
                    mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2)))) {
                /*
                 * if theta_TZ > nsfactor*NSAT(2)   % near sat. now too, so keep adding:
                 */
                if (mclGtBool(
                      mclVv(theta_TZ, "theta_TZ"),
                      mclMtimes(
                        mclVv(nsfactor, "nsfactor"),
                        mclIntArrayRef1(mclVv(NSAT, "NSAT"), 2)))) {
                    /*
                     * tsttznsat=tsttznsat+dt; % track total simulation near sat. time
                     */
                    mlfAssign(
                      &tsttznsat,
                      mclPlus(mclVv(tsttznsat, "tsttznsat"), mclVv(dt, "dt")));
                    /*
                     * tcttznsat=tcttznsat+dt; % continuous near sat. time
                     */
                    mlfAssign(
                      &tcttznsat,
                      mclPlus(mclVv(tcttznsat, "tcttznsat"), mclVv(dt, "dt")));
                /*
                 * else      % means check if maximum value in simulation so far:
                 */
                } else {
                    /*
                     * if tcttznsat>tcttznsatmax; tcttznsatmax=tcttznsat; end; 
                     */
                    if (mclGtBool(
                          mclVv(tcttznsat, "tcttznsat"),
                          mclVv(tcttznsatmax, "tcttznsatmax"))) {
                        mlfAssign(&tcttznsatmax, mclVv(tcttznsat, "tcttznsat"));
                    }
                    /*
                     * % track max cont. near sat. time
                     * tcttznsat=0;  % reset event counter
                     */
                    mlfAssign(&tcttznsat, _mxarray30_);
                /*
                 * end
                 */
                }
            /*
             * end
             */
            }
            /*
             * %         V_et = V_et + ET_avg / 100 * A_perv * dt;
             * %         Stg = V_inf-V_rech-V_et; %TAKE OUT  %Change in water stored in the soil [m^3]
             * DSSW1=(theta_RZ - theta_init(1))*thick(1)/100*A_perv; %TAKE OUT
             */
            mlfAssign(
              &DSSW1,
              mclMtimes(
                mclMrdivide(
                  mclMtimes(
                    mclMinus(
                      mclVv(theta_RZ, "theta_RZ"),
                      mclIntArrayRef1(mclVv(theta_init, "theta_init"), 1)),
                    mclIntArrayRef1(mclVv(thick, "thick"), 1)),
                  _mxarray41_),
                mclVv(A_perv, "A_perv")));
            /*
             * DSSW2=(theta_TZ - theta_init(2))*thick(2)/100*A_perv; %TAKE OUT
             */
            mlfAssign(
              &DSSW2,
              mclMtimes(
                mclMrdivide(
                  mclMtimes(
                    mclMinus(
                      mclVv(theta_TZ, "theta_TZ"),
                      mclIntArrayRef1(mclVv(theta_init, "theta_init"), 2)),
                    mclIntArrayRef1(mclVv(thick, "thick"), 2)),
                  _mxarray41_),
                mclVv(A_perv, "A_perv")));
            /*
             * DSSW3=(theta_CZ - theta_init(3))*thick(3)/100*A_perv; %TAKE OUT
             */
            mlfAssign(
              &DSSW3,
              mclMtimes(
                mclMrdivide(
                  mclMtimes(
                    mclMinus(
                      mclVv(theta_CZ, "theta_CZ"),
                      mclIntArrayRef1(mclVv(theta_init, "theta_init"), 3)),
                    mclIntArrayRef1(mclVv(thick, "thick"), 3)),
                  _mxarray41_),
                mclVv(A_perv, "A_perv")));
            /*
             * %         DSSW=DSSW1+DSSW2+DSSW3; %TAKE OUT
             * %         if Stg < DSSW %TAKE OUT
             * %             fprintf('n= %e\n', n); %TAKE OUT
             * % %         fprintf('Stg = %e\n',Stg); %TAKE OUT
             * %             fprintf('DSSW= %e\n',DSSW); %TAKE OUT
             * %         elseif Stg > DSSW %TAKE OUT
             * %             fprintf('n= %e\n', n); %TAKE OUT
             * % %         fprintf('Stg = %e\n',Stg); %TAKE OUT
             * %             fprintf('DSSW= %e\n',DSSW); %TAKE OUT
             * %         end %TAKE OUT
             * %         fprintf(' %4.0f',n);
             * %                           fprintf('time= %f\n',t); %TAKE OUT
             * %         fprintf('  %3.3f',theta_RZ);
             * %         fprintf(' %e',V_inp-V_roff);
             * %         fprintf('  %4.2f\n',hp);
             * %         qst(n,:)=[t Fcum f hp theta_RZ theta_TZ theta_CZ perc1 perc2 perc3];   % value to record
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * %RECORD WATER BALANCE TERMS
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             * Inflow(s,1)=inp(n)/4;
             */
            mclArrayAssign2(
              &Inflow,
              mclMrdivide(
                mclArrayRef1(mclVv(inp, "inp"), mclVv(n, "n")), _mxarray55_),
              mclVv(s, "s"),
              _mxarray12_);
            /*
             * udrain(s,1)= uflow/4;
             */
            mclArrayAssign2(
              &udrain,
              mclMrdivide(mclVv(uflow, "uflow"), _mxarray55_),
              mclVv(s, "s"),
              _mxarray12_);
            /*
             * infilt(s,1)= f/4;
             */
            mclArrayAssign2(
              &infilt,
              mclMrdivide(mclVv(f, "f"), _mxarray55_),
              mclVv(s, "s"),
              _mxarray12_);
            /*
             * Runoff(s,1)= roff/4;
             */
            mclArrayAssign2(
              &Runoff,
              mclMrdivide(mclVv(roff, "roff"), _mxarray55_),
              mclVv(s, "s"),
              _mxarray12_);
            /*
             * EvapoT(s,1)= et/4;
             */
            mclArrayAssign2(
              &EvapoT,
              mclMrdivide(mclVv(et, "et"), _mxarray55_),
              mclVv(s, "s"),
              _mxarray12_);
            /*
             * Recharge(s,1)=perc3(n,1)/4;
             */
            mclArrayAssign2(
              &Recharge,
              mclMrdivide(
                mclArrayRef2(mclVv(perc3, "perc3"), mclVv(n, "n"), _mxarray12_),
                _mxarray55_),
              mclVv(s, "s"),
              _mxarray12_);
            /*
             * s=s+1;
             */
            mlfAssign(&s, mclPlus(mclVv(s, "s"), _mxarray12_));
            /*
             * if mod(t,1)==0
             */
            if (mclEqBool(mlfMod(mclVv(t, "t"), _mxarray12_), _mxarray30_)) {
                /*
                 * WBSum(x,:)=[t sum(Inflow(:,1)) hp sum(infilt(:,1)) sum(Runoff(:,1)) sum(udrain(:,1)) sum(Recharge(:,1)) sum(EvapoT(:,1)) theta_RZ theta_TZ theta_CZ];
                 */
                mclArrayAssign2(
                  &WBSum,
                  mlfHorzcat(
                    mclVv(t, "t"),
                    mlfSum(
                      mclArrayRef2(
                        mclVv(Inflow, "Inflow"),
                        mlfCreateColonIndex(),
                        _mxarray12_),
                      NULL),
                    mclVv(hp, "hp"),
                    mlfSum(
                      mclArrayRef2(
                        mclVv(infilt, "infilt"),
                        mlfCreateColonIndex(),
                        _mxarray12_),
                      NULL),
                    mlfSum(
                      mclArrayRef2(
                        mclVv(Runoff, "Runoff"),
                        mlfCreateColonIndex(),
                        _mxarray12_),
                      NULL),
                    mlfSum(
                      mclArrayRef2(
                        mclVv(udrain, "udrain"),
                        mlfCreateColonIndex(),
                        _mxarray12_),
                      NULL),
                    mlfSum(
                      mclArrayRef2(
                        mclVv(Recharge, "Recharge"),
                        mlfCreateColonIndex(),
                        _mxarray12_),
                      NULL),
                    mlfSum(
                      mclArrayRef2(
                        mclVv(EvapoT, "EvapoT"),
                        mlfCreateColonIndex(),
                        _mxarray12_),
                      NULL),
                    mclVv(theta_RZ, "theta_RZ"),
                    mclVv(theta_TZ, "theta_TZ"),
                    mclVv(theta_CZ, "theta_CZ"),
                    NULL),
                  mclVv(x, "x"),
                  mlfCreateColonIndex());
                /*
                 * Inflow(:,1)=0;
                 */
                mclArrayAssign2(
                  &Inflow, _mxarray30_, mlfCreateColonIndex(), _mxarray12_);
                /*
                 * udrain(:,1)= 0;
                 */
                mclArrayAssign2(
                  &udrain, _mxarray30_, mlfCreateColonIndex(), _mxarray12_);
                /*
                 * infilt(:,1)= 0;
                 */
                mclArrayAssign2(
                  &infilt, _mxarray30_, mlfCreateColonIndex(), _mxarray12_);
                /*
                 * Runoff(:,1)= 0;
                 */
                mclArrayAssign2(
                  &Runoff, _mxarray30_, mlfCreateColonIndex(), _mxarray12_);
                /*
                 * EvapoT(:,1)= 0;
                 */
                mclArrayAssign2(
                  &EvapoT, _mxarray30_, mlfCreateColonIndex(), _mxarray12_);
                /*
                 * Recharge(:,1)=0;
                 */
                mclArrayAssign2(
                  &Recharge, _mxarray30_, mlfCreateColonIndex(), _mxarray12_);
                /*
                 * x=x+1;
                 */
                mlfAssign(&x, mclPlus(mclVv(x, "x"), _mxarray12_));
                /*
                 * s=1;
                 */
                mlfAssign(&s, _mxarray12_);
            /*
             * end
             */
            }
            /*
             * 
             * 
             * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
             * t = t + dt;
             */
            mlfAssign(&t, mclPlus(mclVv(t, "t"), mclVv(dt, "dt")));
            /*
             * %fprintf('n= %-3.0f| ',n);fprintf('free draining\n',t);  % {DBA} screen print, omit if want to speed up program execution
             * Pcompl = 100*n/(4*stime);
             */
            mlfAssign(
              &Pcompl,
              mclMrdivide(
                mclMtimes(_mxarray41_, mclVv(n, "n")),
                mclMtimes(_mxarray55_, mclVv(stime, "stime"))));
            /*
             * %set(handles.Progress,'String',num2str(Pcompl));
             * n = n + 1;
             */
            mlfAssign(&n, mclPlus(mclVv(n, "n"), _mxarray12_));
        /*
         * end
         */
        }
    /*
     * 
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * end     % END OF SIMULATION
     */
    }
    /*
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * % CHECK FOR LAST VALUES OF NEAR SATURATION AND PONDING TIMES:
     * if tctrznsat>tctrznsatmax; tctrznsatmax=tctrznsat; end; % track max cont. near sat. time
     */
    if (mclGtBool(
          mclVv(tctrznsat, "tctrznsat"), mclVv(tctrznsatmax, "tctrznsatmax"))) {
        mlfAssign(&tctrznsatmax, mclVv(tctrznsat, "tctrznsat"));
    }
    /*
     * if sum(nstimes(:,2))<tstrznsat; inst=inst+1; nstimes(inst,:)=[t tctrznsat]; end;    % keep all times stored in vector
     */
    if (mclLtBool(
          mlfSum(
            mclArrayRef2(
              mclVv(nstimes, "nstimes"), mlfCreateColonIndex(), _mxarray13_),
            NULL),
          mclVv(tstrznsat, "tstrznsat"))) {
        mlfAssign(&inst, mclPlus(mclVv(inst, "inst"), _mxarray12_));
        mclArrayAssign2(
          &nstimes,
          mlfHorzcat(mclVv(t, "t"), mclVv(tctrznsat, "tctrznsat"), NULL),
          mclVv(inst, "inst"),
          mlfCreateColonIndex());
    }
    /*
     * if tcttznsat>tcttznsatmax; tcttznsatmax=tcttznsat; end; % track max cont. near sat. time
     */
    if (mclGtBool(
          mclVv(tcttznsat, "tcttznsat"), mclVv(tcttznsatmax, "tcttznsatmax"))) {
        mlfAssign(&tcttznsatmax, mclVv(tcttznsat, "tcttznsat"));
    }
    /*
     * if tctponded>tctpmax; tctpmax=tctponded; end; % track max cont. ponded t
     */
    if (mclGtBool(mclVv(tctponded, "tctponded"), mclVv(*tctpmax, "tctpmax"))) {
        mlfAssign(tctpmax, mclVv(tctponded, "tctponded"));
    }
    /*
     * if sum(nstimes(:,2))<tstrznsat; inst=inst+1; nstimes(inst,:)=[t tctrznsat]; end;    % keep all times stored in vector
     */
    if (mclLtBool(
          mlfSum(
            mclArrayRef2(
              mclVv(nstimes, "nstimes"), mlfCreateColonIndex(), _mxarray13_),
            NULL),
          mclVv(tstrznsat, "tstrznsat"))) {
        mlfAssign(&inst, mclPlus(mclVv(inst, "inst"), _mxarray12_));
        mclArrayAssign2(
          &nstimes,
          mlfHorzcat(mclVv(t, "t"), mclVv(tctrznsat, "tctrznsat"), NULL),
          mclVv(inst, "inst"),
          mlfCreateColonIndex());
    }
    /*
     * %CHECK MASS BALANCE:
     * DSSW1=(theta_RZ - theta_init(1))*thick(1)/100*A_perv;
     */
    mlfAssign(
      &DSSW1,
      mclMtimes(
        mclMrdivide(
          mclMtimes(
            mclMinus(
              mclVv(theta_RZ, "theta_RZ"),
              mclIntArrayRef1(mclVv(theta_init, "theta_init"), 1)),
            mclIntArrayRef1(mclVv(thick, "thick"), 1)),
          _mxarray41_),
        mclVv(A_perv, "A_perv")));
    /*
     * DSSW2=(theta_TZ - theta_init(2))*thick(2)/100*A_perv;
     */
    mlfAssign(
      &DSSW2,
      mclMtimes(
        mclMrdivide(
          mclMtimes(
            mclMinus(
              mclVv(theta_TZ, "theta_TZ"),
              mclIntArrayRef1(mclVv(theta_init, "theta_init"), 2)),
            mclIntArrayRef1(mclVv(thick, "thick"), 2)),
          _mxarray41_),
        mclVv(A_perv, "A_perv")));
    /*
     * DSSW3=(theta_CZ - theta_init(3))*thick(3)/100*A_perv;
     */
    mlfAssign(
      &DSSW3,
      mclMtimes(
        mclMrdivide(
          mclMtimes(
            mclMinus(
              mclVv(theta_CZ, "theta_CZ"),
              mclIntArrayRef1(mclVv(theta_init, "theta_init"), 3)),
            mclIntArrayRef1(mclVv(thick, "thick"), 3)),
          _mxarray41_),
        mclVv(A_perv, "A_perv")));
    /*
     * DSSW=DSSW1+DSSW2+DSSW3;
     */
    mlfAssign(
      &DSSW,
      mclPlus(
        mclPlus(mclVv(DSSW1, "DSSW1"), mclVv(DSSW2, "DSSW2")),
        mclVv(DSSW3, "DSSW3")));
    /*
     * unitSM=DSSW/A_perv;
     */
    mlfAssign(
      &unitSM, mclMrdivide(mclVv(DSSW, "DSSW"), mclVv(A_perv, "A_perv")));
    /*
     * 
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * % Evaporated water from Priestley & Taylor:
     * %if etrun==1
     * %   Qevap=Arg*tau*ep(n)/100;   % potential soil evaporation
     * %else
     * %   Qevap=0;
     * %end
     * % DO SOMETHING TO WARN IF BEGINNING WITH SUPERSATURATED CONDITIONS
     * % COULD PUT SOMETHING IN THE STORM TIME LOOP TO DETECT WHETHER THERE
     * %    IS A TIME DURING EVENT WHEN PONDING = 0, AND IF SO, TO START OVER AND 
     * %    MAKE THE STORM GAP SMALLER.
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %  SUMMARIZE OUTPUT DATA
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * dVsfce=(hp-hp_0)/100*A_perv; % sfce storage change
     */
    mlfAssign(
      &dVsfce,
      mclMtimes(
        mclMrdivide(
          mclMinus(mclVv(hp, "hp"), mclVv(hp_0, "hp_0")), _mxarray41_),
        mclVv(A_perv, "A_perv")));
    /*
     * Vinsoil=V_inp-V_roff-dVsfce;
     */
    mlfAssign(
      &Vinsoil,
      mclMinus(
        mclMinus(mclVv(V_inp, "V_inp"), mclVv(V_roff, "V_roff")),
        mclVv(dVsfce, "dVsfce")));
    /*
     * infdiff=(V_inf-Vinsoil)/Vinsoil*100;
     */
    mlfAssign(
      &infdiff,
      mclMtimes(
        mclMrdivide(
          mclMinus(mclVv(V_inf, "V_inf"), mclVv(Vinsoil, "Vinsoil")),
          mclVv(Vinsoil, "Vinsoil")),
        _mxarray41_));
    /*
     * Vr=V_inf-V_et-DSSW-V_ud;   % check on Vrecharge by mass balance
     */
    mlfAssign(
      &Vr,
      mclMinus(
        mclMinus(
          mclMinus(mclVv(V_inf, "V_inf"), mclVv(V_et, "V_et")),
          mclVv(DSSW, "DSSW")),
        mclVv(V_ud, "V_ud")));
    /*
     * V_rech = sum(perc3(:,1))/100*A_perv*dt; %{DBA}
     */
    mlfAssign(
      &V_rech,
      mclMtimes(
        mclMtimes(
          mclMrdivide(
            mlfSum(
              mclArrayRef2(
                mclVv(perc3, "perc3"), mlfCreateColonIndex(), _mxarray12_),
              NULL),
            _mxarray41_),
          mclVv(A_perv, "A_perv")),
        mclVv(dt, "dt")));
    /*
     * gwrdiff=(V_rech-Vr)/Vr*100;
     */
    mlfAssign(
      &gwrdiff,
      mclMtimes(
        mclMrdivide(
          mclMinus(mclVv(V_rech, "V_rech"), mclVv(Vr, "Vr")), mclVv(Vr, "Vr")),
        _mxarray41_));
    /*
     * gwravgdepth=0.5*(V_rech+Vr)/((1+Ri+Rp)*A_perv)*100; % [cm]
     */
    mlfAssign(
      &gwravgdepth,
      mclMtimes(
        mclMrdivide(
          mclMtimes(
            _mxarray56_, mclPlus(mclVv(V_rech, "V_rech"), mclVv(Vr, "Vr"))),
          mclMtimes(
            mclPlus(mclPlus(_mxarray12_, mclVv(Ri, "Ri")), mclVv(Rp, "Rp")),
            mclVv(A_perv, "A_perv"))),
        _mxarray41_));
    /*
     * unitgwr=0.5*(V_rech+Vr)/A_perv; %{DBA} [m3/m2] recharge per unit RGarea
     */
    mlfAssign(
      &unitgwr,
      mclMrdivide(
        mclMtimes(
          _mxarray56_, mclPlus(mclVv(V_rech, "V_rech"), mclVv(Vr, "Vr"))),
        mclVv(A_perv, "A_perv")));
    /*
     * %unitgwr=V_rech/A_perv;
     * dVsoil=V_inf-V_rech-V_et;
     */
    mlfAssign(
      &dVsoil,
      mclMinus(
        mclMinus(mclVv(V_inf, "V_inf"), mclVv(V_rech, "V_rech")),
        mclVv(V_et, "V_et")));
    /*
     * unitet=V_et/A_perv; %{DBA} [m3/m2] evaporation per unit RG area
     */
    mlfAssign(
      &unitet, mclMrdivide(mclVv(V_et, "V_et"), mclVv(A_perv, "A_perv")));
    /*
     * unitron=V_inp/A_perv; %{DBA} [m3/m2] runon per unit RG area
     */
    mlfAssign(
      &unitron, mclMrdivide(mclVv(V_inp, "V_inp"), mclVv(A_perv, "A_perv")));
    /*
     * unitroff=V_roff/A_perv; %{DBA} [m3/m2] runoff per unit RG area
     */
    mlfAssign(
      &unitroff, mclMrdivide(mclVv(V_roff, "V_roff"), mclVv(A_perv, "A_perv")));
    /*
     * unitud=V_ud/A_perv; %{DBA} [m3/m2] underdrain vol per unit RG area
     */
    mlfAssign(
      &unitud, mclMrdivide(mclVv(V_ud, "V_ud"), mclVv(A_perv, "A_perv")));
    /*
     * Drunon=(100*Arg*unitron)/(Arg+TribA);
     */
    mlfAssign(
      Drunon,
      mclMrdivide(
        mclMtimes(
          mclMtimes(_mxarray41_, mclVa(Arg, "Arg")), mclVv(unitron, "unitron")),
        mclPlus(mclVa(Arg, "Arg"), mclVa(TribA, "TribA"))));
    /*
     * Drunoff=(100*Arg*unitroff)/(Arg+TribA);
     */
    mlfAssign(
      Drunoff,
      mclMrdivide(
        mclMtimes(
          mclMtimes(_mxarray41_, mclVa(Arg, "Arg")),
          mclVv(unitroff, "unitroff")),
        mclPlus(mclVa(Arg, "Arg"), mclVa(TribA, "TribA"))));
    /*
     * Drecharge=(100*Arg*unitgwr)/(Arg+TribA);
     */
    mlfAssign(
      Drecharge,
      mclMrdivide(
        mclMtimes(
          mclMtimes(_mxarray41_, mclVa(Arg, "Arg")), mclVv(unitgwr, "unitgwr")),
        mclPlus(mclVa(Arg, "Arg"), mclVa(TribA, "TribA"))));
    /*
     * Devap=(100*Arg*unitet)/(Arg+TribA);
     */
    mlfAssign(
      Devap,
      mclMrdivide(
        mclMtimes(
          mclMtimes(_mxarray41_, mclVa(Arg, "Arg")), mclVv(unitet, "unitet")),
        mclPlus(mclVa(Arg, "Arg"), mclVa(TribA, "TribA"))));
    /*
     * Dudrain=(100*Arg*unitud)/(Arg+TribA);
     */
    mlfAssign(
      Dudrain,
      mclMrdivide(
        mclMtimes(
          mclMtimes(_mxarray41_, mclVa(Arg, "Arg")), mclVv(unitud, "unitud")),
        mclPlus(mclVa(Arg, "Arg"), mclVa(TribA, "TribA"))));
    /*
     * Dstore=(100*Arg*unitSM)/(Arg+TribA);
     */
    mlfAssign(
      Dstore,
      mclMrdivide(
        mclMtimes(
          mclMtimes(_mxarray41_, mclVa(Arg, "Arg")), mclVv(unitSM, "unitSM")),
        mclPlus(mclVa(Arg, "Arg"), mclVa(TribA, "TribA"))));
    /*
     * dVs=V_inp-V_roff-V_inf;      % check on sfce storage change by mass balance
     */
    mlfAssign(
      &dVs,
      mclMinus(
        mclMinus(mclVv(V_inp, "V_inp"), mclVv(V_roff, "V_roff")),
        mclVv(V_inf, "V_inf")));
    /*
     * ttcpu=cputime-tcpu;
     */
    mlfAssign(&ttcpu, mclMinus(mlfCputime(), mclVv(tcpu, "tcpu")));
    /*
     * 
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %  CREATE OUTPUT FILE
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * if Summary == 1;
     */
    if (mclEqBool(mclVa(Summary, "Summary"), _mxarray12_)) {
        /*
         * outfile=strcat(outpre,'.txt');
         */
        mlfAssign(
          &outfile, mlfStrcat(mclVa(outpre, "outpre"), _mxarray58_, NULL));
        /*
         * fid=fopen(outfile,'w'); fclose(fid);	% delete previous version & open
         */
        mlfAssign(
          &fid,
          mlfFopen(NULL, NULL, mclVv(outfile, "outfile"), _mxarray60_, NULL));
        mclAssignAns(&ans, mlfFclose(mclVv(fid, "fid")));
        /*
         * fid=fopen(outfile,'at');
         */
        mlfAssign(
          &fid,
          mlfFopen(NULL, NULL, mclVv(outfile, "outfile"), _mxarray62_, NULL));
        /*
         * fprintf(fid,'RECARGA 2.X\n');
         */
        mclAssignAns(
          &ans, mlfNFprintf(0, mclVv(fid, "fid"), _mxarray64_, NULL));
        /*
         * fprintf(fid,'Date:\t %s\n',datestr(clock));
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray66_,
            mlfDatestr(mlfClock(), NULL, NULL),
            NULL));
        /*
         * fprintf(fid,'Output File:\t %s\n',outfile);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray68_,
            mclVv(outfile, "outfile"),
            NULL));
        /*
         * fprintf(fid,'Number of time steps=\t %-4.0f\n',n);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(0, mclVv(fid, "fid"), _mxarray70_, mclVv(n, "n"), NULL));
        /*
         * fprintf(fid,'CPU elapsed time (s)=\t %-10.0f\n\n',ttcpu);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray72_, mclVv(ttcpu, "ttcpu"), NULL));
        /*
         * fprintf(fid,'INPUT TERMS\n');
         */
        mclAssignAns(
          &ans, mlfNFprintf(0, mclVv(fid, "fid"), _mxarray74_, NULL));
        /*
         * fprintf(fid,'%%Impervious=\t %g\n',EIA*100);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray76_,
            mclMtimes(mclVa(EIA, "EIA"), _mxarray41_),
            NULL));
        /*
         * fprintf(fid,'Pervious CN=\t %g\n',CN);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray78_, mclVa(CN, "CN"), NULL));
        /*
         * fprintf(fid,'Facility Area (m2)=\t %g\n',Arg);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray80_, mclVa(Arg, "Arg"), NULL));
        /*
         * fprintf(fid,'Trib. Area (m2)=\t %g\n',TribA);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray82_, mclVa(TribA, "TribA"), NULL));
        /*
         * fprintf(fid,'RATIOimp2rg=\t %g\n',Ri);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray84_, mclVv(Ri, "Ri"), NULL));
        /*
         * fprintf(fid,'RATIOperv2rg=\t %g\n\n',Rp);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray86_, mclVv(Rp, "Rp"), NULL));
        /*
         * fprintf(fid,'LAYER       \t DEPTH (cm)  \t Ksat (cm/hr)\t TEXTURE     \n');
         */
        mclAssignAns(
          &ans, mlfNFprintf(0, mclVv(fid, "fid"), _mxarray88_, NULL));
        /*
         * fprintf(fid,'Depression  \t     %g      \t     ---      \t   ---       \n',dmin);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray90_, mclVa(dmin, "dmin"), NULL));
        /*
         * fprintf(fid,'Root Zone   \t     %g      \t',thick(1));
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray92_,
            mclIntArrayRef1(mclVv(thick, "thick"), 1),
            NULL));
        /*
         * fprintf(fid,'    %4.2f    \t',KSAT(1));
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray94_,
            mclIntArrayRef1(mclVv(KSAT, "KSAT"), 1),
            NULL));
        /*
         * fprintf(fid,'    %g      \n',Rsoilnumber);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray96_,
            mclVa(Rsoilnumber, "Rsoilnumber"),
            NULL));
        /*
         * fprintf(fid,'Storage     \t     %g      \t',thick(2));
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray98_,
            mclIntArrayRef1(mclVv(thick, "thick"), 2),
            NULL));
        /*
         * fprintf(fid,'    %4.2f     \t',KSAT(2));
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray100_,
            mclIntArrayRef1(mclVv(KSAT, "KSAT"), 2),
            NULL));
        /*
         * fprintf(fid,'    %g      \n',Ssoilnumber);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray96_,
            mclVa(Ssoilnumber, "Ssoilnumber"),
            NULL));
        /*
         * fprintf(fid,'Native      \t     ---     \t     %4.2f   \t',KSAT(3));
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray102_,
            mclIntArrayRef1(mclVv(KSAT, "KSAT"), 3),
            NULL));
        /*
         * fprintf(fid,'    %g      \n',Nsoilnumber);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray96_,
            mclVa(Nsoilnumber, "Nsoilnumber"),
            NULL));
        /*
         * fprintf(fid,'\nInit. Cond.: suction head (cm)=\t %5.0f\n',hin);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray104_, mclVv(hin, "hin"), NULL));
        /*
         * fprintf(fid,'Max. Underdrain flow (cm/hr)=\t %5.0f\n',Umax);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray106_, mclVv(Umax, "Umax"), NULL));
        /*
         * fprintf(fid,'Underdrain diam (mm)=\t %5.0f\n',udiam*1000);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray108_,
            mclMtimes(mclVv(udiam, "udiam"), _mxarray40_),
            NULL));
        /*
         * fprintf(fid,'\nWATER BALANCE TERMS\n\n');
         */
        mclAssignAns(
          &ans, mlfNFprintf(0, mclVv(fid, "fid"), _mxarray110_, NULL));
        /*
         * fprintf(fid,'           \t Volume(m3) \t Depth (cm)\t %%of Inflow \n');
         */
        mclAssignAns(
          &ans, mlfNFprintf(0, mclVv(fid, "fid"), _mxarray112_, NULL));
        /*
         * fprintf(fid,'Inflow     \t %10.2f \t',V_inp);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray114_, mclVv(V_inp, "V_inp"), NULL));
        /*
         * fprintf(fid,' %9.2f \t',Drunon);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray116_,
            mclVv(*Drunon, "Drunon"),
            NULL));
        /*
         * fprintf(fid,' %9.2f \n',100*Drunon/Drunon);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray118_,
            mclMrdivide(
              mclMtimes(_mxarray41_, mclVv(*Drunon, "Drunon")),
              mclVv(*Drunon, "Drunon")),
            NULL));
        /*
         * fprintf(fid,'Runoff     \t %10.2f \t',V_roff);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray120_, mclVv(V_roff, "V_roff"), NULL));
        /*
         * fprintf(fid,' %9.2f \t',Drunoff);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray116_,
            mclVv(*Drunoff, "Drunoff"),
            NULL));
        /*
         * fprintf(fid,' %9.2f \n',100*Drunoff/Drunon);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray118_,
            mclMrdivide(
              mclMtimes(_mxarray41_, mclVv(*Drunoff, "Drunoff")),
              mclVv(*Drunon, "Drunon")),
            NULL));
        /*
         * fprintf(fid,'Recharge   \t %10.2f \t',V_rech);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray122_, mclVv(V_rech, "V_rech"), NULL));
        /*
         * fprintf(fid,' %9.2f \t',Drecharge);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray116_,
            mclVv(*Drecharge, "Drecharge"),
            NULL));
        /*
         * fprintf(fid,' %9.2f \n',100*Drecharge/Drunon);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray118_,
            mclMrdivide(
              mclMtimes(_mxarray41_, mclVv(*Drecharge, "Drecharge")),
              mclVv(*Drunon, "Drunon")),
            NULL));
        /*
         * fprintf(fid,'Evaporation\t %10.2f \t',V_et);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray124_, mclVv(V_et, "V_et"), NULL));
        /*
         * fprintf(fid,' %9.2f \t',Devap);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray116_, mclVv(*Devap, "Devap"), NULL));
        /*
         * fprintf(fid,' %9.2f \n',100*Devap/Drunon);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray118_,
            mclMrdivide(
              mclMtimes(_mxarray41_, mclVv(*Devap, "Devap")),
              mclVv(*Drunon, "Drunon")),
            NULL));
        /*
         * fprintf(fid,'Underdrain \t %10.2f \t',V_ud);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray126_, mclVv(V_ud, "V_ud"), NULL));
        /*
         * fprintf(fid,' %9.2f \t',Dudrain);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray116_,
            mclVv(*Dudrain, "Dudrain"),
            NULL));
        /*
         * fprintf(fid,' %9.2f \n',100*Dudrain/Drunon);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray118_,
            mclMrdivide(
              mclMtimes(_mxarray41_, mclVv(*Dudrain, "Dudrain")),
              mclVv(*Drunon, "Drunon")),
            NULL));
        /*
         * fprintf(fid,'Storage    \t %10.2f \t',DSSW);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray128_, mclVv(DSSW, "DSSW"), NULL));
        /*
         * fprintf(fid,' %9.2f \t',Dstore);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray116_,
            mclVv(*Dstore, "Dstore"),
            NULL));
        /*
         * fprintf(fid,' %9.2f \n',100*Dstore/Drunon);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray118_,
            mclMrdivide(
              mclMtimes(_mxarray41_, mclVv(*Dstore, "Dstore")),
              mclVv(*Drunon, "Drunon")),
            NULL));
        /*
         * fprintf(fid,'Infiltrated\t %10.2f \n',V_inf);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray130_, mclVv(V_inf, "V_inf"), NULL));
        /*
         * fprintf(fid,'Precip.    \t            \t %9.2f\n',Dprecip);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray132_,
            mclVv(Dprecip, "Dprecip"),
            NULL));
        /*
         * fprintf(fid,'Imp. Inflow\t %10.2f \t',Vimp);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray134_, mclVv(Vimp, "Vimp"), NULL));
        /*
         * fprintf(fid,' %9.2f \t',Dimp);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray116_, mclVv(*Dimp, "Dimp"), NULL));
        /*
         * fprintf(fid,' %9.2f \n',100*Vimp/V_inp);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray118_,
            mclMrdivide(
              mclMtimes(_mxarray41_, mclVv(Vimp, "Vimp")),
              mclVv(V_inp, "V_inp")),
            NULL));
        /*
         * fprintf(fid,'Perv.Inflow\t %10.2f \t',Vperv);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray136_, mclVv(Vperv, "Vperv"), NULL));
        /*
         * fprintf(fid,' %9.2f \t',Dperv);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray116_, mclVv(*Dperv, "Dperv"), NULL));
        /*
         * fprintf(fid,' %9.2f \n',100*Vperv/V_inp);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray118_,
            mclMrdivide(
              mclMtimes(_mxarray41_, mclVv(Vperv, "Vperv")),
              mclVv(V_inp, "V_inp")),
            NULL));
        /*
         * fprintf(fid,'Basin Stayon(cm)\t\t %9.2f\n\n',Dbasin);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray138_, mclVv(Dbasin, "Dbasin"), NULL));
        /*
         * fprintf(fid,'Mass Balance Checks:\n');
         */
        mclAssignAns(
          &ans, mlfNFprintf(0, mclVv(fid, "fid"), _mxarray140_, NULL));
        /*
         * fprintf(fid,'Vinsoil=input-roff-dvh\n');
         */
        mclAssignAns(
          &ans, mlfNFprintf(0, mclVv(fid, "fid"), _mxarray142_, NULL));
        /*
         * fprintf(fid,'Vinsoil (m3)=\t %e\n',Vinsoil);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray144_,
            mclVv(Vinsoil, "Vinsoil"),
            NULL));
        /*
         * fprintf(fid,'Vinf difference (%%)=\t %4.1f\n',infdiff);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray146_,
            mclVv(infdiff, "infdiff"),
            NULL));
        /*
         * fprintf(fid,'Vrecharg1=inf-et-dsws\n');
         */
        mclAssignAns(
          &ans, mlfNFprintf(0, mclVv(fid, "fid"), _mxarray148_, NULL));
        /*
         * fprintf(fid,'Vrecharg1 (m3)=\t %e\n',Vr);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray150_, mclVv(Vr, "Vr"), NULL));
        /*
         * fprintf(fid,'Vrech difference (%%)=\t %4.1f\n',gwrdiff);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray152_,
            mclVv(gwrdiff, "gwrdiff"),
            NULL));
        /*
         * fprintf(fid,'dVsoil=inf-rech-et\n');
         */
        mclAssignAns(
          &ans, mlfNFprintf(0, mclVv(fid, "fid"), _mxarray154_, NULL));
        /*
         * fprintf(fid,'dVsoil=\t %e\n',dVsoil);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray156_, mclVv(dVsoil, "dVsoil"), NULL));
        /*
         * fprintf(fid,'dVsfce (m3)=\t %e\n',dVsfce);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray158_, mclVv(dVsfce, "dVsfce"), NULL));
        /*
         * fprintf(fid,'dVs2=input-roff-inf\n');
         */
        mclAssignAns(
          &ans, mlfNFprintf(0, mclVv(fid, "fid"), _mxarray160_, NULL));
        /*
         * fprintf(fid,'dVs2 (m3)= %e\n',dVs);%
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray162_, mclVv(dVs, "dVs"), NULL));
        /*
         * fprintf(fid,'\nPLANT SURVIVABILITY TERMS\n');
         */
        mclAssignAns(
          &ans, mlfNFprintf(0, mclVv(fid, "fid"), _mxarray164_, NULL));
        /*
         * fprintf(fid,'total time ponded (h)=\t %5.1f\n',tstponded);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray166_,
            mclVv(*tstponded, "tstponded"),
            NULL));
        /*
         * fprintf(fid,'max time ponded (h)=\t %5.1f\n',tctpmax);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray168_,
            mclVv(*tctpmax, "tctpmax"),
            NULL));
        /*
         * fprintf(fid,'total time RZ saturated (h)=\t %5.1f\n',tstrznsat);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray170_,
            mclVv(tstrznsat, "tstrznsat"),
            NULL));
        /*
         * fprintf(fid,'total time RZ at wilting point (h)=\t %5.1f\n', tstrzwilt);  % {DBA}
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray172_,
            mclVv(tstrzwilt, "tstrzwilt"),
            NULL));
        /*
         * fprintf(fid,'total time TZ saturated (h)=\t %5.1f\n',tsttznsat);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray174_,
            mclVv(tsttznsat, "tsttznsat"),
            NULL));
        /*
         * fprintf(fid,'max time RZ saturated (h)=\t %5.1f\n',tctrznsatmax);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray176_,
            mclVv(tctrznsatmax, "tctrznsatmax"),
            NULL));
        /*
         * fprintf(fid,'max time RZ at wilting point (h)=\t %5.1f\n',tctrzwiltmax);  %{DBA}
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray178_,
            mclVv(tctrzwiltmax, "tctrzwiltmax"),
            NULL));
        /*
         * fprintf(fid,'max time TZ saturated (h)=\t %5.1f\n',tcttznsatmax);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray180_,
            mclVv(tcttznsatmax, "tcttznsatmax"),
            NULL));
        /*
         * fprintf (fid,'near saturated times for simulation=\t %6.0f\n', inst); fprintf(fid,'%6.1f\t %6.1f\n',nstimes');
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray182_, mclVv(inst, "inst"), NULL));
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray184_,
            mlfCtranspose(mclVv(nstimes, "nstimes")),
            NULL));
        /*
         * fprintf (fid,'wilting point times for simulation=\t %6.0f\n',inwt); fprintf(fid,'%6.1f\t %6.1f\n',wtimes'); % {DBA}
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray186_, mclVv(inwt, "inwt"), NULL));
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray184_,
            mlfCtranspose(mclVv(wtimes, "wtimes")),
            NULL));
        /*
         * fprintf (fid,'ponded times for simulation=\t %6.0f\n', inpt); fprintf(fid,'%6.1f\t %6.1f\n',ptimes');
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray188_, mclVv(inpt, "inpt"), NULL));
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray184_,
            mlfCtranspose(mclVv(ptimes, "ptimes")),
            NULL));
        /*
         * fprintf (fid,'overflow times for simulation=\t %6.0f\n',inot); fprintf(fid,'%6.1f\t %6.1f\n',oftimes');
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0, mclVv(fid, "fid"), _mxarray190_, mclVv(*inot, "inot"), NULL));
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray184_,
            mlfCtranspose(mclVv(oftimes, "oftimes")),
            NULL));
        /*
         * fclose(fid);
         */
        mclAssignAns(&ans, mlfFclose(mclVv(fid, "fid")));
    /*
     * end
     */
    }
    /*
     * if Record == 1;
     */
    if (mclEqBool(mclVa(Record, "Record"), _mxarray12_)) {
        /*
         * WBfile = strcat(outpre,'_WB.txt');    
         */
        mlfAssign(
          &WBfile, mlfStrcat(mclVa(outpre, "outpre"), _mxarray192_, NULL));
        /*
         * fid=fopen(WBfile,'w'); fclose(fid);	% delete previous version & open
         */
        mlfAssign(
          &fid,
          mlfFopen(NULL, NULL, mclVv(WBfile, "WBfile"), _mxarray60_, NULL));
        mclAssignAns(&ans, mlfFclose(mclVv(fid, "fid")));
        /*
         * fid=fopen(WBfile,'at');
         */
        mlfAssign(
          &fid,
          mlfFopen(NULL, NULL, mclVv(WBfile, "WBfile"), _mxarray62_, NULL));
        /*
         * fprintf(fid,'Time(hr)\t Runon(cm)\t Ponding(cm)\t Infil(cm)\t Runoff(cm)\t Drain(cm)\t Recharge(cm)\t ET(cm)\t ThetaRZ\t ThetaSZ\t ThetaCZ\n');
         */
        mclAssignAns(
          &ans, mlfNFprintf(0, mclVv(fid, "fid"), _mxarray194_, NULL));
        /*
         * fprintf(fid,'%6.3f\t %6.3f\t %6.3f\t %6.3f\t %6.3f\t %6.3f\t %6.3f\t %6.3f\t %6.3f\t %6.3f\t %.3f\n',WBSum');fclose(fid);
         */
        mclAssignAns(
          &ans,
          mlfNFprintf(
            0,
            mclVv(fid, "fid"),
            _mxarray196_,
            mlfCtranspose(mclVv(WBSum, "WBSum")),
            NULL));
        mclAssignAns(&ans, mlfFclose(mclVv(fid, "fid")));
    /*
     * end
     */
    }
    mclValidateOutput(Dprecip, 1, nargout_, "Dprecip", "recarga");
    mclValidateOutput(*Dimp, 2, nargout_, "Dimp", "recarga");
    mclValidateOutput(*Dperv, 3, nargout_, "Dperv", "recarga");
    mclValidateOutput(*Drunon, 4, nargout_, "Drunon", "recarga");
    mclValidateOutput(*Drunoff, 5, nargout_, "Drunoff", "recarga");
    mclValidateOutput(*Devap, 6, nargout_, "Devap", "recarga");
    mclValidateOutput(*Drecharge, 7, nargout_, "Drecharge", "recarga");
    mclValidateOutput(*Dudrain, 8, nargout_, "Dudrain", "recarga");
    mclValidateOutput(*Dstore, 9, nargout_, "Dstore", "recarga");
    mclValidateOutput(*tctpmax, 10, nargout_, "tctpmax", "recarga");
    mclValidateOutput(*tstponded, 11, nargout_, "tstponded", "recarga");
    mclValidateOutput(*inot, 12, nargout_, "inot", "recarga");
    mxDestroyArray(ans);
    mxDestroyArray(g);
    mxDestroyArray(Sd_imp);
    mxDestroyArray(Si_perv);
    mxDestroyArray(Si_rg);
    mxDestroyArray(Evap_imp);
    mxDestroyArray(Evap_perv);
    mxDestroyArray(Cpan);
    mxDestroyArray(soilsdata);
    mxDestroyArray(NSATsz);
    mxDestroyArray(NRESsz);
    mxDestroyArray(NSATrz);
    mxDestroyArray(NRESrz);
    mxDestroyArray(AA);
    mxDestroyArray(Ss);
    mxDestroyArray(he);
    mxDestroyArray(b);
    mxDestroyArray(rzwilt);
    mxDestroyArray(nwiltfactor);
    mxDestroyArray(nsfactor);
    mxDestroyArray(time);
    mxDestroyArray(dt);
    mxDestroyArray(stime);
    mxDestroyArray(Scz);
    mxDestroyArray(N);
    mxDestroyArray(h0);
    mxDestroyArray(NSATcz);
    mxDestroyArray(NREScz);
    mxDestroyArray(hwf);
    mxDestroyArray(NSAT);
    mxDestroyArray(NRES);
    mxDestroyArray(KSAT);
    mxDestroyArray(thick);
    mxDestroyArray(hin);
    mxDestroyArray(i);
    mxDestroyArray(theta_init);
    mxDestroyArray(Sd_perv);
    mxDestroyArray(Rb);
    mxDestroyArray(Ri);
    mxDestroyArray(Rp);
    mxDestroyArray(Umax);
    mxDestroyArray(udiam);
    mxDestroyArray(hmax);
    mxDestroyArray(Co);
    mxDestroyArray(pi);
    mxDestroyArray(hp_0);
    mxDestroyArray(A_perv);
    mxDestroyArray(Storm_gap);
    mxDestroyArray(t);
    mxDestroyArray(theta_RZ);
    mxDestroyArray(theta_TZ);
    mxDestroyArray(theta_CZ);
    mxDestroyArray(curve_area);
    mxDestroyArray(ponding);
    mxDestroyArray(F);
    mxDestroyArray(Fprev);
    mxDestroyArray(ttp);
    mxDestroyArray(tp);
    mxDestroyArray(f);
    mxDestroyArray(flim);
    mxDestroyArray(Fcum);
    mxDestroyArray(hp);
    mxDestroyArray(roff);
    mxDestroyArray(perc1);
    mxDestroyArray(perc2);
    mxDestroyArray(perc3);
    mxDestroyArray(WBSum);
    mxDestroyArray(Inflow);
    mxDestroyArray(udrain);
    mxDestroyArray(infilt);
    mxDestroyArray(Runoff);
    mxDestroyArray(EvapoT);
    mxDestroyArray(Recharge);
    mxDestroyArray(x);
    mxDestroyArray(s);
    mxDestroyArray(perc1lim);
    mxDestroyArray(perc2lim);
    mxDestroyArray(perc3lim);
    mxDestroyArray(V_inp);
    mxDestroyArray(V_inf);
    mxDestroyArray(V_roff);
    mxDestroyArray(V_rech);
    mxDestroyArray(V_ud);
    mxDestroyArray(V_et);
    mxDestroyArray(V_pond);
    mxDestroyArray(tctponded);
    mxDestroyArray(tctoverflow);
    mxDestroyArray(tstrznsat);
    mxDestroyArray(tsttznsat);
    mxDestroyArray(tctrznsatmax);
    mxDestroyArray(tcttznsatmax);
    mxDestroyArray(tctrznsat);
    mxDestroyArray(tcttznsat);
    mxDestroyArray(tstrzwilt);
    mxDestroyArray(tctrzwilt);
    mxDestroyArray(tctrzwiltmax);
    mxDestroyArray(inst);
    mxDestroyArray(nstimes);
    mxDestroyArray(inpt);
    mxDestroyArray(ptimes);
    mxDestroyArray(oftimes);
    mxDestroyArray(prevof);
    mxDestroyArray(inwt);
    mxDestroyArray(wtimes);
    mxDestroyArray(qs);
    mxDestroyArray(es);
    mxDestroyArray(pdfile);
    mxDestroyArray(len);
    mxDestroyArray(tk);
    mxDestroyArray(CummRain);
    mxDestroyArray(DARain);
    mxDestroyArray(intensity);
    mxDestroyArray(intfile);
    mxDestroyArray(datafile);
    mxDestroyArray(irain);
    mxDestroyArray(evap);
    mxDestroyArray(runon1);
    mxDestroyArray(Vimp);
    mxDestroyArray(rgrain);
    mxDestroyArray(rain);
    mxDestroyArray(Gap_P);
    mxDestroyArray(pdatafile);
    mxDestroyArray(prain);
    mxDestroyArray(runon2);
    mxDestroyArray(Vperv);
    mxDestroyArray(Dbasin);
    mxDestroyArray(runon);
    mxDestroyArray(inp);
    mxDestroyArray(indstorm);
    mxDestroyArray(ttstorm);
    mxDestroyArray(n);
    mxDestroyArray(tcpu);
    mxDestroyArray(test2);
    mxDestroyArray(test3);
    mxDestroyArray(te);
    mxDestroyArray(thetainiRZ);
    mxDestroyArray(thetainiTZ);
    mxDestroyArray(thetainiCZ);
    mxDestroyArray(T);
    mxDestroyArray(test);
    mxDestroyArray(FGA);
    mxDestroyArray(Fp);
    mxDestroyArray(A);
    mxDestroyArray(F1);
    mxDestroyArray(hu);
    mxDestroyArray(uarea);
    mxDestroyArray(uflowmax);
    mxDestroyArray(uflow);
    mxDestroyArray(theta_RZold);
    mxDestroyArray(theta_TZold);
    mxDestroyArray(theta_CZold);
    mxDestroyArray(deficit);
    mxDestroyArray(leftover);
    mxDestroyArray(xswater);
    mxDestroyArray(hpold);
    mxDestroyArray(aw);
    mxDestroyArray(Pevap);
    mxDestroyArray(et);
    mxDestroyArray(DSSW1);
    mxDestroyArray(DSSW2);
    mxDestroyArray(DSSW3);
    mxDestroyArray(Pcompl);
    mxDestroyArray(DSSW);
    mxDestroyArray(unitSM);
    mxDestroyArray(dVsfce);
    mxDestroyArray(Vinsoil);
    mxDestroyArray(infdiff);
    mxDestroyArray(Vr);
    mxDestroyArray(gwrdiff);
    mxDestroyArray(gwravgdepth);
    mxDestroyArray(unitgwr);
    mxDestroyArray(dVsoil);
    mxDestroyArray(unitet);
    mxDestroyArray(unitron);
    mxDestroyArray(unitroff);
    mxDestroyArray(unitud);
    mxDestroyArray(dVs);
    mxDestroyArray(ttcpu);
    mxDestroyArray(outfile);
    mxDestroyArray(fid);
    mxDestroyArray(WBfile);
    mxDestroyArray(TribA);
    mxDestroyArray(Arg);
    mxDestroyArray(EIA);
    mxDestroyArray(CN);
    mxDestroyArray(Urate);
    mxDestroyArray(Kcz);
    mxDestroyArray(Kst);
    mxDestroyArray(Krz);
    mxDestroyArray(Nsoilnumber);
    mxDestroyArray(Ssoilnumber);
    mxDestroyArray(Rsoilnumber);
    mxDestroyArray(St);
    mxDestroyArray(Srz);
    mxDestroyArray(dmin);
    mxDestroyArray(Record);
    mxDestroyArray(Summary);
    mxDestroyArray(outpre);
    mxDestroyArray(rundays);
    mxDestroyArray(avet);
    mxDestroyArray(Sd);
    mxDestroyArray(nprecip);
    mclSetCurrentLocalFunctionTable(save_local_function_table_);
    return Dprecip;
    /*
     * 
     * 
     */
}
