# -*- coding: utf-8 -*-

#******************************************************************************
#
# Freewat
# ---------------------------------------------------------
#
#
# Copyright (C) 2014 - 2015 Iacopo Borsi (iacopo.borsi@tea-group.com)
#
# This source is free software; you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free
# Software Foundation, either version 2 of the License, or (at your option)
# any later version.
#
# This code is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# A copy of the GNU General Public License is available on the World Wide Web
# at <http://www.gnu.org/licenses/>. You can also obtain it by writing
# to the Free Software Foundation, 51 Franklin Street, Suite 500 Boston,
# MA 02110-1335 USA.
#
#******************************************************************************
import os
import time
from PyQt4 import QtGui, uic
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from PyQt4.QtCore import pyqtSlot,SIGNAL,SLOT, pyqtSignal
from qgis.core import *
import numpy as np
from collections import OrderedDict

#from ui_modelBuilder import Ui_ModelBuilderDialog
from freewat.freewat_utils import getVectorLayerByName, getVectorLayerNames, getModelsInfoLists, getModelInfoByName, getModelNlayByName, getModelNspecByName, getTransportModelsByName, fileDialog, getModelUnitsByName, getGroupLayerByName, getUniqueValuesArray, pop_message, ComboStyledItemDelegate
from freewat.sqlite_utils import getTableData, getTableArrayData, checkIfTableExists
import sys
import subprocess as sub

# load flopy and grid utils
from flopy.modflow import *
from flopy.mt3d import *
from flopy.seawat import *

# manage options according to FloPy Version
from flopy import version as vs
flopy_old_vs = True
flopy_latest_vs = False
if vs.__version__ > '3.2.4':
    flopy_latest_vs = True
if vs.__version__ >= '3.2.4':
    from flopy.pakbase import Package
    flopy_old_vs = False
if vs.__version__ == '3.2.3':
    from flopy.mbase import Package
# ---
from flopy.utils import *
import freewat.createGrid_utils
# FloPy Add-On
from freewat.flopyaddon import ModflowFmp, ModflowHob, ModflowFlwob, UcodeWriter, Params, TemplateWriter, ModflowLak, Mt3dmsUSGS, Mt3dUzt, ModflowLmt
#


FORM_CLASS, _ = uic.loadUiType(os.path.join(
    os.path.dirname(__file__), 'ui/ui_modelBuilder.ui'))

class ModelBuilderDialog(QtGui.QDialog, FORM_CLASS):
#class ModelBuilderDialog(QDialog, Ui_ModelBuilderDialog):
    def __init__(self, iface):
        QDialog.__init__(self)
        self.iface = iface
        self.setupUi(self)

        # link the flow model with related transport model(s) - if any:
        self.cmbModelName.currentIndexChanged.connect(self.reloadFields)
        self.boxViscosity.clicked.connect(self.reloadTransportInfo)
        # self.cmbTransportName.currentIndexChanged.connect(self.reloadTransportInfo)
        #
        # self.buttonBoxMF.rejected.connect(self.reject)
        self.buttonBoxMF.button(QDialogButtonBox.Ok).clicked.connect(self.buildModflow)
        self.buttonBoxMF.button(QDialogButtonBox.Ok).setText("Run")
        self.buttonBoxMF.button(QDialogButtonBox.Open).clicked.connect(self.openReportMF)
        self.buttonBoxMF.button(QDialogButtonBox.Open).setText("Open Report")
        #
        self.buttonBoxMT3D.rejected.connect(self.reject)
        self.buttonBoxMT3D.button(QDialogButtonBox.Ok).clicked.connect(self.buildMt3d)
        self.buttonBoxMT3D.button(QDialogButtonBox.Ok).setText("Run")
        self.buttonBoxMT3D.button(QDialogButtonBox.Open).clicked.connect(self.openReportMT)
        self.buttonBoxMT3D.button(QDialogButtonBox.Open).setText("Open Report")
        #
        self.buttonBoxFMP.rejected.connect(self.reject)
        self.buttonBoxFMP.button(QDialogButtonBox.Ok).clicked.connect(self.buildFmp)
        self.buttonBoxFMP.button(QDialogButtonBox.Ok).setText("Run")
        self.buttonBoxFMP.button(QDialogButtonBox.Open).clicked.connect(self.openReportMF)
        self.buttonBoxFMP.button(QDialogButtonBox.Open).setText("Open Report")
        # self.buttonCropModule.clicked.connect(self.cropModule)
        #
        self.buttonBoxUcode.button(QDialogButtonBox.Yes).clicked.connect(self.writeTemplate)
        self.buttonBoxUcode.button(QDialogButtonBox.Yes).setText("Write Template File")
        self.buttonBoxUcode.button(QDialogButtonBox.Save).clicked.connect(self.runUcode)
        self.buttonBoxUcode.button(QDialogButtonBox.Save).setText("Run UCODE")
        self.buttonBoxUcode.button(QDialogButtonBox.Ok).clicked.connect(self.writeUcodeInput)
        self.buttonBoxUcode.button(QDialogButtonBox.Ok).setText("Write Input File")
        self.buttonBoxUcode.button(QDialogButtonBox.Open).setText("Open Report")
        self.buttonBoxUcode.button(QDialogButtonBox.Open).clicked.connect(self.openReportUCODE)
        #

        #
        self.manageGui()
##
##
    def manageGui(self):
        self.cmbModelName.clear()
        layerNameList = getVectorLayerNames()

        layerNameList.sort()

        # Remark: here we assume the name of models table starts with "modeltable"
        (modelNameList, pathList) =  getModelsInfoLists(layerNameList)

# TO DO: inserire qui un messaggio di errore in caso di eccezione
##            else:
##                QMessageBox.warning(self, self.tr('No model'),
##                                self.tr('There is no model table in TOC '
##                                        'You have to create a MODEL before '
##                                        'running Model Layer Creation ' ))

        self.cmbModelName.addItems(modelNameList)
        self.model_name = self.cmbModelName.currentText()
        #
        chd_layers = []
        for nametemp in layerNameList:
            if '_chd' in nametemp:
                chd_layers.append(nametemp)
        self.cmbCHD.addItems(chd_layers)
        self.cmbSSM_CHD.addItems(chd_layers)
        #
        # wel
        wel_layers = []
        for nametemp in layerNameList:
            if '_wel' in nametemp:
                wel_layers.append(nametemp)
        self.cmbWEL.addItems(wel_layers)
        self.cmbSSM_WEL.addItems(wel_layers)
        # mnw
        mnw_layers = []
        for nametemp in layerNameList:
            if '_mnw' in nametemp:
                mnw_layers.append(nametemp)
        self.cmbMNW.addItems(mnw_layers)
        # riv
        riv_layers = []
        for nametemp in layerNameList:
            if '_riv' in nametemp:
                riv_layers.append(nametemp)
        self.cmbRIV.addItems(riv_layers)
        self.cmbSSM_RIV.addItems(riv_layers)

        lak_layers = []
        for nametemp in layerNameList:
            if 'lak_layer_group' in nametemp:
                lak_layers.append(nametemp)
        self.cmbLAK.addItems(lak_layers)

        # ghb
        ghb_layers = []
        for nametemp in layerNameList:
            if '_ghb' in nametemp:
                ghb_layers.append(nametemp)
        self.cmbGHB.addItems(ghb_layers)
        self.cmbSSM_GHB.addItems(ghb_layers)
        # sfr
        sfr_layers = []
        for nametemp in layerNameList:
            if '_sfr' in nametemp:
                sfr_layers.append(nametemp)
        self.cmbSFR.addItems(sfr_layers)
        self.cmbSfrTable.addItems(sfr_layers)
        # rch
        rch_layers = []
        for nametemp in layerNameList:
            if '_rch' in nametemp:
                rch_layers.append(nametemp)
        self.cmbRCH.addItems(rch_layers)
        #
        self.rchopt_list = ['Recharge to top grid', 'Recharge layer defined in irch', 'Recharge to highest active cell']
        self.cmbRchOp.addItems(self.rchopt_list)
        # drn
        drn_layers = []
        for nametemp in layerNameList:
            if '_drn' in nametemp:
                drn_layers.append(nametemp)
        self.cmbDRN.addItems(drn_layers)
        #
        # evt
        evt_layers = []
        for nametemp in layerNameList:
            if '_evt' in nametemp:
                evt_layers.append(nametemp)
        self.cmbEVT.addItems(evt_layers)
        self.evtopt_list = ['ET to top grid', 'ET layer defined in ievt', 'ET to highest active cell']
        self.cmbEvtOp.addItems(self.evtopt_list)
        # uzf
        uzf_layers = []
        surf_layers = []
        for nametemp in layerNameList:
            if '_uzf' in nametemp:
                uzf_layers.append(nametemp)
            if '_sml' in nametemp:
                surf_layers.append(nametemp)
        self.cmbUZF.addItems(uzf_layers)
        self.cmbSurf.addItems(surf_layers)
        #
        # Observations
        # hob
        hob_layers = []
        for nametemp in layerNameList:
            if '_hob' in nametemp:
                hob_layers.append(nametemp)
        self.cmbHOB.addItems(hob_layers)

        # rvob
        rvob_layers = []
        for nametemp in layerNameList:
            if '_rvob' in nametemp:
                rvob_layers.append(nametemp)
        self.cmbRVOB.addItems(rvob_layers)

        # drob
        drob_layers = []
        for nametemp in layerNameList:
            if '_drob' in nametemp:
                drob_layers.append(nametemp)
        self.cmbDROB.addItems(drob_layers)

        # ghob
        gbob_layers = []
        for nametemp in layerNameList:
            if '_gbob' in nametemp:
                gbob_layers.append(nametemp)
        self.cmbGBOB.addItems(gbob_layers)

        # chob
        chob_layers = []
        for nametemp in layerNameList:
            if '_chob' in nametemp:
                chob_layers.append(nametemp)
        self.cmbCHOB.addItems(chob_layers)

        # For MT3DMS
        zoneSSM_layers = []
        for nametemp in layerNameList:
            if '_ssm' in nametemp:
                zoneSSM_layers.append(nametemp)
        self.cmbSSM_Distrib.addItems(zoneSSM_layers)
        self.cmbSSM_Mass.addItems(zoneSSM_layers)
        self.cmbSSM_Constant.addItems(zoneSSM_layers)

        rct_layers = []
        for nametemp in layerNameList:
            if '_rct' in nametemp:
                rct_layers.append(nametemp)
        self.cmbReaction.addItems(rct_layers)

        uzt_layers = []
        for nametemp in layerNameList:
            if '_uzt' in nametemp:
                uzt_layers.append(nametemp)
        self.cmbUZT.addItems(uzt_layers)

        #
        self.cmbAlgorithm.addItems(['Modified Incomplete Cholesky', 'Polynomial'])
        self.cmbMUTPCG.addItems(['Print maximum head change and residual', 'Print only total number of iterations ', 'No printing', 'Print only if convergence fails'])

        # Water Management and crop Modeling  - Farm Process
        FarmIDList = []
        WaterUnitTableList = []
        FarmWellsList = []
        SoilsPropertiesList = []
        CropSoilIDList = []
        CropPropertiesList = []
        Kc_and_rootList = []
        WaterPipelineslist = []
        NRDelivList = []
        GWAllotList = []
        SWAllotList = []
        ClimateList = []
        for nametemp in layerNameList:
            if '_farm_id' in nametemp :
                FarmIDList.append(nametemp)
            if 'WaterDemandSites_' in nametemp:
                WaterUnitTableList.append(nametemp)
            if 'ExternalDeliveries' in nametemp:
                NRDelivList.append(nametemp)
            if 'GW_Allotments' in nametemp:
                GWAllotList.append(nametemp)
            if 'SW_Allotments' in nametemp:
                SWAllotList.append(nametemp)
            if "_farm_well" in nametemp:
                FarmWellsList.append(nametemp)
            if '_soil_crop' in nametemp:
                CropSoilIDList.append(nametemp)
            if "CropsTable_" in nametemp :
                CropPropertiesList.append(nametemp)
            if "CropCoefficients_" in nametemp :
                Kc_and_rootList.append(nametemp)
            if  "_div" in nametemp:
                WaterPipelineslist.append(nametemp)
            if  "ClimateData_" in nametemp:
                ClimateList.append(nametemp)
            if  "SoilsTable_" in nametemp:
                SoilsPropertiesList.append(nametemp)

        self.cmbFarmID.addItems(FarmIDList)
        self.cmbWaterUnitTable.addItems(WaterUnitTableList)
        self.cmbNRDeliv.addItems(NRDelivList)
        self.cmbGWAllot.addItems(GWAllotList)
        self.cmbSWAllot.addItems(SWAllotList)
        self.cmbFarmWells.addItems(FarmWellsList)
        self.cmbCropSoilID.addItems(CropSoilIDList)
        self.cmbCropProperties.addItems(CropPropertiesList)
        self.cmbCropCoefficients.addItems(Kc_and_rootList)
        self.cmbWaterPipelines.addItems(WaterPipelineslist)
        self.cmbClimate.addItems(ClimateList)
        self.cmbSoilsProperties.addItems(SoilsPropertiesList)

        # Use SLM layer already filtered for UZF
        self.cmbEt.addItems(surf_layers)

        # Observation and parameters under UCODE part:
        param_layers = []
        for nametemp in layerNameList:
            if '_params' in nametemp:
                param_layers.append(nametemp)

        self.cmbHob.addItems(hob_layers)
        self.cmbRvob_ucode.addItems(rvob_layers)
        self.cmbGbob_ucode.addItems(gbob_layers)
        self.cmbDrob_ucode.addItems(drob_layers)
        self.cmbChob_ucode.addItems(chob_layers)
        self.cmbParameters.addItems(param_layers)
        self.cmb2DParameters.addItems(param_layers)


##
    def reloadFields(self):
        # Insert list of transport model(s) according with selected flow model:
        self.cmbTransportName.clear()
        layerNameList = getVectorLayerNames()

        self.model_name = self.cmbModelName.currentText()

        try:
            list_of_models = getTransportModelsByName(self.cmbModelName.currentText())
        except:
            list_of_models = []
        self.cmbTransportName.addItems(list_of_models)

    def reloadTransportInfo(self):

        # nspec
        (nspec, mspec, munit) = getModelNspecByName(self.cmbModelName.currentText(), self.cmbTransportName.currentText())

        list_species = [str(i) for i in range(1,nspec+1)]
        # Insert in ComboBox
        model = QtGui.QStandardItemModel(nspec, 1)# nspec rows, 1 col
        it0 = QtGui.QStandardItem('Select Species')
        for i, sp  in enumerate(list_species):
            item = QtGui.QStandardItem(str(sp))
            item.setFlags(Qt.ItemIsUserCheckable | Qt.ItemIsEnabled)
            item.setData(Qt.Checked, Qt.CheckStateRole)
            model.setItem(i+1, 0, item)
        model.setItem(0,it0)

        # Insert Stress Periods in combos
        self.cmbSpecTemp.setModel(model)
        self.cmbSpecTemp.setItemDelegate(ComboStyledItemDelegate())
        self.modelSPEC = model

        ##
    def getProgramLocation(self):
        try:
            # Retrieve programs location
            layerNameList = getVectorLayerNames()
            for layname in layerNameList:
                if layname == 'prg_locations_'+ self.cmbModelName.currentText():
                    locslayer = getVectorLayerByName(layname)
                    dirdict = {}
                    for ft in locslayer.getFeatures():
                        dirdict[ft['code']] = ft['executable']

            return dirdict
        except:
            message  = ''' You didn't enter any location for Executable files !!
                        Open Program Locations and fill in the right path
                        to your codes  '''
            QtGui.QMessageBox.information(None, 'Warning !!', message  )

##
    def reject(self):
        QDialog.reject(self)
##
    def launchProgressBar(self):
        # Define a pop-up Progress Dialog containing a Progress Bar
        self.progressBarDialog = QtGui.QProgressDialog("Collecting Data ...", "Cancel", 0, 0, self)
        self.progressBarDialog.resize(320,50)
        self.progressBarDialog.move(400,200)
        self.progressBar = QtGui.QProgressBar(self.progressBarDialog)
        self.progressBar.setAlignment(Qt.AlignCenter | Qt.AlignVCenter)
        # I set 0.5 h as maximum time for preparing simulation... but this is arbitray!!
        self.maxProgressValue = 1800
        self.progressBar.setRange(0,self.maxProgressValue)
        self.progressBar.setValue(0)
        self.progressBar.move(100,50)
        # Express progress in sec.
        self.progressBar.setFormat("%v {}".format('sec.'))
        self.progressBar.show()
        self.progressBarDialog.setBar(self.progressBar)
        self.progressBarDialog.show()
        self.progressBar.setValue(1)
##
    def updateProgressBar(self, start_time, time_now, string_info = None):
        if string_info is not None :
            self.progressBarDialog.setLabelText(string_info)
        increment = (time_now - start_time)
        self.progressBar.setValue(self.progressBar.value() + increment)

##
    def closeProgressBar(self):
        self.progressBar.setMaximum(self.maxProgressValue)
        self.progressBar.close()
        self.progressBarDialog.reject()
##
    def displayMessage(self, normal_msg, result):
        box = QtGui.QMessageBox()
        box.setWindowTitle('Information')
        box.setDetailedText(result)
        lista = [item.strip() for item in result.split('\n')]

        if normal_msg in lista[-2] :
            txtInfo = '------------------------------------------- \n'+ lista[-2] + '\n------------------------------------------- '
        else:
            txtInfo = 'Simulation was not successful !!\n Click on Show Details to get more information'

        box.setText(txtInfo)
        box.setIcon(1)
        box.setSizePolicy(QtGui.QSizePolicy.Expanding, QtGui.QSizePolicy.Expanding)
        box.exec_()

##
    def execute_query(self, query, params=None):
        """
            Query the database
        """
        from pyspatialite import dbapi2 as db
        self.conn = db.connect(self.dbpath)

        try:
            if params:
                res = self.conn.execute(query, params).fetchall()
            else:
                res = self.conn.execute(query).fetchall()
        except Exception as e:
            #print "Exception"
            raise e

        self.conn.commit()
        self.conn.close()
        self.conn = None

        return res
##
    def buildModflow(self, noRun = False):
        # Define starting time for tracking the process duration,
        #        and launch the progress dialog
        start_time = time.time()
        self.launchProgressBar()

        # ------------ Load input MDO  ------------
        modelName = self.cmbModelName.currentText()
        pathfile, _ = getModelInfoByName(self.model_name)
        self.dbpath = os.path.join(pathfile, '{}.sqlite'.format(self.model_name))
        #
        # Load MDOs of selected Packages
        #
        if self.ckChd.isChecked():
            chdlayer = getVectorLayerByName(self.cmbCHD.currentText())
        #
        if self.ckWel.isChecked():
            wellayer = getVectorLayerByName(self.cmbWEL.currentText())
        #
        if self.ckMnw.isChecked():
            mnwlayer = getVectorLayerByName(self.cmbMNW.currentText())
            tableName =  'mnwtable_' + modelName
            mnwtable  = getVectorLayerByName(tableName)
            if mnwtable == None:
                messageMNWtable = 'You have to create the MNW Table before running MNW! \n' \
                                  'Go to MODFLOW Boundary Conditions > Create MNW Layer'

                QtGui.QMessageBox.warning(None, 'Error in MNW Package', messageMNWtable )
        #
        if self.ckRch.isChecked():
            rchlayer = getVectorLayerByName(self.cmbRCH.currentText())
        #
        if self.ckRiv.isChecked():
            rivlayer = getVectorLayerByName(self.cmbRIV.currentText())
        #
        if self.ckDrn.isChecked():
            drnlayer = getVectorLayerByName(self.cmbDRN.currentText())
        #
        if self.ckGhb.isChecked():
            ghblayer = getVectorLayerByName(self.cmbGHB.currentText())
        #
        if self.ckEvt.isChecked():
            evtlayer = getVectorLayerByName(self.cmbEVT.currentText())
        #
        #
        if self.boxUZF.isChecked():
            uzflayer = getVectorLayerByName(self.cmbUZF.currentText())
        #
        if self.boxSFR.isChecked():
            sfrlayer = getVectorLayerByName(self.cmbSFR.currentText())
            sfrtable = getVectorLayerByName(self.cmbSfrTable.currentText())

        ###--------------------LAK module ------------------------###
        if self.ckLAK.isChecked():
            lak_group = getGroupLayerByName(self.cmbLAK.currentText())
        # ----------
        #
        layerNameList = getVectorLayerNames()

        # ------------ Build a Model ------------
        layNameList = []
        isok = 0
        for mName in layerNameList:
            # Retrieve the model table
            if mName == 'modeltable_'+modelName:
                isok = 1
                modelNameTable = getVectorLayerByName(mName)
                for f in modelNameTable.getFeatures():
                    pathFile =  f['working_dir']
                    lengthString  = f['length_unit']
                    timeString = f["time_unit"]
                    # Convert into MODFLOW flags
                    if lengthString == 'undefined':
                        lenuni  = 0
                    elif lengthString == 'ft':
                        lenuni  = 1
                    elif lengthString == 'm':
                        lenuni  = 2
                    elif lengthString == 'cm':
                        lenuni  = 3

                    if timeString == 'undefined':
                        itmuni = 0
                    elif timeString == 'sec':
                        itmuni = 1
                    elif timeString == 'min':
                        itmuni = 2
                    elif timeString == 'hour':
                        itmuni = 3
                    elif timeString == 'day':
                        itmuni = 4
                    elif timeString == 'year':
                        itmuni = 5

        ## Message Error if no model is found
##            if isok == 0:
##               QMessageBox.warning(self, self.tr('No model found!!'),
##                                self.tr('There is no model table in TOC '
##                                        'You have to create a MODEL before '
##                                        'running Model Layer Creation ' ))
##            # --
            # Retrieve the model time table and data for temporal discretization
            if mName == 'timetable_'+modelName:
                timetable = getVectorLayerByName(mName)

                # number of stress periods
                nper = 0 #timetable.featureCount()

                # Create lists of layers properties
                perlenList = []
                nstpList = []
                tsmultList =[]
                steadyList = []

                for ft in timetable.getFeatures():
                    nper = nper + 1
                    perlenList.append(ft['length'])
                    nstpList.append(ft['ts'])
                    tsmultList.append(ft['multiplier'])
                    state = ft['state']
                    if state == 'SS':
                        steadyList.append(True)
                    else:
                        steadyList.append(False)

                # --
            # Retrieve LPF table, and from there model layers and data needed for LPF
            if mName ==  "lpf_"+ modelName:
                lpftable = getVectorLayerByName(mName)
                # Number of layers
                nlay = 0

                # Get layers name from LPF table
                dpLPF = lpftable.dataProvider()

                # Create lists of layers properties
                layTypeList = []
                layAvgList = []
                layChaniList =[]
                layWetList = []

                for ft in lpftable.getFeatures():
                    attrs = ft.attributes()
                    layNameList.append(attrs[0])
                    layTypeList.append(attrs[1])
                    layAvgList.append(attrs[2])
                    layChaniList.append(attrs[3])
                    layWetList.append(attrs[4])
                    nlay = nlay + 1


                # Translate from keywords to MODFLOW flags, if needed:
                for i, kwtype in enumerate(layTypeList):
                    if  kwtype == 'convertible':
                        layTypeList[i] = 1
                    else:
                        layTypeList[i] = 0
                # --
                for i, kwavg in enumerate(layAvgList):
                    if kwavg == 'harmonic':
                        layAvgList[i] = 0
                    elif kwavg == 'logarithmic':
                        layAvgList[i] = 1
                    else:
                        layAvgList[i] = 2
                # --
                for i, kwet in enumerate(layWetList):
                    if kwet == 'No':
                        layWetList[i] = 0
                    else:
                        layWetList[i] = 1
                # --

                layersList = [layNameList, layTypeList, layAvgList, layChaniList, layWetList]

        # -------- Load MODFLOW Directory
        dirdict = self.getProgramLocation()
        modflowdir = dirdict['MF2005']
        mfnwtdir = dirdict['MF-NWT']
        # pop message if the value of the dictionary is empty
        if not self.radioLmt_unsat.isChecked():
            if dirdict['MF2005'] == '':
                pop_message(self.tr('MODFLOW executable is needed!\nPlease enter the path of such executable in the prg_locations table!'), self.tr('warning'))
                return

        if self.radioLmt_unsat.isChecked() and dirdict['MF-NWT']=='':
            errorMsg = '''
                You are attempting to save groundwater flow also for the
                vadoze zone, but you didn't enter
                MODFLOW NEWTON (MF-NWT) executable in Program Locations,
                to solve your flow model!
                Please enter the path of such executable in the prg_locations table.  '''
            pop_message(self.tr(errorMsg), self.tr('warning'))
            return

        if len(layNameList) == 0:
            pop_message(self.tr('There are no layers defined for the simulation.\nEither the lpf_{} table is missing from the project, or it\'s empty'.format(modelName)), self.tr('warning'))
            return

        # Spatial discretization:
        # Number of rows (along width)
        # new method to get directly nrow and ncol: it works also with rotated grids
        nrow, ncol =  freewat.createGrid_utils.get_row_col(getVectorLayerByName(layNameList[0]))

        # delc, delrow
        delr, delc = freewat.createGrid_utils.get_rgrid_delr_delc(getVectorLayerByName(layNameList[0]))

        # Set TOP and Bottom(s), IBOUND, STRT, KX, KY, KZ
        self.updateProgressBar(start_time, time.time(), string_info = 'Getting data on model layers')
        # Get data from DB
        # check if table exists:
        tbl_exists = checkIfTableExists(self.dbpath,layNameList[0], popupNoExists = True)
        if tbl_exists:
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data on model TOP')
            modeltop = getTableArrayData(self.dbpath, layNameList[0] , col = 'TOP', nrow = nrow, ncol = ncol)
        else:
            return
        # --
        botm = np.zeros(shape = (nlay, nrow, ncol))
        ibound = np.zeros(shape = (nlay, nrow, ncol))
        hstart = np.zeros(shape = (nlay, nrow, ncol))
        kx = np.zeros(shape = (nlay, nrow, ncol))
        ky = np.zeros(shape = (nlay, nrow, ncol))
        kz = np.zeros(shape = (nlay, nrow, ncol))
        ss = np.zeros(shape = (nlay, nrow, ncol))
        sy = np.zeros(shape = (nlay, nrow, ncol))
        wetdry = np.zeros(shape = (nlay, nrow, ncol))

        for i in range(0,nlay):
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data on model layer {}'.format(i+1))
            tbl_exists = checkIfTableExists(self.dbpath,layNameList[i], popupNoExists = True)
            if tbl_exists:
                btemp = getTableArrayData(self.dbpath, layNameList[i] , col = 'BOTTOM' , nrow = nrow, ncol = ncol)
                botm[i, : ,: ] = btemp
                #ibtemp = freewat.createGrid_utils.get_param_array(getVectorLayerByName(layNameList[i]), fieldName = 'ACTIVE')
                ibtemp = getTableArrayData(self.dbpath, layNameList[i] , col = 'ACTIVE' , nrow = nrow, ncol = ncol)
                ibound[i, : ,: ] = ibtemp
                #hstemp = freewat.createGrid_utils.get_param_array(getVectorLayerByName(layNameList[i]), fieldName = 'STRT')
                hstemp = getTableArrayData(self.dbpath, layNameList[i] , col = 'STRT' , nrow = nrow, ncol = ncol)
                hstart[i, :, :] = hstemp
                #kxtemp = freewat.createGrid_utils.get_param_array(getVectorLayerByName(layNameList[i]), fieldName = 'KX')
                kxtemp = getTableArrayData(self.dbpath, layNameList[i] , col = 'KX' , nrow = nrow, ncol = ncol)
                kytemp = getTableArrayData(self.dbpath, layNameList[i] , col = 'KY' , nrow = nrow, ncol = ncol)
                kztemp = getTableArrayData(self.dbpath, layNameList[i] , col = 'KZ' , nrow = nrow, ncol = ncol)
                #kytemp = freewat.createGrid_utils.get_param_array(getVectorLayerByName(layNameList[i]), fieldName = 'KY')
                #kztemp = freewat.createGrid_utils.get_param_array(getVectorLayerByName(layNameList[i]), fieldName = 'KZ')
                #wetdrytemp = freewat.createGrid_utils.get_param_array(getVectorLayerByName(layNameList[i]), fieldName = 'WETDRY')
                wetdrytemp = getTableArrayData(self.dbpath, layNameList[i] , col = 'WETDRY' , nrow = nrow, ncol = ncol)
                #sstemp = freewat.createGrid_utils.get_param_array(getVectorLayerByName(layNameList[i]), fieldName = 'SS')
                sstemp = getTableArrayData(self.dbpath, layNameList[i] , col = 'SS' , nrow = nrow, ncol = ncol)
                #sytemp = freewat.createGrid_utils.get_param_array(getVectorLayerByName(layNameList[i]), fieldName = 'SY')
                sytemp = getTableArrayData(self.dbpath, layNameList[i] , col = 'SY' , nrow = nrow, ncol = ncol)
                # update 3D arrays
                kx[i, :, : ] = kxtemp
                ky[i, :, : ] = kytemp
                kz[i, :, : ] = kztemp
                ss[i,:,:] = sstemp
                sy[i,:,:] = sytemp
                wetdry[i,:,:] = wetdrytemp
            else:
                return

        # Transfor list in array, if necessary
        laytyp = np.zeros(nlay)
        layavg = np.zeros(nlay)
        chani = np.zeros(nlay)
        laywet = np.zeros(nlay)
        layvka = np.zeros(nlay)

        for i in range (0,nlay):
            laytyp[i] = layTypeList[i]
            layavg[i] = layAvgList[i]
            chani[i] = layChaniList[i]
            laywet[i] = layWetList[i]
        #
        # Retrieving layers data finished... update progress bar
        self.updateProgressBar(start_time, time.time())
        # ----
        # Search for MDOs in Legend:
        # ----
        # Retrieve data for CHD: a dictionary
        self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for CHD')
        #
        if self.ckChd.isChecked():
            # Retrieve data for CHD: a triple list
            layer_row_column_CHD = {}
            for i in range(1,nper+1):
                bn = []
                for f in chdlayer.getFeatures():
                    fromlay =  f['from_lay']-1
                    tolay = f['to_lay']-1

                    shead = f[str(i) + '_shead']
                    ehead = f[str(i) + '_ehead']
                    nr = f['row'] - 1
                    nc = f['col'] - 1

                    for k in range(fromlay, tolay+1):
                        bcn = [k, nr, nc, shead, ehead]
                        bn.append(bcn)

                layer_row_column_CHD[i-1] = bn

            #layer_row_column_data = {0 : [ [1, 1, 1, 10, 10 ], [1, 2, 1, 10, 10 ], [1, 3, 1, 10, 10 ] ] }
        # ----
        # Retrieve data for WEL: a triple list
        #
        if self.ckWel.isChecked():
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for WEL')
            layer_row_column_Q = {}
            for i in range(1,nper+1):
                wl = []
                for f in wellayer.getFeatures():
                    fromlay =  f['from_lay'] -1
                    tolay = f['to_lay'] - 1
                    # layer(s) where wel is applied

                    qwel = f['sp_' + str(i) ]
                    nr = f['row'] - 1
                    nc = f['col'] - 1

                    for k in range(fromlay, tolay + 1 ):
                        wlcn = [k, nr, nc, qwel]
                        wl.append(wlcn)

                layer_row_column_Q[i-1] = wl
        # ----
        # ----
        # Retrieve data for MNW:
        #
        if self.ckMnw.isChecked():
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for MNW')
            wellidlist = []
            rowlist = []
            collist = []
            nodeslist = []

            for f in mnwlayer.getFeatures():
                wellidlist.append(f['WELLID'])
                rowlist.append(f['ROW'])
                collist.append(f['COL'])

            mnwmax = len(wellidlist)

            # Define lists as 1D array
            # well ID
            wellidarray = np.empty((mnwmax), dtype='S25')
            for wi, ww in enumerate(wellidlist):
                wellidarray[wi] = ww

            # array of nnodes - initialize
            nnodes = np.zeros((mnwmax), dtype=np.int32)
            # List of arrays for well properties - initialize
            mnw_layer_row_column  = []
            qdes = [0 for i in range(nper)]

            for p in range(1,nper +1):
                # wqdes is a structured array: 1st column a STRING, 2ns column FLOAT
                wqdes = np.empty(shape = (mnwmax, 2) , dtype = 'object')
                wqdes[:,0] = [ww for ww in wellidlist]
                wj = 0
                for f in mnwlayer.getFeatures():
                    wqdes[wj,1] = f['Qw_' + str(p)]
                    wj += 1
                # Updates qdes, list of arrays
                qdes[p-1] = wqdes

            #
            Rwarray = np.zeros(mnwmax )
            Barray = np.zeros(mnwmax )
            Carray = np.zeros(mnwmax )
            Parray = np.zeros(mnwmax )

            for h, w in enumerate(wellidlist):
                ntemp = 0
                # compute nnodes
                for g in mnwtable.getFeatures():
                    if g['well_id'] == w:
                        ntemp += 1
                nnodes[h] = ntemp

                # retreieve data for this well
                wellarray = np.zeros(shape = (ntemp, 3), dtype='int32')

                # Initialize node iterator ni
                ni = 0
                for g in mnwtable.getFeatures():
                    #
                    if g['well_id'] == w:

                        lr = g['layer']
                        rw = rowlist[h]
                        cl = collist[h]
                        #print 'ora --- ', g['well_id'], lr, rw, cl
                        #print 'ni numero nodo ', ni
                        wellarray[ni,0] = lr - 1
                        wellarray[ni,1] = rw - 1
                        wellarray[ni,2] = cl - 1
                        # This is not a smart way of getting parameters: they don't change for each node !!
                        # However, they are repeated in the table, so this methods works. To DO: improving the efficiency!!
                        Rwarray[h] = g['Rw']

                        Barray[h] = g['B']
                        Carray[h] = g['C']
                        Parray[h] = g['P']
                        # Update node iterator
                        ni += 1

                mnw_layer_row_column.append(wellarray)

            # Prepare data to be passed to FloPy:
            losstype = np.empty(mnwmax, dtype='S25')
            losstype[:] = 'GENERAL'
            pumploc = np.zeros(mnwmax, dtype='int32')
            qlimit = np.zeros(mnwmax, dtype='int32')
            ppflag = np.zeros(mnwmax, dtype='int32')
            pumpcap = np.zeros(mnwmax, dtype='int32')
            itmp = np.zeros(nper, dtype='int32')
            itmp[:] = mnwmax
        # ----
        # Retrieve data for RCH: a triple list
        if self.ckRch.isChecked():
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for RCH')
            #
            rech_dict = {}
            irch_dict = {}
            # Select IF rchlayer is in DB (Spatialite) or not:
            if checkIfTableExists(self.dbpath,rchlayer.name()) :
                # This means that RCH MDO is in DB... so Use the Fast Method
                for i in range(1,nper+1):
                    rc  = getTableArrayData(self.dbpath, rchlayer.name() , col = 'sp_' + str(i) + '_rech' , nrow = nrow, ncol = ncol)
                    irc = getTableArrayData(self.dbpath, rchlayer.name() , col = 'sp_' + str(i) + '_irch' , nrow = nrow, ncol = ncol)
                    rech_dict[i-1] = rc
                    irch_dict[i-1] = irc
            else:
                # This means that RCH MDO is NOT in DB... so Use the Old Slower Method
                for i in range(1,nper+1):
                    rc = np.zeros(shape = (nrow,ncol))
                    irc = np.zeros(shape = (nrow,ncol))
                    for f in rchlayer.getFeatures():
                        nr = f['row'] - 1
                        nc = f['col'] - 1
                        rc[nr,nc] = f['sp_' + str(i) + '_rech']
                        irc[nr,nc] = f['sp_' + str(i) + '_irch']

                    rech_dict[i-1] = rc
                    irch_dict[i-1] = irc
            # Retrieve rch_option from GUI:
            optiontext = self.cmbRchOp.currentText()
            if optiontext == self.rchopt_list[0] :
                rch_option = 1
            if optiontext == self.rchopt_list[1] :
                rch_option = 2
            if optiontext == self.rchopt_list[2]:
                rch_option = 3

        # ----
        # ----
        # Retrieve data for RIV:
        #
        if self.ckRiv.isChecked():
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for RIV')
            riv_dict = {}
            #
            for i in range(1,nper+1):
                rv = []
                for f in rivlayer.getFeatures():
                    layrv =  f['layer'] - 1
                    stg = f['stage_' + str(i) ]
                    rb = f['rbot_' + str(i) ]
                    cnd = f['cond_' + str(i) ]
                    nr = f['row'] - 1
                    nc = f['col'] - 1

                    rvlst = [layrv, nr, nc, stg, cnd, rb]
                    rv.append(rvlst)
                riv_dict[i-1] = rv
        # ----
        # ----
        # Retrieve data for DRN:
        #
        if self.ckDrn.isChecked():
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for DRN')
            drn_dict = {}
            for i in range(1,nper+1):
                dr = []
                for f in drnlayer.getFeatures():
                    layrv = f['layer'] - 1
                    elev = f['elev_' + str(i) ]
                    cnd = f['cond_' + str(i) ]
                    nr = f['row'] - 1
                    nc = f['col'] - 1

                    drnlst = [layrv, nr, nc, elev, cnd]
                    dr.append(drnlst)
                drn_dict[i-1] = dr
        #
        # ----
        # ----
        # Retrieve data for GHB:
        #
        if self.ckGhb.isChecked():
            #
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for GHB')
            ghb_dict = {}
            for i in range(1,nper+1):
                gh = []
                for f in ghblayer.getFeatures():
                    fromlay =  f['from_lay'] - 1
                    tolay =  f['to_lay'] - 1
                    bhead = f['bhead_' + str(i) ]
                    cnd = f['cond_' + str(i) ]
                    nr = f['row'] - 1
                    nc = f['col'] - 1

                    for k in range(fromlay, tolay + 1 ):
                        ghblst = [k, nr, nc, bhead, cnd]
                        gh.append(ghblst)

                ghb_dict[i-1] = gh

                # ----
        #

        # ----
        # Retrieve data for EVT:
        #
        if self.ckEvt.isChecked():
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for EVT')
            surf_dict = {}
            evtr_dict = {}
            exdp_dict = {}
            # Select IF evtlayer is in DB (Spatialite) or not:
            if checkIfTableExists(self.dbpath,evtlayer.name()) :
                # This means that EVT MDO is in DB... so Use the Fast Method
                for i in range(1,nper+1):
                    sf  = getTableArrayData(self.dbpath, evtlayer.name() , col = 'surf_' + str(i) , nrow = nrow, ncol = ncol)
                    etr = getTableArrayData(self.dbpath, evtlayer.name() , col = 'evtr_' + str(i) , nrow = nrow, ncol = ncol)
                    edp = getTableArrayData(self.dbpath, evtlayer.name() , col = 'exdp_' + str(i) , nrow = nrow, ncol = ncol)
                    surf_dict[i-1] = sf
                    evtr_dict[i-1] = etr
                    exdp_dict[i-1] = edp
            else:
                # This means that RCH MDO is NOT in DB... so Use the Old Slower Method
                for i in range(1,nper+1):
                    sf = np.zeros(shape = (nrow,ncol))
                    etr = np.zeros(shape = (nrow,ncol))
                    edp = np.zeros(shape = (nrow,ncol))
                    for f in evtlayer.getFeatures():
                        nr = f['row'] - 1
                        nc = f['col'] - 1
                        sf[nr,nc] = f['surf_' + str(i) ]
                        etr[nr,nc] = f['evtr_' + str(i) ]
                        edp[nr,nc] = f['exdp_' + str(i) ]
                    surf_dict[i-1] = sf
                    evtr_dict[i-1] = etr
                    exdp_dict[i-1] = edp

            # Retrieve evt_option from GUI:
            optiontext = self.cmbEvtOp.currentText()
            if optiontext == self.evtopt_list[0] :
                evt_option = 1
            if optiontext == self.evtopt_list[1] :
                evt_option = 2
            if optiontext == self.evtopt_list[2]:
                evt_option = 3

        #

        # ----
        # Retrieve data for UZF:
        #
        if self.boxUZF.isChecked():
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for UFZ')
            #
            iuzfbound = freewat.createGrid_utils.get_param_array(uzflayer, fieldName = 'iuzfbnd')
            eps = freewat.createGrid_utils.get_param_array(uzflayer, fieldName = 'eps')
            thts = freewat.createGrid_utils.get_param_array(uzflayer, fieldName = 'thts')
            thti = freewat.createGrid_utils.get_param_array(uzflayer, fieldName = 'thti')

            finf_list = []
            pet_list = []
            extdp_list = []
            extwc_list = []
            for i in range(1,nper+1):
                ftemp = freewat.createGrid_utils.get_param_array(uzflayer, fieldName = 'finf_'+ str(i))
                finf_list.append(ftemp)
                pttemp = freewat.createGrid_utils.get_param_array(uzflayer, fieldName = 'pet_'+ str(i))
                pet_list.append(pttemp)
                dptemp = freewat.createGrid_utils.get_param_array(uzflayer, fieldName = 'extdp_'+ str(i))
                extdp_list.append(dptemp)
                wctemp = freewat.createGrid_utils.get_param_array(uzflayer, fieldName = 'extwc_'+ str(i))
                extwc_list.append(wctemp)
        # ----
        # ----
        # Retrieve data for SFR2:
        #
        if self.boxSFR.isChecked():
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for SFR2')
            #
            nreach = sfrlayer.featureCount()
            # Count number of segment
            nseglist = [1]
            for fs in sfrlayer.getFeatures():
                if int(fs['seg_id']) not in nseglist:
                    nseglist.append(int(fs['seg_id']))
            nseg = len(nseglist)

            # Initialize reach_data recarray and its columns
            reach_data = np.recarray((nreach,), dtype=[('k', '<i4'), ('i', '<i4'), ('j', '<i4'), ('iseg', '<i4'), ('ireach', '<i4'), ('rchlen', '<f4')])

            for i, f in enumerate(sfrlayer.getFeatures()):
                    reach_data['k'][i]  =  f['layer'] - 1
                    reach_data['i'][i] = f['row'] - 1
                    reach_data['j'][i] = f['col'] - 1
                    reach_data['iseg'][i] = f['seg_id']
                    reach_data['ireach'][i] = f['ireach']
                    reach_data['rchlen'][i] = f['length']

            # Initialize segment_data dictionary
            segment_data = {}
            #sp_list = [n for n in range(nsp)]

            for j in range(1,nper+1):
                # Initialize rec_array for current stress period
                oneseg_data = np.recarray((nseg,),  dtype= [('nseg', '<i4'), ('icalc', '<i4'), ('outseg', '<i4'), ('iupseg', '<i4'), \
                            ('iprior', '<i4'), ('flow', '<f4'), ('runoff', '<f4'), ('etsw', '<f4'), ('pptsw', '<f4'), ('roughch', '<f4'), \
                            ('hcond1', '<f4'), ('thickm1', '<f4'), ('elevup', '<f4'), ('width1', '<f4'), ('thts1', '<f4'), \
                            ('thti1', '<f4'), ('eps1', '<f4'), ('hcond2', '<f4'), ('thickm2', '<f4'), \
                            ('elevdn', '<f4'), ('width2', '<f4'), ('thts2', '<f4'), ('thti2', '<f4'), \
                            ('eps2', '<f4')])
                # ICALC = 1 ,   as default values.
                # DA CAPIRE WIDTH !!!
                oneseg_data['icalc'] = np.ones(shape = (1,nseg))
                oneseg_data['iprior'] = np.ones(shape = (1,nseg))*4
                for i, f in enumerate(sfrtable.getFeatures()):
                    # get data only for current stress period
                    if f['SP'] == j:
                        seg_id = f['SEG_ID'] - 1
                        oneseg_data['nseg'][seg_id] = f['SEG_ID']
                        oneseg_data['outseg'][seg_id] = f['OUT_SEG']
                        oneseg_data['iupseg'][seg_id] = f['UP_SEG']
                        # IPRIOR is necessary only if iupseg NOT 0.... how we can menage this point?
                        oneseg_data['iprior'][seg_id] = f['IPRIOR']
                        oneseg_data['flow'][seg_id] = f['FLOW']
                        oneseg_data['runoff'][seg_id] = f['RUNOFF']
                        oneseg_data['etsw'][seg_id] = f['ETSW']
                        oneseg_data['pptsw'][seg_id] = f['PPTSW']
                        oneseg_data['roughch'][seg_id] = f['ROUGHCH']
                        oneseg_data['hcond1'][seg_id] = f['HCOND1']
                        oneseg_data['thickm1'][seg_id] = f['THICKM1']
                        oneseg_data['elevup'][seg_id] = f['ELEVUP']
                        oneseg_data['width1'][seg_id] = f['WIDTH1']
                        #oneseg_data['depth1'] = f['DEPTH1']
                        oneseg_data['thts1'][seg_id] = f['THTS1']
                        oneseg_data['thti1'][seg_id] = f['THTI1']
                        oneseg_data['eps1'][seg_id] = f['eps1']
                        oneseg_data['hcond2'][seg_id] = f['HCOND2']
                        oneseg_data['thickm2'][seg_id] = f['THICKM2']
                        oneseg_data['elevdn'][seg_id] = f['ELEVDN']
                        oneseg_data['width2'][seg_id] = f['WIDTH2']
                        #oneseg_data['depth2'] = f['DEPTH2']
                        oneseg_data['thts2'][seg_id] = f['THTS2']
                        oneseg_data['thti2'][seg_id] = f['THTI2']
                        oneseg_data['eps2'][seg_id] = f['eps2']

                segment_data[j-1] = oneseg_data
        # end of SFR2
        # ----
        ### ---------- Read lak params from db ----------###
        if self.ckLAK.isChecked():
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for LAK')
            #
            pathfile, _ = getModelInfoByName(self.model_name)
            self.dbpath = os.path.join(pathfile, '{}.sqlite'.format(self.model_name))

            sql = "SELECT * FROM lake_{} WHERE name=?".format(self.model_name)
            params = (self.model_name, )

            res = self.execute_query(sql, params)[0]

            lak_theta = res[1]
            lak_nssitr = res[2]
            lak_sscncr = res[3]

            sql = "SELECT * FROM lak_{} ".format(self.model_name)

            lakes = self.execute_query(sql)

            lak_nlakes = len(lakes)

            lak_surfdep = lakes[0][1]
            lak_stages = []
            lak_stage_range = []

            lak_flux_data = {}
            #lak_leakance = lakes[0][-1]
            for i in range(0, nper):
                lak_flux_data[i] = {}

            for idx, lak in enumerate(lakes):
                lak_stages.append(lak[2])
                lak_stage_range.append([lak[3], lak[4]])

                for i in range(0, nper):
                    sql = "SELECT * FROM laksp_{} WHERE lake_id=? and sp=?".format(self.model_name)

                    sp_data = self.execute_query(sql, (lak[0], i + 1))[0]

                    lak_flux_data[i][idx] = [sp_data[2], sp_data[3], sp_data[4], sp_data[5]]
                    if nper > 1:
                        lak_flux_data[i][idx].append(lak[3])
                        lak_flux_data[i][idx].append(lak[4])

            lak_array = {}
            lak_bdlknc = {}

            for i in range(0, nper):

                lak_array[i] = np.zeros(shape=(nlay, nrow, ncol), dtype=np.int8)
                lak_bdlknc[i] = np.zeros(shape=(nlay, nrow, ncol), dtype=np.float64)

                for l in lak_group.findLayers():
                    layer = l.layer()

                    if not str(layer.name()).endswith('_lak'):
                        continue

                    for f in layer.getFeatures():

                        i_layer = f['layer']
                        i_row = f['row'] - 1
                        i_col = f['col'] - 1

                        lak_array[i][i_layer][i_row][i_col] = f['lake']
                        if lak_array[i][i_layer][i_row][i_col] != 0:
                            lak_bdlknc[i][i_layer][i_row][i_col] = lakes[f['lake'] - 1][-1]  # float(1.0)

        #
        self.updateProgressBar(start_time, time.time(), string_info = 'Building the model')
        # ----
        # -------- New model
        ml = Modflow(modelName, exe_name=modflowdir, version='mf2005', model_ws=pathFile)

        # DIS (Discretization) package
        laycbd = 0

        discret = ModflowDis(ml,nlay=nlay, nrow=nrow, ncol=ncol, nper=nper, delr=delr, delc=delc, laycbd = laycbd, perlen=perlenList, top= modeltop, botm= botm, nstp=nstpList, tsmult=tsmultList, steady=steadyList, itmuni= itmuni , lenuni=lenuni)

        # BAS package
        bas = ModflowBas(ml,ibound=ibound,strt=hstart, hnoflo=-9999.00)

        # LPF (Layer-property Flow) package
        # Remark: by default, layvka = 0.0, so that VKA = KZ, a 3D array od vertical hydraulic conductivity, so vka = kz
        # Retrieve re-wetting parameters from UI

        wetfct= float(self.txtWetfct.text())
        iwetit= int(self.txtIwetit.text())
        ihdwet= int(self.txtIhdwet.text())

        unit_cbcfile = 953
        lpf = ModflowLpf(ml, ipakcb = unit_cbcfile, laytyp = laytyp, layavg = layavg, chani= 1.0, layvka = layvka, laywet= laywet, wetfct=wetfct, iwetit=iwetit, ihdwet=ihdwet, hk=kx, hani= ky, vka= kz, ss= ss, sy= sy, vkcb=0.0, wetdry= wetdry, storagecoefficient= True)

        # CHD (Fixed head boundary condition)
        if self.ckChd.isChecked():
            chd = ModflowChd(ml, stress_period_data = layer_row_column_CHD, cosines=None, extension='chd', unitnumber=24)

        # WEL (WELL package)
        if self.ckWel.isChecked():
            wel = ModflowWel(ml, stress_period_data = layer_row_column_Q )

        # MNW2 (MNW2 Package)
        if self.ckMnw.isChecked():
            # FIXED Values, for each well: Losstype = GENERAL, pumploc = 0, qlimit = 0, pumpcap = 0, ppflag = 0
            mnw = ModflowMnw2(ml, mnwmax= mnwmax , wellid= wellidarray , nnodes= nnodes, losstype= losstype,
                             pumploc= pumploc, qlimit= qlimit, ppflag= ppflag, pumpcap= pumpcap,
                             lay_row_col= mnw_layer_row_column, rw= Rwarray, b= Barray, c= Carray, p= Parray,
                             itmp= itmp, wellid_qdes= qdes)

        # RCH (RCH package)
        if self.ckRch.isChecked():
            rch = ModflowRch(ml, nrchop=rch_option, ipakcb=0, rech = rech_dict, irch = irch_dict, extension='rch', unitnumber=19)

        # RIV (RIVER package)
        if self.ckRiv.isChecked():
            riv = ModflowRiv(ml, stress_period_data = riv_dict )

        ### ---------- Add lak to model ----------###
        if self.ckLAK.isChecked():
            ModflowLak(ml, nlakes=lak_nlakes, ipakcb=0, theta=lak_theta, nssitr=lak_nssitr, sscncr=lak_sscncr, surfdep=lak_surfdep, stages=lak_stages,
                    stage_range=lak_stage_range, flux_data=lak_flux_data, lakarr=lak_array, bdlknc=lak_bdlknc, sill_data={})

        # DRN (DRAIN package)
        if self.ckDrn.isChecked():
            drn = ModflowDrn(ml, stress_period_data = drn_dict )

        # GHB (GHB package)
        if self.ckGhb.isChecked():
            ghb = ModflowGhb(ml, stress_period_data = ghb_dict )

        # EVT (EVAPORATION package)
        if self.ckEvt.isChecked():
            # default value for ievt = 1
            evt = ModflowEvt(ml, nevtop = evt_option , surf =  surf_dict, evtr = evtr_dict, exdp = exdp_dict )


        # UZF package (Unsaturated Zone)
        if self.boxUZF.isChecked():
            #
            nuztopstr = self.cmbUzfOpt.currentText()
            if 'Top' in nuztopstr :
                nuztop = 1
            if 'specified' in nuztopstr :
                nuztop = 2
            if 'Active' in nuztopstr :
                nuztop = 3
            #
            if self.chkUzfSfr.isChecked():
                irunflg = 1
                surflayer = getVectorLayerByName(self.cmbSurf.currentText())
                irunbnd = freewat.createGrid_utils.get_param_array(surflayer, fieldName = 'irunbnd')
            else:
                irunflg = 0
                irunbnd = 0
            #
            if self.chkUzfEvt.isChecked():
                ietflg = 1
            else:
                ietflg = 0
            #
            # FIXED Values: iuzfopt= 2, iuzfcb1=57, iuzfcb2=0, ntrail2=10, nsets=20, , nuzgag=0, surfdep=0.2,
            # So far GAG package is not include .... TO DO !!!
            from flopy import version as vs
            if vs.__version__ < '3.2.6':
                uzf = ModflowUzf1(ml, nuztop=  nuztop, iuzfopt= 2, irunflg= irunflg, ietflg= ietflg, iuzfcb1=57, iuzfcb2=0, ntrail2=10, nsets=20, nuzgag=0, surfdep=0.2,
                                 iuzfbnd= iuzfbound, irunbnd = irunbnd, eps = eps, thts = thts, thtr=0.15, thti = thti, finf= finf_list, pet= pet_list, extdp= extdp_list , extwc= extwc_list )
            else:
                uzf = ModflowUzf1(ml, nuztop=  nuztop, iuzfopt= 2, irunflg= irunflg, ietflg= ietflg, iuzfcb2=0, ntrail2=10, nsets=20, nuzgag=0, surfdep=0.2,
                                 iuzfbnd= iuzfbound, irunbnd = irunbnd, eps = eps, thts = thts, thtr=0.15, thti = thti, finf= finf_list, pet= pet_list, extdp= extdp_list , extwc= extwc_list )



            #uzf = iuzfbnd=1, uzfbud_ext=[], extension='uzf', unitnumber=19)

        # SFR2 package (Stream Flow Routing)
        if self.boxSFR.isChecked():
            if self.radioSfrUnsat.isChecked():
                isfropt = 4
            else:
                isfropt = 0

            # We assume that stream flow routing is always active (nstream < 0)
            nstream = - nreach
            # Inputs from the GUI of ModelBuilder
            const = float(self.txtConstant.text())
            dleak = float(self.txtDleak.text())
            numtim = int(self.txtNumtim.text())
            weight = float(self.txtWeight.text())
            flwtol = float(self.txtFlwtol.text())

            # Data set 5: ITMP IRDFLG IPTFLG (for each stress period)
            dataset_5 = {}
            for i in range(nper):
                dataset_5[i] = [nseg, 0, 0]

            # ATTENTION!!! In FloPy code there is an error (line 825 of mfsr2.py)
            #              irtflg > 0 implies writing rouiting parameters (NUMSTRIM, etc.)
            #              and NOT irtflg < 0 !!! For the moment, we switch irtflg to -1 in our code
            # ATTENTION !!! Starting from flopy 3.2.6, unitnumber attribute has been changed to unit_number
            #              Select this difference according to flopy version:
            from flopy import version as vs
            if vs.__version__ >= '3.2.4':
                sfr = ModflowSfr2(ml,  nstrm = nstream, nss= nseg , const= const, dleak= dleak, isfropt= isfropt,
                         nstrail=10, nsfrsets=30, irtflg = -1, numtim= numtim, weight= weight, flwtol= flwtol,
                         reach_data= reach_data, segment_data= segment_data, dataset_5 = dataset_5, unit_number= 35 )
            else:
                sfr = ModflowSfr2(ml,  nstrm = nstream, nss= nseg , const= const, dleak= dleak, isfropt= isfropt,
                         nstrail=10, nsfrsets=30, irtflg = -1, numtim= numtim, weight= weight, flwtol= flwtol,
                         reach_data= reach_data, segment_data= segment_data, dataset_5 = dataset_5, unitnumber= 35 )

            # Change package name from SFR2 to SFR (this is a bug already solved in dev. version of FloPy!!)
            if vs.__version__ < '3.2.4':
                # FloPy < 3.2.3 bug: in NAM the SFR2 package written SFR2 instead of SFR!
                sfr.name = ['SFR']

        # ---
        # Observations
        # HOB Package (Head Observations)
        if self.ckHob.isChecked():
            hoblayer = getVectorLayerByName(self.cmbHOB.currentText())
            nh = int(hoblayer.featureCount())
            obsnam=[]
            obslayer=[]
            obsrow=[]
            obscolumn=[]
            irefsp=[]
            toffset=[]
            roff=[]
            coff=[]
            observed=[]
            obsfromlay=[]
            obstolay=[]
            maxm = 0
            mlay= nlay*np.ones(shape = (nh))
            pr= [[0 for l in range(nlay)] for i in range(nh)]
            ipr = 0
            mobs = 0
            itt=[1 for i in range(nh) ]
            for fh in hoblayer.getFeatures():
                obsnam.append(fh['OBSNAM'])
                obslayer.append(fh['from_lay'])
                obsfromlay.append(fh['from_lay'])
                obstolay.append(fh['to_lay'])
                obsrow.append(int(fh['ROW']))
                obscolumn.append(int(fh['COL']))
                irefsp.append(int(fh['IREFSP']))
                toffset.append(float(fh['TOFFSET']))
                roff.append(float(fh['ROFF']))
                coff.append(float(fh['COFF']))
                observed.append(float(fh['HOBS']))
                # set up the leayer-percentage array
                if fh['to_lay'] > fh['from_lay']:
                    # update mobs
                    mobs += 1
                    for il, l in enumerate(range(nlay)):
                        obstring = 'PR_%s'%str(il + 1)
                        pr[ipr][il] = float(fh[obstring])
                # update ipr
                ipr += 1

            # update maxm
            if mobs > 0:
                maxm = nlay
            else:
                pr = []

            hob = ModflowHob(ml, nh=nh, mobs=mobs, maxm= maxm, hobdry=0,
                         tomulth=1.0, obsnam= obsnam, layer= obslayer, row= obsrow, column= obscolumn,
                         irefsp= irefsp, toffset= toffset, roff= roff, coff= coff, hob= observed,
                         fromlay= obsfromlay ,tolay= obstolay, mlay= mlay, pr= pr, itt= itt)

        # RVOB Package (RIV Flow Observations)
        if self.ckRvob.isChecked():
            rvoblayer = getVectorLayerByName(self.cmbRVOB.currentText())

            nqcfb_rvob = int(rvoblayer.featureCount()) # This is over estimated, nut nqcfb is only used for max array size so this will work
            # Initialize the recarray representing the attributes table
            data_rvob = np.recarray(shape=(nqcfb_rvob), dtype =[('idcell', int), ('gagename', object ),
                                    ('obsnam', object), ('flwobs', float), ('irefsp', int), ('toffset', float), ('layer', int), ('row', int), ('column', int) ])

            for i, fh in enumerate(rvoblayer.getFeatures()):
                data_rvob['idcell'][i] = i+1
                data_rvob['gagename'][i] = fh['GAGENAME']
                data_rvob['obsnam'][i] = fh['OBSNAM']
                data_rvob['flwobs'][i] = float(fh['FLOWOBS'])
                data_rvob['irefsp'][i] = int(fh['IREFSP'])
                data_rvob['toffset'][i] = float(fh['TOFFSET'])
                data_rvob['layer'][i] = int(fh['layer'])
                data_rvob['row'][i] = int(fh['ROW'])
                data_rvob['column'][i] = int(fh['COL'])

            # Get uniqe values, preserving order
            gagename_rvob = getUniqueValuesArray(data_rvob['gagename'])
            obsnam_rvob   = getUniqueValuesArray(data_rvob['obsnam'])
            zz, rvob_idx   =  np.unique(data_rvob['obsnam'], return_index=True)
            flwobs_rvob   = data_rvob['flwobs'][rvob_idx]
            irefsp_rvob   = data_rvob['irefsp'][rvob_idx]
            toffset_rvob  = data_rvob['toffset'][rvob_idx]

            nqfb_rvob =len(gagename_rvob)
            nqtfb_rvob = len(obsnam_rvob)  # this is correct: total number of ALL observations!!

            # Initialize:
            nqobfb_rvob = [0 for i in range(nqfb_rvob)]
            nqclfb_rvob = [0 for i in range(nqfb_rvob)]

            # Initilize array size
            gg_dict = {}
            for i,gg in enumerate(gagename_rvob):
                gg_dict[gg] = []
                # Count cells and store lay, row, col
                countCells = 0
                for j, ll in enumerate(data_rvob['layer']):
                    if data_rvob['gagename'][j] == gg and [data_rvob['layer'][j], data_rvob['row'][j], data_rvob['column'][j]] not in gg_dict[gg] :
                        countCells += 1
                        gg_dict[gg].append([data_rvob['layer'][j], data_rvob['row'][j], data_rvob['column'][j]])
                nqclfb_rvob[i] = countCells
                # Count obs and store flowbs
                countObs = 0
                oo_list = []
                for k, ff in enumerate(data_rvob['obsnam']):
                    if data_rvob['gagename'][k] == gg and data_rvob['obsnam'][k] not in oo_list:
                        oo_list.append(ff)
                        countObs += 1
                nqobfb_rvob[i] = countObs

            layer_rvob = np.zeros((nqfb_rvob, max(nqclfb_rvob)), dtype='int32')
            row_rvob = np.zeros((nqfb_rvob, max(nqclfb_rvob)), dtype='int32')
            column_rvob = np.zeros((nqfb_rvob, max(nqclfb_rvob)), dtype='int32')
            factor_rvob = np.ones((nqfb_rvob, max(nqclfb_rvob)), dtype='int32')

            for i,gg in enumerate(gg_dict.keys()):
                # Count cells and store lay, row, col
                for j, cell in enumerate(gg_dict[gg]):
                    layer_rvob[i,j] = cell[0]
                    row_rvob[i,j] = cell[1]
                    column_rvob[i,j] = cell[2]

            # Finally, we are in position to
            # define RVOB Package : factor is set 1.0 for ALL cells
            rvob = ModflowFlwob(ml, nqfb= nqfb_rvob, nqcfb= nqcfb_rvob, nqtfb= nqtfb_rvob, tomultfb=1.0, nqobfb= nqobfb_rvob, nqclfb= nqclfb_rvob, obsnam= obsnam_rvob , irefsp= irefsp_rvob ,
                     toffset= toffset_rvob, flwobs= flwobs_rvob, layer= layer_rvob, row= row_rvob, column= column_rvob, factor= factor_rvob, flowtype= 'RIV')

        # DROB Package (DRN Flow Observations)
        if self.ckDrob.isChecked():
            droblayer = getVectorLayerByName(self.cmbDROB.currentText())

            nqcfb_drob = int(droblayer.featureCount()) # This is over estimated, nut nqcfb is only used for max array size so this will work
            # Initialize the recarray representing the attributes table
            data_drob = np.recarray(shape=(nqcfb_drob), dtype =[('idcell', int), ('gagename', object ),
                                    ('obsnam', object), ('flwobs', float), ('irefsp', int), ('toffset', float), ('layer', int), ('row', int), ('column', int) ])

            for i, fh in enumerate(droblayer.getFeatures()):
                data_drob['idcell'][i] = i+1
                data_drob['gagename'][i] = fh['GAGENAME']
                data_drob['obsnam'][i] = fh['OBSNAM']
                data_drob['flwobs'][i] = float(fh['FLOWOBS'])
                data_drob['irefsp'][i] = int(fh['IREFSP'])
                data_drob['toffset'][i] = float(fh['TOFFSET'])
                data_drob['layer'][i] = int(fh['layer'])
                data_drob['row'][i] = int(fh['ROW'])
                data_drob['column'][i] = int(fh['COL'])

            # Get uniqe values, preserving order
            gagename_drob = getUniqueValuesArray(data_drob['gagename'])
            obsnam_drob   = getUniqueValuesArray(data_drob['obsnam'])
            zz, drob_idx   =  np.unique(data_drob['obsnam'], return_index=True)
            flwobs_drob   = data_drob['flwobs'][drob_idx]
            irefsp_drob   = data_drob['irefsp'][drob_idx]
            toffset_drob  = data_drob['toffset'][drob_idx]

            nqfb_drob =len(gagename_drob)
            nqtfb_drob = len(obsnam_drob)  # this is correct: total number of ALL observations!!

            # Initialize:
            nqobfb_drob = [0 for i in range(nqfb_drob)]
            nqclfb_drob = [0 for i in range(nqfb_drob)]

            # Initilize array size
            gg_dict = {}
            for i,gg in enumerate(gagename_drob):
                gg_dict[gg] = []
                # Count cells and store lay, row, col
                countCells = 0
                for j, ll in enumerate(data_drob['layer']):
                    if data_drob['gagename'][j] == gg and [data_drob['layer'][j], data_drob['row'][j], data_drob['column'][j]] not in gg_dict[gg] :
                        countCells += 1
                        gg_dict[gg].append([data_drob['layer'][j], data_drob['row'][j], data_drob['column'][j]])
                nqclfb_drob[i] = countCells
                # Count obs and store flowbs
                countObs = 0
                oo_list = []
                for k, ff in enumerate(data_drob['obsnam']):
                    if data_drob['gagename'][k] == gg and data_drob['obsnam'][k] not in oo_list:
                        oo_list.append(ff)
                        countObs += 1
                nqobfb_drob[i] = countObs

            layer_drob = np.zeros((nqfb_drob, max(nqclfb_drob)), dtype='int32')
            row_drob = np.zeros((nqfb_drob, max(nqclfb_drob)), dtype='int32')
            column_drob = np.zeros((nqfb_drob, max(nqclfb_drob)), dtype='int32')
            # factor is set 1.0 for ALL cells
            factor_drob = np.ones((nqfb_drob, max(nqclfb_drob)), dtype='int32')

            for i,gg in enumerate(gg_dict.keys()):
                # Count cells and store lay, row, col
                for j, cell in enumerate(gg_dict[gg]):
                    layer_drob[i,j] = cell[0]
                    row_drob[i,j] = cell[1]
                    column_drob[i,j] = cell[2]

            # Finally, we are in position to
            # define DROB Package : factor is set 1.0 for ALL cells
            drob = ModflowFlwob(ml, nqfb= nqfb_drob, nqcfb= nqcfb_drob, nqtfb= nqtfb_drob, tomultfb=1.0, nqobfb= nqobfb_drob, nqclfb= nqclfb_drob, obsnam= obsnam_drob , irefsp= irefsp_drob ,
                     toffset= toffset_drob, flwobs= flwobs_drob, layer= layer_drob, row= row_drob, column= column_drob, factor= factor_drob, flowtype= 'DRN')

        # GBOB Package (GHB Flow Observations)
        if self.ckGbob.isChecked():
            gboblayer = getVectorLayerByName(self.cmbGBOB.currentText())

            nqcfb_gbob = int(gboblayer.featureCount()) # This is over estimated, nut nqcfb is only used for max array size so this will work
            # Initialize the recarray representing the attributes table
            data_gbob = np.recarray(shape=(nqcfb_gbob), dtype =[('idcell', int), ('gagename', object ),
                                    ('obsnam', object), ('flwobs', float), ('irefsp', int), ('toffset', float), ('layer', int), ('row', int), ('column', int) ])

            for i, fh in enumerate(gboblayer.getFeatures()):
                data_gbob['idcell'][i] = i+1
                data_gbob['gagename'][i] = fh['GAGENAME']
                data_gbob['obsnam'][i] = fh['OBSNAM']
                data_gbob['flwobs'][i] = float(fh['FLOWOBS'])
                data_gbob['irefsp'][i] = int(fh['IREFSP'])
                data_gbob['toffset'][i] = float(fh['TOFFSET'])
                data_gbob['layer'][i] = int(fh['layer'])
                data_gbob['row'][i] = int(fh['ROW'])
                data_gbob['column'][i] = int(fh['COL'])

            # Get uniqe values, preserving order
            gagename_gbob = getUniqueValuesArray(data_gbob['gagename'])
            obsnam_gbob   = getUniqueValuesArray(data_gbob['obsnam'])
            zz, gbob_idx   = np.unique(data_gbob['obsnam'], return_index=True)
            flwobs_gbob   = data_gbob['flwobs'][gbob_idx]
            irefsp_gbob   = data_gbob['irefsp'][gbob_idx]
            toffset_gbob  = data_gbob['toffset'][gbob_idx]

            nqfb_gbob =len(gagename_gbob)
            nqtfb_gbob = len(obsnam_gbob)  # this is correct: total number of ALL observations!!

            # Initialize:
            nqobfb_gbob = [0 for i in range(nqfb_gbob)]
            nqclfb_gbob = [0 for i in range(nqfb_gbob)]

            # Initilize array size
            gg_dict = {}
            for i,gg in enumerate(gagename_gbob):
                gg_dict[gg] = []
                # Count cells and store lay, row, col
                countCells = 0
                for j, ll in enumerate(data_gbob['layer']):
                    if data_gbob['gagename'][j] == gg and [data_gbob['layer'][j], data_gbob['row'][j], data_gbob['column'][j]] not in gg_dict[gg] :
                        countCells += 1
                        gg_dict[gg].append([data_gbob['layer'][j], data_gbob['row'][j], data_gbob['column'][j]])
                nqclfb_gbob[i] = countCells
                # Count obs and store flowbs
                countObs = 0
                oo_list = []
                for k, ff in enumerate(data_gbob['obsnam']):
                    if data_gbob['gagename'][k] == gg and data_gbob['obsnam'][k] not in oo_list:
                        oo_list.append(ff)
                        countObs += 1
                nqobfb_gbob[i] = countObs

            layer_gbob = np.zeros((nqfb_gbob, max(nqclfb_gbob)), dtype='int32')
            row_gbob = np.zeros((nqfb_gbob, max(nqclfb_gbob)), dtype='int32')
            column_gbob = np.zeros((nqfb_gbob, max(nqclfb_gbob)), dtype='int32')
            # factor is set 1.0 for ALL cells
            factor_gbob = np.ones((nqfb_gbob, max(nqclfb_gbob)), dtype='int32')

            for i,gg in enumerate(gg_dict.keys()):
                # Count cells and store lay, row, col
                for j, cell in enumerate(gg_dict[gg]):
                    layer_gbob[i,j] = cell[0]
                    row_gbob[i,j] = cell[1]
                    column_gbob[i,j] = cell[2]

            # Finally, we are in position to
            # define GBOB Package : factor is set 1.0 for ALL cells
            gbob = ModflowFlwob(ml, nqfb= nqfb_gbob, nqcfb= nqcfb_gbob, nqtfb= nqtfb_gbob, tomultfb=1.0, nqobfb= nqobfb_gbob, nqclfb= nqclfb_gbob, obsnam= obsnam_gbob , irefsp= irefsp_gbob ,
                     toffset= toffset_gbob, flwobs= flwobs_gbob, layer= layer_gbob, row= row_gbob, column= column_gbob, factor= factor_gbob, flowtype= 'GHB')


        # CHOB Package (CHD Flow Observations)
        if self.ckChob.isChecked():
            choblayer = getVectorLayerByName(self.cmbCHOB.currentText())

            nqcfb_chob = int(choblayer.featureCount()) # This is over estimated, nut nqcfb is only used for max array size so this will work
            # Initialize the recarray representing the attributes table
            data_chob = np.recarray(shape=(nqcfb_chob), dtype =[('idcell', int), ('gagename', object ),
                                    ('obsnam', object), ('flwobs', float), ('irefsp', int), ('toffset', float), ('layer', int), ('row', int), ('column', int) ])

            for i, fh in enumerate(choblayer.getFeatures()):
                data_chob['idcell'][i] = i+1
                data_chob['gagename'][i] = fh['GAGENAME']
                data_chob['obsnam'][i] = fh['OBSNAM']
                data_chob['flwobs'][i] = float(fh['FLOWOBS'])
                data_chob['irefsp'][i] = int(fh['IREFSP'])
                data_chob['toffset'][i] = float(fh['TOFFSET'])
                data_chob['layer'][i] = int(fh['layer'])
                data_chob['row'][i] = int(fh['ROW'])
                data_chob['column'][i] = int(fh['COL'])

            # Get uniqe values, preserving order
            gagename_chob = getUniqueValuesArray(data_chob['gagename'])
            obsnam_chob   = getUniqueValuesArray(data_chob['obsnam'])
            zz, chob_idx = np.unique(data_chob['obsnam'], return_index=True)
            flwobs_chob   = data_chob['flwobs'][chob_idx]
            irefsp_chob   = data_chob['irefsp'][chob_idx]
            toffset_chob  = data_chob['toffset'][chob_idx]

            nqfb_chob =len(gagename_chob)
            nqtfb_chob = len(obsnam_chob)  # this is correct: total number of ALL observations!!

            # Initialize:
            nqobfb_chob = [0 for i in range(nqfb_chob)]
            nqclfb_chob = [0 for i in range(nqfb_chob)]

            # Initilize array size
            gg_dict = {}
            for i,gg in enumerate(gagename_chob):
                gg_dict[gg] = []
                # Count cells and store lay, row, col
                countCells = 0
                for j, ll in enumerate(data_chob['layer']):
                    if data_chob['gagename'][j] == gg and [data_chob['layer'][j], data_chob['row'][j], data_chob['column'][j]] not in gg_dict[gg] :
                        countCells += 1
                        gg_dict[gg].append([data_chob['layer'][j], data_chob['row'][j], data_chob['column'][j]])
                nqclfb_chob[i] = countCells
                # Count obs and store flowbs
                countObs = 0
                oo_list = []
                for k, ff in enumerate(data_chob['obsnam']):
                    if data_chob['gagename'][k] == gg and data_chob['obsnam'][k] not in oo_list:
                        oo_list.append(ff)
                        countObs += 1
                nqobfb_chob[i] = countObs

            layer_chob = np.zeros((nqfb_chob, max(nqclfb_chob)), dtype='int32')
            row_chob = np.zeros((nqfb_chob, max(nqclfb_chob)), dtype='int32')
            column_chob = np.zeros((nqfb_chob, max(nqclfb_chob)), dtype='int32')
            # factor is set 1.0 for ALL cells
            factor_chob = np.ones((nqfb_chob, max(nqclfb_chob)), dtype='int32')

            for i,gg in enumerate(gg_dict.keys()):
                # Count cells and store lay, row, col
                for j, cell in enumerate(gg_dict[gg]):
                    layer_chob[i,j] = cell[0]
                    row_chob[i,j] = cell[1]
                    column_chob[i,j] = cell[2]

            # Finally, we are in position to
            # define GBOB Package : factor is set 1.0 for ALL cells
            chob = ModflowFlwob(ml, nqfb= nqfb_chob, nqcfb= nqcfb_chob, nqtfb= nqtfb_chob, tomultfb=1.0, nqobfb= nqobfb_chob, nqclfb= nqclfb_chob, obsnam= obsnam_chob , irefsp= irefsp_chob ,
                     toffset= toffset_chob, flwobs= flwobs_chob, layer= layer_chob, row= row_chob, column= column_chob, factor= factor_chob, flowtype= 'CHD')


        # ------------------------
        # Output control package
##        try:
##            oc = ModflowOc88(ml, unitnumber=[314, 51, 52, unit_cbcfile], item2=[[0,1,1,1]], item3=[[1,0,1,1]])
##        except:
##            # Load the ModflowOc88 if it is not in Flopy anymore (strating from 3.2.6)
##            #from freewat.flopyaddon import ModflowOc88
##            from freewat.flopyaddon.mfoc88freewat import ModflowOc88
##            oc = ModflowOc88(ml, unitnumber=[314, 51, 52, unit_cbcfile], item2=[[0,1,1,1]], item3=[[1,0,1,1]])
        # New for FloPy 3.2.6
        try:
            oc = ModflowOc(ml, unitnumber=[314, 51, 52, 53], compact=True,
               save_every = 1, save_types = ['save head', 'print head','print budget', 'save budget'])
        except:
            from freewat.flopyaddon.mfoc88freewat import ModflowOc88
            oc = ModflowOc88(ml, unitnumber=[314, 51, 52, unit_cbcfile], item2=[[0,1,1,1]], item3=[[1,0,1,1]])

        # PCG package
        mxiter = self.txtOuter.text()
        iter1 = self.txtInner.text()
        hclose = self.txtHclose.text()
        rclose = self.txtRclose.text()
        relax = self.txtRelax.text()

        npcond = 1
        if self.cmbAlgorithm.currentText() == 'Polynomial':
            npocnd = 2

        iprpcg = self.txtIprpcg.text()
        damp = self.txtDamp.text()
        # TO DO: retrieve also MUTPCG from UI
        mutpcg = 3

        pcg = ModflowPcg(ml, mxiter= mxiter, iter1=iter1, npcond=npcond, hclose=hclose, rclose=rclose, relax= relax, nbpol=0, iprpcg= iprpcg, mutpcg=mutpcg, damp=damp, dampt=1.0, ihcofadd=0, extension='pcg', unitnumber=27)

        # Get the MODFLOW code exe:
        modflowexe = modflowdir

        # Create LMT package. Standard flow name is modelname.ftl
        if self.ckLmt.isChecked():
            package_flows = []
            if self.radioLmt_unsat.isChecked():
                # Check if additional package flows are active
                # This uses the new LMT object created by E. Morway, compliant with LMT8
                package_list = ml.get_package_list()
##              if 'sfr' in [x.lower() for x in package_list]:
##                  package_flows.append('SFR')
##              if 'lak' in [x.lower() for x in package_list]:
##                  package_flows.append('LAK')
                if 'uzf' in [x.lower() for x in package_list]:
                    package_flows.append('UZF')
                # And select also MF-NWT as code for the flow the model
                modflowexe = mfnwtdir


            lmt = ModflowLmt(ml,output_file_name= modelName + '.ftl', package_flows = package_flows )

        # End of process to prepare input


        # Proceed in model running or input writing, otherwise only Model Object is saved
        if noRun == False :
            # Write input files
            self.updateProgressBar(start_time, time.time(), string_info = 'Writing input files ...')
            #
            ml.write_input()
            # Close the progress bar for input preparation
            self.updateProgressBar(start_time, time.time(), 'Pre-processing completed!')
            # Sleep for a while ... so the User knows from the ProgressBar that everything was done!
            time.sleep(3)
            self.closeProgressBar()
            #
            if self.chkOnlyInputMF.isChecked():
                QtGui.QMessageBox.information(None, 'Information', 'Input files have been written in: \n' + pathFile )

            # Run model
            # Remark: here the subprocess library is used directly, overpassing
            # the flopy method: ml.run_model(), which seems to fail.
            if not self.chkOnlyInputMF.isChecked():
                QtGui.QMessageBox.information(None, 'Information', 'You are attempting to run: \n' + modflowexe +  '\n' + 'Click OK to start!'  )
                namefile = modelName + '.nam'
                proc = sub.Popen([modflowexe, namefile], stdin = sub.PIPE, stdout = sub.PIPE, stderr = sub.PIPE, cwd=pathFile)
                result = proc.communicate()[0]

                # Write results on LOG file.
                logfileName  = os.path.join(pathFile, modelName + '_MODFLOW_log.txt')
                logfile = open(logfileName, 'w')
                logfile.write(result)
                logfile.close()

                #txtInfo = 'Simulation done!! \n \nSee details on simulation run in log file: \n \n ' + logfileName + '\n \nYou can see the report of results by clicking Open Report'
                # Show message box:
                normal_msg = 'Normal'
                self.displayMessage(normal_msg,result)
                #QtGui.QMessageBox.information(None, 'Information', txtInfo  )

        # Save FloPy Model Object for using it in other parts (Farm Process)
        self.modflowmodel = ml

##
    def buildMt3d(self):
        # Define starting time for tracking the process duration,
        #        and launch the progress dialog
        start_time = time.time()
        self.launchProgressBar()

        # retrieve program location
        # -------- Load MODFLOW Directory
        dirdict = self.getProgramLocation()
        modflowdir = dirdict['MF2005']
        mfnwtdir = dirdict['MF-NWT']
        mt3ddir = dirdict['MT3DMS']
        mtusgsdir = dirdict['MT3D-USGS']
        swtdir = dirdict['SEAWAT']

        # ------------ Load input  Flow Model  ------------
        layerNameList = getVectorLayerNames()
        #
        modelName = self.cmbModelName.currentText()
        transportName = self.cmbTransportName.currentText()
        # Retrieve number of stress periods from FlowModel
        (pathfile, nsp ) = getModelInfoByName(modelName)
        self.dbpath = os.path.join(pathfile, '{}.sqlite'.format(modelName))

        namefile = modelName + '.nam'
        ml = Modflow.load(namefile, version='mf2005', model_ws=pathfile)


        # Retrieve info for compulsory packages BTN and ADV
        # nspec
        (nspec, mspec, munit) = getModelNspecByName(modelName, transportName)

        # nrow, ncol, nlay
        (nrow, ncol, nlay, nper) = ml.get_nrow_ncol_nlay_nper()
        # Units
        (lenuni, itmuni) = getModelUnitsByName(modelName)
        # Convert in string for MT3DMS
        #print '---------- unit of measure', lenuni
        if lenuni == 'ft':
            lunit = 'FT'
        elif lenuni == 'm':
            lunit = 'M'
        else:
            lunit = 'M'
            messageUnit = ''' You are using Length Unit different from m or ft
                              This is not compliant with MT3DMS. Switched to METER
                              Consider revising '''
            QtGui.QMessageBox.warning(None, 'Error in Units of Measure', messageUnit )



        if itmuni == 'hour':
            tunit = 'H'
        elif itmuni == 'day':
            tunit = 'D'
        else:
            tunit = 'D'
            messageUnit = ''' You are using Time Unit different from Hours or Day
                              This is not compliant with MT3DMS. Switched to DAY
                              Consider revising '''
            QtGui.QMessageBox.warning(None, 'Error in Units of Measure', messageUnit )


        # -- flow layers
        for mName in layerNameList:
            if mName == 'lpf_'+modelName :
                lpftable = getVectorLayerByName(mName)
                # Get layers name from LPF table
                dpLPF = lpftable.dataProvider()

                # Create lists of layers properties
                layFlowList = []
                for ft in lpftable.getFeatures():
                    attrs = ft.attributes()
                    layFlowList.append(attrs[0])

        # -- transport layers are named: transportmodelname_lay_n, where n is the layer number
        self.updateProgressBar(start_time, time.time(), string_info = 'Getting data on transport layers')
        layTranspList = []
        for i in range (1, nlay +1):
            layTranspList.append(transportName+'_lay_'+str(i))

        # retrieve ICBUND and PRSITY - also AL and other params for (possible DSP package)
        icbund = np.zeros(shape = (nlay, nrow, ncol))
        prsity = np.zeros(shape = (nlay, nrow, ncol))
        longd  = np.zeros(shape = (nlay, nrow, ncol))
        trpt   = np.zeros(shape = nlay )
        trpv   = np.zeros(shape = nlay)

        for i in range(0,nlay):
            # Get data from DB
            # check if table exists:
            layerT = layTranspList[i]
            tbl_exists = checkIfTableExists(self.dbpath, layerT, popupNoExists = True)
            if tbl_exists:
                ictemp = getTableArrayData(self.dbpath, layerT , col = 'ACTIVE', nrow = nrow, ncol = ncol)
                altemp = getTableArrayData(self.dbpath, layerT , col = 'LONG_D', nrow = nrow, ncol = ncol)
                tttemp = getTableArrayData(self.dbpath, layerT , col = 'TRPT', nrow = nrow, ncol = ncol)
                tvtemp = getTableArrayData(self.dbpath, layerT , col = 'TRPV', nrow = nrow, ncol = ncol)
                prtemp = getTableArrayData(self.dbpath, layFlowList[i] , col = 'NE', nrow = nrow, ncol = ncol)
            else:
                return
            #layerT = getVectorLayerByName(layTranspList[i])
            #
            #ictemp = freewat.createGrid_utils.get_param_array(layerT, fieldName = 'ACTIVE')
            icbund[i, : ,: ] = ictemp
            #
            # altemp = freewat.createGrid_utils.get_param_array(layerT, fieldName = 'LONG_D')
            longd[i, : ,: ] = altemp
            #
            # Remark: aT and aV are stored for each cell (nrow, ncol) BUT actually
            #         there is a unique value for all the layer. We take the MAX just to convert to a scalr
            #tttemp = freewat.createGrid_utils.get_param_array(layerT, fieldName = 'TRPT')
            trpt[i ] = tttemp[0,0]
            #
            #tvtemp = freewat.createGrid_utils.get_param_array(layerT, fieldName = 'TRPV')
            trpv[i ] = tvtemp[0,0]

            # POROSITY from NE in corresponding MODEL LAYER of Flow Model
            #prtemp = freewat.createGrid_utils.get_param_array(getVectorLayerByName(layFlowList[i]), fieldName = 'NE')
            prsity[i, :, : ] = prtemp

        # retrieve SCONC (for each species, for each layer) - also DMCOEF for (possible DSP package)
        sconc_list = []
        dm_list = []
        for j  in range(0,nspec):
            scarray = np.zeros(shape = (nlay, nrow, ncol))
            dmarray = np.zeros(shape = (nlay, nrow, ncol))
            for i in range(0,nlay):
                layerT = layTranspList[i]
                #layerT = getVectorLayerByName(layTranspList[i])
                sctemp = getTableArrayData(self.dbpath, layerT , col = 'SCONC_'+str(j+1), nrow = nrow, ncol = ncol)
                #sctemp = freewat.createGrid_utils.get_param_array(layerT, fieldName = 'SCONC_'+str(j+1))
                dmtemp = getTableArrayData(self.dbpath, layerT, col = 'DMCOEF_'+str(j+1), nrow = nrow, ncol = ncol)
                #dmtemp = freewat.createGrid_utils.get_param_array(layerT, fieldName = 'DMCOEF_'+str(j+1))
                scarray[i, :, : ] = sctemp
                dmarray[i, :, : ] = dmtemp

            sconc_list.append(scarray)
            dm_list.append(dmarray)


        # Retrieve parameters for ADV package
        #
        strMixelm = self.cmbMixelm.currentText()
        if 'Standard' in strMixelm :
            mixelm = 0
        elif  'Forward' in strMixelm:
            mixelm = 1
        elif 'Backward' in strMixelm :
            mixelm = 2
        elif 'Hybrid' in strMixelm :
            mixelm = 3
        elif 'Third' in strMixelm :
            mixelm = -1
        #
        strItrack = self.cmbItrack.currentText()
        # Default: itrack = 1, First-order Eulerian
        itrack = 1
        if strItrack == 'Fourth-order Runge-Kutta':
            itrack = 2
        elif strItrack == 'Mixed First- and Fourth-order':
            itrack = 3

        #
        percel = float(self.txtPercel.text())
        mxpart = float(self.txtMxpart.text())


        # Check if FTL file was created running MODFLOW: if not, MODFLOW is run once again
        # -- here we are assuming that FTL file is named "modelname.ftl"

        if not os.path.isfile(os.path.join(pathfile, modelName + '.ftl')):
        # if not os.path.isfile(pathfile +'\\'+ modelName +'.ftl')  :
            messageNotFTL = ''' You didn't create the Flow Link File!
                                You need to run MODFLOW once again
                                activating the link to LMT package.
                                Please, do it, and come back here later!!  '''

            QtGui.QMessageBox.information(None, 'Information', messageNotFTL )
            return
            # # Create LMT package. Standard flow name is modelname.ftl
            # lmt = ModflowLmt(ml,output_file_name= modelName + '.ftl' )
            # # Write input files
            # # It seems that after loading the model, flopy writes the output files (DATA) twice
            # # ... so, remove the external file names to avoid this bug
            # ml.external_fnames = []
            # #
            # ml.write_name_file()
            # ml.write_input()
            #
            # # Run model
            # # Remark: here the subprocess library is used directly, overpassing
            # # the flopy method: ml.run_model(), which seems to fail.
            #
            # namefile = modelName + '.nam'
            # sub.Popen([modflowdir, namefile], cwd=pathfile )

        # Create MT3DMS model
        mt = Mt3dmsUSGS(modelname= transportName, namefile_ext= 'nam_mt3dms', modflowmodel=ml, ftlfilename= modelName + '.ftl', model_ws= pathfile, exe_name = mt3ddir )

        # Write external Binary files for output
        # mt.external = True
        mt.external_units = [3001]
        mt.external_fnames= [transportName + '.ucn']

        # Add BTN and ADV to Transport Model


        # BTN
        sconc1 = sconc_list[0]
        sconc_dict = {}

        if len(sconc_list) > 1:
            # Multi-component case:
            #sconc2 = sconc_list[1:]
            btn = Mt3dBtn(mt, ncomp= nspec, mcomp= mspec, tunit = tunit, lunit= lunit, munit= munit,
                              prsity= prsity, icbund=icbund, sconc= sconc1, ifmtcn=12, ifmtnp=5,
                              ifmtrf=12, ifmtdp=12, savucn=True, nprs = -1 ,  **sconc_dict )
        else:
            sconc2 = 0
            # Mono-component case
            btn = Mt3dBtn(mt, ncomp= nspec, mcomp= mspec, tunit = tunit, lunit= lunit, munit= munit,
                              prsity= prsity, icbund=icbund, sconc= sconc1, ifmtcn=12, ifmtnp=5,
                              ifmtrf=12, ifmtdp=12, savucn=True, nprs = -1 )

        # ADV
        adv = Mt3dAdv(mt, mixelm= mixelm, percel= percel, mxpart = mxpart, itrack = itrack)


        # Add DSP to Transport Model - IF active
        if self.chkDispersion.isChecked():
            dsp = Mt3dDsp(mt, al= longd, trpt= trpt, trpv= trpv, multiDiff = True, dmcoef= dmarray )


        # Retrieve info for SSM Package - IF active
        flgSSM = 0  # flag for SSM point sources activation
        flgDist = 0 # flag for SSM distributed sources activation
        flgDistMC = 0 # flag for SSM Multi- Component distributed sources activation
        itype = Mt3dSsm.itype_dict()

        # Source in SSM
        ssm_data = {}
        crch_data = {}
        cevt_data = {}

        # If RCH and/or EVT package are active, you have to input the RCH and/or EVT anyway, as 0.0
        # or put incrch = -1 and/or incevt = -1.
        package_list = ml.get_package_list()
        if 'RCH' in package_list:
            crch = 0.0
        else:
            crch = None
        if 'EVT' in package_list:
            cevt = 0.0
        else:
            cevt = None

        # Initialize a list for each stress period
        for i in range(0,nper):
            # Distributed sources  - CRCH e CEVT
            if self.chkSsmDistrib.isChecked():
                self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for SSM (RCH and EVT)')
                flgDist = 1
                distlayer = self.cmbSSM_Distrib.currentText()
                tbl_exists = checkIfTableExists(self.dbpath, distlayer, popupNoExists = True)
                if tbl_exists:
                #distlayer = getVectorLayerByName(self.cmbSSM_Distrib.currentText())

                #crch_array = np.zeros(shape = (nrow, ncol))
                #cevt_array = np.zeros(shape = (nrow, ncol))
                    try:
                        if crch is not None:
                            crch_array = getTableArrayData(self.dbpath, distlayer , col = 'CRCH_sp_%i_spec_1'%(i+1), nrow = nrow, ncol = ncol)
                        if cevt is not None:
                            cevt_array = getTableArrayData(self.dbpath, distlayer , col = 'CEVT_sp_%i_spec_1'%(i+1), nrow = nrow, ncol = ncol)
                    except:
                        pass

##                for f in distlayer.getFeatures():
##                    nr = f['row'] - 1
##                    nc = f['col'] - 1
##                    try:
##                        if crch is not None:
##                            crch_array[nr,nc] = f['CRCH_sp_%i_spec_1'%(i+1)]
##                        if cevt is not None:
##                            cevt_array[nr,nc] = f['CEVT_sp_%i_spec_1'%(i+1)]
##                    except:
##                        pass

##                    try f['CEVT_sp_%i_spec_1'%(i+1)]:
##                        cevt_data[nr,nc] = f['CEVT_sp_%i_spec_1'%(i+1)]
##                    except:
##                        pass


                crch_data[i] = crch_array
                cevt_data[i] = cevt_array

            # Point sources
            ssm_data[i] = []
            # CHD
            if self.chkSsmChd.isChecked():
                self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for SSM (CHD)')
                flgSSM = 1
                chdlayer = getVectorLayerByName(self.cmbSSM_CHD.currentText())
                ssm_chd_list = []
                for f in chdlayer.getFeatures():
                    fromlay =  f['from_lay_'+ str(i+1)] -1
                    tolay =  f['to_lay_'+ str(i+1)] -1
                    nr = f['row'] - 1
                    nc = f['col'] - 1
                    value = []
                    for j in range(1, nspec +1):
                        vl = f['CSSMS_sp_'+ str(i+1) + '_spec_' + str(j) ]
                        value.append(vl)
                    # layer(s) where wel is applied

                    for lay in range(fromlay, tolay +1):
                        tupletemp = (lay, nr, nc, value[0], itype['CHD'] )
                        if nspec > 1:
                            for v in value:
                                tupletemp = tupletemp + (v,)
                        ssm_data[i].append(tupletemp)
            # ---
            # WEL
            if self.chkSsmWel.isChecked():
                self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for SSM (WEL)')
                flgSSM = 1
                wellayer = getVectorLayerByName(self.cmbSSM_WEL.currentText())
                ssm_wel_list = []
                for f in wellayer.getFeatures():
                    try:
                        fromlay =  f['from_lay_'+ str(i+1)] -1
                        tolay =  f['to_lay_'+ str(i+1)] -1
                        nr = f['row'] - 1
                        nc = f['col'] - 1
                        value = []
                        for j in range(1, nspec +1):
                            vl = float(f['CSSMS_sp_'+ str(i+1) + '_spec_' + str(j) ])
                            value.append(vl)
                        # layer(s) where wel is applied
                        for lay in range(fromlay, tolay +1):
                            tupletemp = (lay, nr, nc, value[0], itype['WEL'] )
                            if nspec > 1:
                                for v in value:
                                    tupletemp = tupletemp + (v,)
                            ssm_data[i].append(tupletemp)
                    except:
                        pass
            # ---
            # RIV
            if self.chkSsmRiv.isChecked():
                self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for SSM (RIV)')
                flgSSM = 1
                rivlayer = getVectorLayerByName(self.cmbSSM_RIV.currentText())
                ssm_riv_list = []
                for f in rivlayer.getFeatures():
                    fromlay =  f['from_lay_'+ str(i+1)] -1
                    tolay =  f['to_lay_'+ str(i+1)] -1
                    nr = f['row'] - 1
                    nc = f['col'] - 1
                    value = []
                    for j in range(1, nspec +1):
                        vl = f['CSSMS_sp_'+ str(i+1) + '_spec_' + str(j) ]
                        value.append(vl)
                    # layer(s) where wel is applied
                    for lay in range(fromlay, tolay +1):
                        tupletemp = (lay, nr, nc, value[0], itype['RIV'] )
                        if nspec > 1:
                            for v in value:
                                tupletemp = tupletemp + (v,)
                        ssm_data[i].append(tupletemp)
            # ---
            # GHB
            if self.chkSsmGhb.isChecked():
                self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for SSM (GHB)')
                flgSSM = 1
                ghblayer = getVectorLayerByName(self.cmbSSM_GHB.currentText())
                ssm_ghb_list = []
                for f in ghblayer.getFeatures():
                    fromlay =  f['from_lay_'+ str(i+1)] -1
                    tolay =  f['to_lay_'+ str(i+1)] -1
                    nr = f['row'] - 1
                    nc = f['col'] - 1
                    value = []
                    for j in range(1, nspec +1):
                        vl = f['CSSMS_sp_'+ str(i+1) + '_spec_' + str(j) ]
                        value.append(vl)
                    # layer(s) where wel is applied
                    for lay in range(fromlay, tolay +1):
                        tupletemp = (lay, nr, nc, value[0], itype['GHB'] )
                        if nspec > 1:
                            for v in value:
                                tupletemp = tupletemp + (v,)
                        ssm_data[i].append(tupletemp)
            # ---
            # Mass Loading
            if self.chkSsmMass.isChecked():
                self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for SSM (Mass Loading)')
                flgSSM = 1
                masslayer = getVectorLayerByName(self.cmbSSM_Mass.currentText())
                ssm_mass_list = []
                for f in masslayer.getFeatures():
                    fromlay =  f['from_lay_'+ str(i+1)] -1
                    tolay =  f['to_lay_'+ str(i+1)] -1
                    nr = f['row'] - 1
                    nc = f['col'] - 1
                    value = []
                    for j in range(1, nspec +1):
                        vl = f['CSSMS_sp_'+ str(i+1) + '_spec_' + str(j) ]
                        value.append(vl)
                    # layer(s) where wel is applied
                    for lay in range(fromlay, tolay +1):
                        tupletemp = (lay, nr, nc, value[0], itype['MAS'] )
                        if nspec > 1:
                            for v in value:
                                tupletemp = tupletemp + (v,)
                        ssm_data[i].append(tupletemp)
            # ---
            # Constant Concentration - Dirichelet Condition
            if self.chkSsmConst.isChecked():
                self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for SSM (Const. Conc.)')
                flgSSM = 1
                constlayer = getVectorLayerByName(self.cmbSSM_Constant.currentText())
                ssm_const_list = []
                for f in constlayer.getFeatures():
                    fromlay =  f['from_lay_'+ str(i+1)] -1
                    tolay =  f['to_lay_'+ str(i+1)] -1
                    nr = f['row'] - 1
                    nc = f['col'] - 1
                    value = []
                    for j in range(1, nspec +1):
                        vl = f['CSSMS_sp_'+ str(i+1) + '_spec_' + str(j) ]
                        value.append(vl)
                    # layer(s) where wel is applied
                    for lay in range(fromlay, tolay +1):
                        tupletemp = (lay, nr, nc, value[0], itype['CC'] )
                        if nspec > 1:
                            for v in value:
                                tupletemp = tupletemp + (v,)
                        ssm_data[i].append(tupletemp)

        # For Multi-component (nspec >1)
        # in case of distributed sources
        # we have to define a dictionary of dictionaries for crch and cevt
        if self.chkSsmDistrib.isChecked():
            if nspec > 1:
                flgDistMC = 1
                # Initialize empty dict
                dist_dict = {}
                for j in range(2, nspec + 1):
                    if crch is not None:
                        dist_dict['crch%i'%j] = {}
                    if cevt is not None:
                        dist_dict['cevt%i'%j] = {}
                    for i in range(0,nper):
                        distrch_array = getTableArrayData(self.dbpath, distlayer , col = 'CRCH_sp_%i_spec_%i'%(i+1, j) , nrow = nrow, ncol = ncol)
                        dist_dict['crch%i'%j][i] = distrch_array
                        distevt_array = getTableArrayData(self.dbpath, distlayer , col = 'CEVT_sp_%i_spec_%i'%(i+1, j) , nrow = nrow, ncol = ncol)
                        dist_dict['cevt%i'%j][i] = distevt_array

##                    distrch_array = np.zeros(shape = (nrow, ncol))
##                    distevt_array = np.zeros(shape = (nrow, ncol))
##                    for i in range(0,nper):
##                        for f in distlayer.getFeatures():
##                            nr = f['row'] - 1
##                            nc = f['col'] - 1
##                            try:
##                                distrch_array[nr,nc] = f['CRCH_sp_%i_spec_%i'%(i+1, j)]
##                                dist_dict['crch%i'%j][i] = distrch_array
##                            except:
##                                pass
##                            try:
##                                distevt_array[nr,nc] = f['CEVT_sp_%i_spec_%i'%(i+1, j)]
##                                dist_dict['cevt%i'%j][i] = distevt_array
##                            except:
##                                pass

        # Add SSM package if distributed or point sources are present
        if (flgSSM == 1) or (flgDist == 1):
            mxss = 2*nrow*ncol

            sp_data = {}
            #
            # It seems that FloPy SSM fails IF stress_period_data = None
            # So, in case no point source is present, I create a Zero source
            if nspec == 1:
                for i in range(0, nper ):
                    sp_data[i] = [(0, 1, 1, 0.0, 1)]
            else:
                for i in range(0,nper):
                    mytuple = (0, 1, 1, 0.0, 1, 0.0)
                    for k in range(1,nspec):
                        mytuple += (0.0,)
                    sp_data[i] = [mytuple]
            # if distributed source are present and inoput by User, update crch and cevt
            if flgDist == 1:
                if crch is not None:
                    crch = crch_data
                if cevt is not None:
                    cevt = cevt_data
            # if point sources are present
            if flgSSM == 1:
                sp_data = ssm_data

            # if 1 comp
            #print 'stampa sp_data from SSM ------- ', sp_data, 'flag = ', flgDistMC
            if flgDistMC == 0:
                ssm = Mt3dSsm(mt, crch = crch, cevt = cevt , stress_period_data = sp_data , mxss = mxss)
            # if Multi-component distributed sources are active:
            else:
                ssm = Mt3dSsm(mt, crch = crch, cevt = cevt, stress_period_data = sp_data, mxss = mxss , **dist_dict)

        # ---
        # Add UZT (Unsaturated Zone Transport)
        uztActive = False
        if self.chkUzt.isChecked():
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for UZT')

            # Check that UZF is present in flow model:
            if 'UZF' not in package_list:
                # message error:
                errorMsg = '''
                You are attempting to simulate solute transport in '
                vadoze zone, but you didn't apply UZF package
                in the flow model !!!
                You need to run MODFLOW once again
                activating UZF.
                Please, do it, and come back here later!!  '''

                QtGui.QMessageBox.warning(None, 'No UZF package found ', errorMsg )
                return

            if mfnwtdir == '' :
                # message error:
                errorMsg = '''
                You are attempting to simulate solute transport in
                vadoze zone, but it seems that you didn't enter
                MODFLOW NEWTON (MF-NWT) executable in Program Locations,
                and probably you used standard MODFLOW-2005
                to solve your flow model!
                This will cause an error in your execution  '''

                QtGui.QMessageBox.information(None, 'check MF-NEWTON Executable in Program Locations ', errorMsg )
                # Here no return: in principle the User could delete the executable line AFTER having run the model !!

            if mtusgsdir == '' :
                # message error:
                errorMsg = ''' You are attempting to simulate solute transport in '
                                vadoze zone, but it seems that you didn't enter
                                MT3D-USGS executable in Program Locations!
                                Please, do it, and come back here later!!  '''

                QtGui.QMessageBox.warning(None, 'No MT3D-USGS executable in Program Locations ', errorMsg )
                return

            # If yes, proceed...
            uztActive = True

            # Get info derived from the model
            baspack = ml.get_package('BAS6')
            dispack = ml.get_package('DIS')
            lpfpack = ml.get_package('LPF')

            top = dispack.gettop()
            bottom = dispack.getbotm()
            thickness = dispack.thickness.get_value()
            strt = baspack.strt.get_value()

            # Build WC and SDH, initial water content and saturated thickness, for each layer
            wc = np.zeros(shape = (nlay, nrow, ncol))
            sdh = np.zeros(shape = (nlay, nrow, ncol))

            # Get info from MDO
            uztlayer = getVectorLayerByName(self.cmbUZT.currentText())
            #iuztbnd = np.zeros(shape = (nlay, nrow, ncol))

            iuztbnd = freewat.createGrid_utils.get_param_array(uztlayer, fieldName = 'IUZFBND')
            wctemp = freewat.createGrid_utils.get_param_array(uztlayer, fieldName = 'WC')

            # CUZINF and CUZET (an array for each stress period, for each species.
            # For 1st species: cuzinf is a dictionary runnin on stress periods: {0: array, 1: array, ...}
            # In case of multi-species, the same dictionary, for each species (spec 1, spec 2, etc.)
            #       is stored in a dictionary {'cuzinf2': dictionary, 'cuzinf3': dictionary, ... }
            #       and this is passed as **kwargs to the FloPy object.
            # THE SAME holds true for CUZET

            cuzinf = {}
            cuzet = {}

            for k in range(1,nper+1):
                cuzinf[k-1] = freewat.createGrid_utils.get_param_array(uztlayer, fieldName = 'CUZINF_sp_%s_spec_1'%str(k))
                cuzet[k-1] = freewat.createGrid_utils.get_param_array(uztlayer, fieldName = 'CUZET_sp_%s_spec_1'%str(k))

            # If nspec > 1, Build the dictionary for multi-component, to be passed as kwargs to the FloPy object
            if nspec > 1:
                cuz_multiComponent = {}
                for i in range(2, nspec +1):
                    infMultiComponent = {}
                    etMultiComponent = {}
                    for k in range(1,nper+1):
                        infMultiComponent[k-1] = freewat.createGrid_utils.get_param_array(uztlayer, fieldName = 'CUZINF_sp_%s_spec_%s'%(str(k),str(i)))
                        etMultiComponent[k-1] = freewat.createGrid_utils.get_param_array(uztlayer, fieldName = 'CUZET_sp_%s_spec_%s'%(str(k),str(i)))
                    # include the dictionary in the dictionary:
                    cuz_multiComponent['cuzinf' + str(i)] = infMultiComponent
                    cuz_multiComponent['cuzet' + str(i)] = etMultiComponent



            # Transform in 3D arrays (IUZTBND, WC)
            for k in range(nlay):
                #iuztbnd[k,:,:] = ibndtemp
                wc[k,:,:] = wctemp

                # Compute values for SDH
                for i in range(nrow):
                    for j in range(ncol):
                        if k == 0: # 1st layer:
                            if strt[k][i,j] > bottom[k][i,j]:
                                if strt[k][i,j] >= top[i,j]:
                                    sdh[k,i,j] = thickness[k][i,j]
                                else:
                                    sdh[k,i,j] = strt[k][i,j] - bottom[k][i,j]
                        else:
                            if strt[k][i,j] > bottom[k][i,j]:
                                if strt[k][i,j] >= bottom[k-1][i,j]:
                                    sdh[k][i,j] = thickness[k][i,j]
                                else:
                                    sdh[k][i,j] = strt[k][i,j] - bottom[k][i,j]

            # build the UZT package:
            if nspec > 1:
                uzt = Mt3dUzt( mt, mxuzcon=0, icbcuz=0, iet=0, iuzfbnd= iuztbnd, wc= wc, sdh= sdh,
                               cuzinf = cuzinf, cuzet=cuzet, unitnumber=47, **cuz_multiComponent )
            else:
                uzt = Mt3dUzt( mt, mxuzcon=0, icbcuz=0, iet=0, iuzfbnd= iuztbnd, wc= wc, sdh= sdh,
                               cuzinf = cuzinf, cuzet=cuzet,unitnumber=47)

        # Add RCT (Reaction)
        if self.chkReaction.isChecked():
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for RCT')

            # Retrieve ISOTHM and IREACT flags
            therm = self.cmbSorption.currentText()
            react = self.cmbRate.currentText()

            if therm == 'No sorption':
                isothm = 0
            elif therm == 'Linear isotherm':
                isothm = 1
            elif therm == 'Freundlich isotherm':
                isothm = 2
            elif therm == 'Langmuir isotherm':
                isothm = 3
            elif therm == 'First-order kinetic':
                isothm = 4
            elif therm == 'Dual-porosity without sorption':
                isothm = 5
            elif therm == 'Dual-porosity with sorption':
                isothm = 6

            if react == 'No kinetic rate reaction':
                ireact = 0
            elif react == 'First-order irreversible reaction':
                ireact = 1

            rctlayer = getVectorLayerByName(self.cmbReaction.currentText())

            #rhob_dict = {}
            rhob_array = np.zeros(shape = (nlay, nrow,ncol))
            prsty2_array = np.zeros(shape = (nlay, nrow,ncol))

            #
            #for j in range(1, nspec):
            # Species: SO FAR RCT runs only for 1 COMPONENT !!!!!!
            j = 1
            for i in range(1,nlay+1):
                #rhob = np.zeros(shape = (nrow,ncol))

                sp1_array = np.zeros(shape = (nlay,nrow,ncol))
                sp2_array = np.zeros(shape = (nlay,nrow,ncol))
                rc1_array = np.zeros(shape = (nlay, nrow,ncol))
                rc2_array = np.zeros(shape = (nlay, nrow,ncol))
                # Property for each species is a DICTIONARY of DICTIONARIES
                #sp1_dict[i-1] = {}
##                sp2_dict[i-1] = {}
##                rc1_dict[i-1] = {}
##                rc2_dict[i-1] = {}

##                    sp1_dict = {}
                for f in rctlayer.getFeatures():
                    nr = f['row'] - 1
                    nc = f['col'] - 1
                    if (isothm > 0 ) and (isothm != 5):
                        #rhob[nr,nc] = f['RHOB_lay_'+ str(i)]
                        rhob_array[i-1,nr,nc] = f['RHOB_lay_'+ str(i)]

                    if (isothm == 5 ) or (isothm == 6):
                        prsty2_array[i-1, nr,nc] = f['PRSITY2_lay_' + str(i)]

                    if isothm > 0:
                        sp1_array[i-1,nr,nc] = f['SP1_lay_' + str(i) + '_spec_' + str(j)]
                        sp2_array[i-1,nr,nc] = f['SP2_lay_' + str(i) + '_spec_' + str(j)]

                    if ireact > 0:
                        rc1_array[i-1,nr,nc] = f['RC1_lay_' + str(i) + '_spec_' + str(j)]
                        rc2_array[i-1,nr,nc] = f['RC2_lay_' + str(i) + '_spec_' + str(j)]

                # Update Dictionary for each species
##                sp1_dict[i-1][j-1] = sp1
##                sp2_dict[i-1][j-1] = sp2
##                rc1_dict[i-1][j-1] = rc1
##                rc2_dict[i-1][j-1] = rc2

                # Update dictionary for each layer
                #rhob_dict[i-1] = rhob
                #prsty2_dict[i-1] = prsty2

            # Activate RCT package
            rhob = 0.0
            prsty2 = 0.0
            sp1 = 0.0
            sp2 = 0.0
            rc1 = 0.0
            rc2 = 0.0

            if (isothm > 0 ) and (isothm != 5):
                rhob  = rhob_array
            if (isothm == 5 ) or (isothm == 6):
                prsty2 = prsty2_array
            if isothm > 0:
                sp1 = sp1_array
                sp2 = sp2_array
            if ireact > 0:
                rc1 = rc1_array
                rc2 = rc2_array

            rct = Mt3dRct(mt, isothm= isothm, ireact= ireact, igetsc= 0, rhob= rhob , prsity2 = prsty2, sp1= sp1, sp2 = sp2, rc1 = rc1, rc2 = rc2)


        # Solver GCG
        mxiter = int(self.txtMxiter.text())
        iter1 = int(self.txtIter1.text())
        #
        isolvestr = self.cmbIsolve.currentText()
        if isolvestr == 'Jacobi':
            isolve = 1
        if isolvestr == 'SSOR':
            isolve = 2
        if isolvestr == 'Modified Incomlete Cholesky':
            isolve = 3
        #
        ncrsstr = self.cmbNcrs.currentText()
        if ncrsstr == 'Lumped':
            ncrs =  0
        if ncrsstr == 'Full Dispersion':
            ncrs = 1
        #
        accl = float(self.txtAccl.text())
        #
        cclose = float(self.txtCclose.text())
        #
        iprgcg = int(self.txtIprgcg.text())

        gcg = Mt3dGcg(mt, mxiter= mxiter, iter1= iter1, isolve= isolve, ncrs= ncrs, accl= accl, cclose=1e-05, iprgcg= iprgcg, extension='gcg')

        # Check if SEAWAT is active and proceed accordingly
        seawatFlag = 0

        if self.boxDensity.isChecked() or self.boxViscosity.isChecked():
            seawatFlag = 1
            #
            # manage options according to FloPy Version
            # -- this is needed to deal with a bug in older version of FloPy
            #    concerning format of MT3DMS input files
            #    ONLY WORKING FOR MONO-COMPONENT:
            #    Other limitations: dmcoeff set as default = 1e-9; No Multidiff option
            if flopy_old_vs :
                # Update MT3DMS variables as scalar
                # btn
                newdelr = max(ml.dis.delr)
                newdelc = max(ml.dis.delc)
                htop = max(ml.dis.top[0,:])
                dz = htop - max(ml.dis.botm[0][0,:])
                newprsity = max(mt.btn.prsity[0][0,:])
                newicbund = max(mt.btn.icbund[0][0,:])
                newsconc = sconc1[0][0,0]

                # --- remove the old package and update
                mt.remove_package('BTN')
                btn = Mt3dBtn(mt, delr = newdelr, delc = newdelc, htop = htop, dz = dz, ncomp= nspec, mcomp= mspec, tunit = tunit, lunit= lunit, munit= munit,
                              prsity= newprsity, icbund= newicbund, sconc= newsconc, ifmtcn=12, ifmtnp=5,
                              ifmtrf=12, ifmtdp=12, savucn=True, nprs = -1 )
                # dsp
                newal = max(mt.dsp.al[0][0,:])
                newtrpt = max(mt.dsp.trpt)
                newtrpv = max(mt.dsp.trpv)
                try:
                    mt.remove_package('DSP')
                    dsp = Mt3dDsp(mt, al= newal, trpt= newtrpt, trpv= newtrpv )
                except:
                    pass
            #
            if flopy_latest_vs:
                # For flopy version > 3.2.4, mt3dmsmodel keyword is changed in mt3dmodel !!
                swt = Seawat(modflowmodel= ml, mt3dmodel=mt, modelname= transportName,
                        namefile_ext='nam_swt', model_ws= pathfile, exe_name = swtdir )
            else:
                swt = Seawat(modflowmodel= ml, mt3dmsmodel=mt, modelname= transportName,
                        namefile_ext='nam_swt', model_ws= pathfile, exe_name = swtdir )

        if self.boxDensity.isChecked():
            # Add VDF package
            if self.cmbMfnadvfd.currentText() == 'Central-in-space':
                mfnadvfd = 2
            else:
                mfnadvfd = 1
            vdf = SeawatVdf(swt, mtdnconc= nspec , mfnadvfd= mfnadvfd, nswtcpl=0,
                            iwtable=0, densemin= 0.0, densemax= 0.0,
                            denseref= float(self.txtDenseref.text()), denseslp= float(self.txtDrhodc.text() ))


        if self.boxViscosity.isChecked():
            # Add VSC package
            if self.radioTemperature.isChecked():
                mt3dmuflg = -1
                nsmueos = nspec - 1
                mutempopt = 2
                mtmuspec = []
                for i in range(1,nspec +1):
                    # Looad species number representing Temperature
                    if self.modelSPEC.item(i).checkState() != 2:
                        ntemp = int(self.modelSPEC.item(i).data(0))
                        mtmuspec.append(i)
                viscref = float(self.txtViscref.text())
                dmudc = float(self.txtDmudc.text())
                # Add VSC package to Seawat model object
                vsc = SeawatVsc(swt, mt3dmuflg= mt3dmuflg, viscref= viscref, mutempopt= mutempopt, nsmueos= nsmueos, dmudc= dmudc,
                               cmuref=0, amucoeff=[0.001, 1, 0.015512, -20, -1.572], invisc=-1, visc=-1 )

            else:
                mt3dmuflg = nspec
                mutempopt = 0
                nsmueos= 0
                viscref = float(self.txtViscref.text())
                dmudc = float(self.txtDmudc.text())
                # Add VSC package to Seawat model object
                vsc = SeawatVsc(swt, mt3dmuflg= mt3dmuflg, viscref= viscref, mutempopt= mutempopt, nsmueos= nsmueos, dmudc= dmudc, cmuref=0)

        # Write Input files
        # Set DATA OUTPUT FILES
        #mt.external_fnames = [transportName+'.ucn']
        #mt.external_units = [201]
        # If the new MT3DMS USGS is used, we need to set 0 as UNIT for each package, so that the
            #                            the code assigns automatically the UNITS number
        package_mt_list = mt.get_package_list()
        if 'UZT' in package_mt_list:
            ftl = ml.get_package('FTL')
            #btn.unit_number = [1]
            #btn.unitnumber = [1]
            btn.prsity.locat = 1

        self.updateProgressBar(start_time, time.time(), string_info = 'Writing input files...')
        mt.write_input()
        mt.write_name_file()
        if seawatFlag == 1 :
            swt.write_input()
        self.updateProgressBar(start_time, time.time(), 'Pre-processing completed!')
        # Sleep for a while ... so the User knows from the ProgressBar that everything was done!
        time.sleep(3)
        self.closeProgressBar()

        if self.chkOnlyInputMT.isChecked():
            QtGui.QMessageBox.information(None, 'Information', 'Input files saved in \n' + pathfile )

        # and run

        # MT3DMS
        if not self.chkOnlyInputMT.isChecked() and seawatFlag == 0:

            if uztActive :
                if mtusgsdir == '' :
                    # pop message if the value of the dictionary is empty
                    pop_message(self.tr('MT3D-USGS executable is needed!\nPlease enter the path of such executable in the prg_locations table!'), self.tr('warning'))
                    return
                else:
                    mt3dexe = mtusgsdir
                # Check that TVD (Ultimate) is selected as "Solution Option" (MIXELM) in ADV Package:
                if mixelm > 0:
                    errorMsg = '''
                    You are running a transport model including vadoze zone,
                    but you selected for ADV package the solution option (MIXELM) = %i !
                    This is not allowed in the code MT3D-USGS.
                    Please, select "Standard Finite-Difference" or "Third-order TVD" as Solution option (MIXELM)
                    and run again the model'''%mixelm
                    QtGui.QMessageBox.critical(None, 'Wrong solver selected!!', errorMsg)
                    return

            else:
                if mt3ddir == '' :
                    # pop message if the value of the dictionary is empty
                    pop_message(self.tr('MT3DMS executable is needed!\nPlease enter the path of such executable in the prg_locations table!'), self.tr('warning'))
                    return
                else:
                    mt3dexe = mt3ddir


            QtGui.QMessageBox.information(None, 'Information', 'You are attempting to run: \n' + mt3dexe +  '\n' + 'Click OK to start!' )
            mt_namefile = transportName + '.nam_mt3dms'

            proc = sub.Popen([mt3dexe, mt_namefile], stdin = sub.PIPE, stdout = sub.PIPE, stderr = sub.PIPE, cwd=pathfile)
            result = proc.communicate()[0]
            # Write results on LOG file.
            logfileName  = os.path.join(pathfile, modelName + '_MT3DMS_log.txt')
            logfile = open(logfileName, 'w')
            logfile.write(result)
            logfile.close()

            # Show message box:
            #txtInfo = 'Simulation done!! \n \nSee details on simulation run in log file: \n \n ' + logfileName + '\n \nYou can see the report of results by clicking Open Report'
            normal_msg = 'completed'
            self.displayMessage(normal_msg,result)
            #QtGui.QMessageBox.information(None, 'Information', txtInfo  )


        # or SEAWAT
        if not self.chkOnlyInputMT.isChecked() and seawatFlag == 1:
            if swtdir == '' :
                # pop message if the value of the dictionary is empty
                pop_message(self.tr('SEAWAT executable is needed!\nPlease enter the path of such executable in the prg_locations table!'), self.tr('warning'))
                return

            QtGui.QMessageBox.information(None, 'Information', 'You are attempting to run: \n' + swtdir +  '\n Click OK to start!'  )
            swt_namefile = transportName + '.nam_swt'
            proc = sub.Popen([swtdir, swt_namefile], stdin = sub.PIPE, stdout = sub.PIPE, stderr = sub.PIPE, cwd=pathfile)
            result = proc.communicate()[0]
            # Write results on LOG file.
            logfileName  = os.path.join(pathfile, modelName + '_SEAWAT_log.txt')
            logfile = open(logfileName, 'w')
            logfile.write(result)
            logfile.close()

            #txtInfo = 'Simulation done!! \n \nSee details on simulation run in log file: \n \n ' + logfileName + '\n \nYou can see the report of results by clicking Open Report'
            normal_msg = 'Normal'
            self.displayMessage(normal_msg,result)
            #QtGui.QMessageBox.information(None, 'Information', result  )
            #QtGui.QMessageBox.information(None, 'Information', txtInfo  )



##
    def buildFmp(self):
        # Define starting time for tracking the process duration,
        #        and launch the progress dialog
        start_time = time.time()
        self.launchProgressBar()

        # -------- Load *.EXE Directory
        dirdict = self.getProgramLocation()
        # pop message if the value of the dictionary is empty
        if dirdict['MFOWHM'] == '':
            pop_message(self.tr('MODFLOW-OWHM executable is needed! Please enter the path of such executable in the prg_locations table!'), self.tr('warning'))
            return
        else:
            owhmdir = dirdict['MFOWHM']

        # ------------ Load input  Flow Model  ------------
        layerNameList = getVectorLayerNames()
        #
        modelName = self.cmbModelName.currentText()
        # Retrieve number of stress periods from FlowModel
        (pathfile, nsp ) = getModelInfoByName(modelName)
        namefile = modelName + '.nam'
        # Get the FloPy Modflow Model Object
        ml = Modflow.load(namefile, version='mf2005', model_ws=pathfile)

        # Model properties
        nrow, ncol, nlay, nper = ml.get_nrow_ncol_nlay_nper()

        dis = ml.get_package('DIS')
        gse = dis.gettop()

        # Get input GIS-layers
        farmidlayer = getVectorLayerByName(self.cmbFarmID.currentText())
        waterunitlayer = getVectorLayerByName(self.cmbWaterUnitTable.currentText())
        farmwellslayer = getVectorLayerByName(self.cmbFarmWells.currentText())
        nrdlayer = getVectorLayerByName(self.cmbNRDeliv.currentText())
        gwalotlayer = getVectorLayerByName(self.cmbGWAllot.currentText())
        swalotlayer = getVectorLayerByName(self.cmbSWAllot.currentText())
        cropsoillayer =  getVectorLayerByName(self.cmbCropSoilID.currentText())
        soilsproperties = getVectorLayerByName(self.cmbSoilsProperties.currentText())
        cropslayer = getVectorLayerByName(self.cmbCropProperties.currentText())
        kc_and_rootTable = getVectorLayerByName(self.cmbCropCoefficients.currentText())
        etlayer = getVectorLayerByName(self.cmbEt.currentText())

        climatelayer = getVectorLayerByName(self.cmbClimate.currentText())

        if self.chkPipes.isChecked():
            pipelineslayer = getVectorLayerByName(self.cmbWaterPipelines.currentText())

        # Farms or Water Demand Sites
        nfarms = int(waterunitlayer.featureCount())

        # Define Farm ID array
        self.updateProgressBar(start_time, time.time(), string_info = 'Getting data Water Units distribution')
        fid = np.zeros(shape = (nrow, ncol))
        for f in farmidlayer.getFeatures():
            nr = f['ROW'] - 1
            nc = f['COL'] - 1
            fid[nr, nc] = f['FARM_ID']

        # Count NSOILS and define Soil ID and Crop ID array
        self.updateProgressBar(start_time, time.time(), string_info = 'Getting data on Soil and Crops distribution')
        sid = np.zeros(shape = (nrow, ncol))
        cid = np.zeros(shape = (nrow, ncol))
        #  Initialize soil list dictionary
        soillist = {}
        #  Initialize crops ID list
        croplist = []
        for f in cropsoillayer.getFeatures():
            if int(f['soil_id']) not in soillist.keys() :
                soillist[int(f['soil_id'])] =  [int(f['soil_id']), 0, 'soilname']
            if int(f['crop_id']) not in croplist :
                croplist.append(int(f['crop_id']))
            nr = f['ROW'] - 1
            nc = f['COL'] - 1
            sid[nr, nc] = f['soil_id']
            cid[nr,nc]  = f['crop_id']

        ncrops = len(croplist)
        nsoils = len(soillist)

        # Load data for Water Units Properties
        self.updateProgressBar(start_time, time.time(), string_info = 'Getting data on Water Units properties')
        ofe = {}
        farmcost = {}
        farmid_list = []
        iallotsw = 0
        iallotgw = 0
        inrdfl = 0
        mxnrdt = 0
        for g in waterunitlayer.getFeatures():
            fidtmp = g['Farm_ID']
            farmid_list.append(fidtmp)
            ofe[fidtmp] = g['OFE']
            farmcost[fidtmp] = [float(g['GWBaseCost']), float(g["GWPumpCost"]), float(g["GWLiftCost"]),  float(g["GWDeliveryCost"]),  float(g["SWFixedPrice1"]), float(g["SWLiftCost"]), float(g["SWDeliveryCost"]), float(g["SWFixedPrice2"])]

        # Non-routed deliveries:
        if self.chkNRDeliv.isChecked():
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data on Non-routed Deliv.')
            # dictionary : keys are SP arguments: list of 3-item lists (1 for FARM)
            nonrouted = {}
            inrdfl = 1
            mxnrdt = 1 #nsp*len(farmid_list)
            for cc in nrdlayer.getFeatures():
                sp = int(cc['sp'])
                nrlist = []
                for j, cpid in enumerate(farmid_list):
                    nrlist.append([float(cc['VolNRDeliv_%i'%(j+1)]), float(cc['RankNRDeliv_%i'%(j+1)]), float(cc['TypeNRDeliv_%i'%(j+1)])] )

                nonrouted[sp-1] = nrlist
        else:
            nonrouted = None

        # GW Allotments:
        if self.chkGWAllot.isChecked():
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data on GW Allotments')
            # Initialize dictionary: keys are SP; arguments: list of gw allots (1 for FARM)
            allotgw = {}
            iallotgw = 2
            for cc in gwalotlayer.getFeatures():
                sp = int(cc['sp'])
                gwalist = []
                for j, cpid in enumerate(farmid_list):
                    gwalist.append(float(cc['GW_Allotment_%i'%(j+1)]) )
                allotgw[sp-1] = gwalist
        else:
            allotgw = None

        # SW Allotments
        if self.chkSWAllot.isChecked():
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data on SW Allotments')
            # Initialize dictionary: keys are SP; arguments: list of sw allots (1 for FARM)
            allotsw = {}
            iallotsw = 2
            for cc in swalotlayer.getFeatures():
                sp = int(cc['sp'])
                swalist = []
                for j, cpid in enumerate(farmid_list):
                    swalist.append(float(cc['SW_Allotment_%i'%(j+1)]) )
                allotsw[sp-1] = swalist
        else:
            allotsw = None

        # Load data for Soils Properties
        self.updateProgressBar(start_time, time.time(), string_info = 'Getting data on Soil Properties')
        for s in soilsproperties.getFeatures():
            try:
                soillist[int(s['soil_ID'])][0] = int(s['soil_id'])
                soillist[int(s['soil_ID'])][1] = float(s['CapFringe'])
                soillist[int(s['soil_ID'])][2] = s['soil_name']
            except:
                pass

        # Load data for Crops Properties
        self.updateProgressBar(start_time, time.time(), string_info = 'Getting data on Crops Properties')
        root = {}
        cropevap = {}
        cropfies = {}
        psi = {}
        # cropgrowth = {}
        cropirr = {}
        ifallow = {}
        cropcost = {}
        for g in cropslayer.getFeatures():
            if g['crop_id'] in croplist :
                crpid = int(g['crop_id'])
                # root[crpid] =  float(g['root'])
                cropevap[crpid] = [float(g['FTR']), float(g['FEP']), float(g['FEI'])]
                cropfies[crpid] = [float(g['FIESWP']), float(g['FIESWI'])]
                psi[crpid] = [float(g['psi1']), float(g['psi2']), float(g['psi3']), float(g['psi4'])]
                # cropgrowth[crpid] = [float(g['BaseT']), float(g['MinCutT']), float(g['MaxCutT']),\
                #                      float(g['C0']), float(g['C1']), float(g['C2']), float(g['C3']), float(g['BegRootD']), float(g['MaxRootD']) , float(g['RootGC']), float(g['NONIRR']) ]
                cropirr[crpid] = [int(g['NONIRR'])]

                ifallow[crpid] = float(g['IFALLOW'])
                cropcost[crpid] = [float(g["WPFSlope"]), float(g["WPFInt"]) , float(g["CropPrice"]) ]

        # Initialize dictionaries for kc and rootdepth
        self.updateProgressBar(start_time, time.time(), string_info = 'Getting data on Crop Coefficients')
        kc = {}
        root = {}

        # kc and root
        for cc in kc_and_rootTable.getFeatures():
            sp = int(cc['sp'])
            kclist = []
            rlist = []
            for j, cpid in enumerate(croplist):
                kclist.append(float(cc['kc_%i'%(j+1)]))
                rlist.append(float(cc['root_%i'%(j+1)]))

            kc[sp-1] = kclist
            root[sp-1] = rlist

        # Reference ET (an array for each stress period)
        self.updateProgressBar(start_time, time.time(), string_info = 'Getting data on Ref. ET')
        etref = {}
        for i in range(1,nper+1):
            et = np.zeros(shape = (nrow,ncol))
            for f in etlayer.getFeatures():
                nr = f['row'] - 1
                nc = f['col'] - 1
                et[nr,nc] = f['et_' + str(i) ]
            etref[i-1] = et


        # Semi-Routed Locations - for each farm [Farm-ID Row Column Segment Reach]
        if self.chkPipes.isChecked():
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data on Semi-Routed deliv.')
            isrdfl = 1
            srswlocation = {}
            # Initialize as empty, since we need to set 0s also for farm without Stream connections
            for ff in farmid_list:
                srswlocation[ff] = [ 0, 0 , 0, 0]
            # write real data
            for f in pipelineslayer.getFeatures():
                srid = int(f['farm_id'] )
                srswlocation[srid] = [int(f['ROW']), int(f['COL']), int(f['segment']), int(f['ireach'])]
            #print 'srswlocation = ', srswlocation
        else:
            isrdfl = 0
            srswlocation = None
            # Change iallotsw beacuse it is possible ONLY if SR delivery are defined:
            iallotsw = 0

        # # Climate data -  (ts_step, MaxT, MinT, Precip, ETref)
        # nclimate = climatelayer.featureCount()
        # climatedata = np.recarray((nclimate,), dtype=[('ts_step', float), ('MaxT', float), ('MinT', float), ('Precip', float), ('ETref', float), ('Rad', float)])
        #
        # clid = 0
        # for fc in climatelayer.getFeatures():
        #     climatedata[clid] = ( float(fc['ts_step']), float(fc['MaxT']), float(fc['MinT']), float(fc['Precip']), float(fc['ETref']),float(fc['Rad']) )
        #     clid += 1

        # Precipitation data (constant for each stress period)
        self.updateProgressBar(start_time, time.time(), string_info = 'Getting data on Precipitation')
        precip = []
        for fc in climatelayer.getFeatures():
            precip.append(float(fc['precip']))

        # Optimization options
        ideffl = 0
        if 'Crop' in self.cmbIdeffl.currentText():
            ideffl = -2
        elif 'Economic' in self.cmbIdeffl.currentText():
            ideffl = 1

        # Load data for Farm Wells (for each stress period )
        #
        self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for Farm Wells')
        layer_row_column_QFMP = {}
        mxactw = farmwellslayer.featureCount()
        itmp_farmwells = [mxactw for h in range(nper)]
        for i in range(1,nper+1):
            wlfm = []
            for f in farmwellslayer.getFeatures():
                fromlay =  f['from_lay'] -1
                tolay = f['to_lay'] - 1
                # layer(s) where wel is applied
                qmaxwel = float(f['Qmax_'+ str(i) ])
                wellid = int(f['Well_ID'])
                farmid = int(f['Farm_ID'])
                nr = f['row'] - 1
                nc = f['col'] - 1

                for k in range(fromlay, tolay + 1 ):
                    wlfmtemp = [k, nr, nc, wellid, farmid, qmaxwel  ]
                    wlfm.append(wlfmtemp)

            layer_row_column_QFMP [i-1] = wlfm
        # ----
        # fmp = ModflowFmp(ml, npfwl = 0, mxl = 0, mxactw= mxactw, nfarms= nfarms, ncrops= ncrops, nsoils= nsoils,
        #              ifrmfl = 1, irtfl = 3, icufl = 3, ipfl = 3, iftefl = 1, iieswfl = 1, ieffl = 1,
        #              iebfl = 0, irotfl = 0, ideffl = ideffl, iben = 1, icost = 1, iallotgw = 1,
        #              iccfl = 1, inrdfl = 0, mxnrdt = 0, isrdfl = isrdfl, irdfl = 0, isrrfl = 0, irrfl = -1, iallotsw = 0,
        #              pclose = 0, ifwlcb = 1, ifnrcb = 1, isdpfl = 1, ifbpfl = 2, ietpfl = 2, irtpfl = 0, iopfl = 1, ipapfl = 1, itmp = itmp_farmwells,
        #              farmwells = layer_row_column_QFMP, gse = gse, fid = fid, ofe = ofe, farmcost = farmcost , allotgw = allotgw, srswlocation = srswlocation,
        #              sid = sid, soillist = soillist , cid = cid, root = root, cropevap = cropevap, cropfies = cropfies,
        #              psi = psi, cropgrowth = cropgrowth, ifallow = ifallow, cropcost = cropcost, climate = climatedata)
        # Fixed values for PCLOSE , i.e. SW - allotment water right closure criteria:
        pclose = 0.1
        fmp = ModflowFmp(ml, npfwl = 0, mxl = 0, mxactw= mxactw, nfarms= nfarms, ncrops= ncrops, nsoils= nsoils, ifrmfl = 1, irtfl = 2, icufl = -1, ipfl = 2, iftefl = 1, iieswfl = 1, ieffl = 1, iebfl = 0, irotfl = 0, ideffl = ideffl, iben = 1, icost = 1, iallotgw = iallotgw,
                     iccfl = 1, inrdfl = inrdfl, mxnrdt = mxnrdt, isrdfl = isrdfl, irdfl = 0, isrrfl = 0, irrfl = -1, iallotsw = iallotsw,
                     pclose = pclose, ifwlcb = 1, ifnrcb = 1, isdpfl = 1, ifbpfl = 2, ietpfl = 2, irtpfl = 0, iopfl = 1, ipapfl = 1, itmp = itmp_farmwells,
                     farmwells = layer_row_column_QFMP, gse = gse, fid = fid, ofe = ofe, farmcost = farmcost , nonrouted = nonrouted, allotgw = allotgw, allotsw = allotsw, srswlocation = srswlocation,
                     sid = sid, soillist = soillist , cid = cid, root = root, cropevap = cropevap, cropfies = cropfies,
                     psi = psi, cropirr = cropirr, kc = kc, ifallow = ifallow, etref =etref, precip = precip, cropcost = cropcost)



        # Run the model including FMP
        # Write input files
        self.updateProgressBar(start_time, time.time(), string_info = 'Writing input files...')
        fmp.write_file()

        # It seems that after loading the model, flopy writes the output files (DATA) twice
        # ... so, remove the external file names to avoid this bug
        ml.external_fnames = []
        from flopy import version as vs
        if vs.__version__ < '3.2.4':
            # FloPy < 3.2.3 bug: in NAM the SFR2 package written SFR2 instead of SFR!
            sfr = ml.get_package('SFR2')
            sfr.name = ['SFR']

        # and remove and upload again the OC file, since there is a Unit Number conflict
        ml.remove_package('OC')

        # ------------------------
        # Output control package
##        try:
##            oc = ModflowOc88(ml, unitnumber=[314, 51, 52, 53], item2=[[0,1,1,1]], item3=[[1,0,1,1]])
##        except:
##            # Load the ModflowOc88 if it is not in Flopy anymore (strating from 3.2.6)
##            #from freewat.flopyaddon import ModflowOc88
##            from freewat.flopyaddon.mfoc88freewat import ModflowOc88
##            oc = ModflowOc88(ml, unitnumber=[314, 51, 52, 53], item2=[[0,1,1,1]], item3=[[1,0,1,1]])
        try:
            oc = ModflowOc(ml, unitnumber=[314, 51, 52, 53], compact=True,
               save_every = 1, save_types = ['save head', 'print head','print budget', 'save budget'])
        except:
            from freewat.flopyaddon.mfoc88freewat import ModflowOc88
            oc = ModflowOc88(ml, unitnumber=[314, 51, 52, unit_cbcfile], item2=[[0,1,1,1]], item3=[[1,0,1,1]])

        # then write the new NAM file, with FMP included

        ml.write_name_file()

        self.updateProgressBar(start_time, time.time(), 'Pre-processing completed!')
        # Sleep for a while ... so the User knows from the ProgressBar that everything was done!
        time.sleep(3)
        self.closeProgressBar()


        if self.chkOnlyInputFMP.isChecked():
            QtGui.QMessageBox.information(None, 'Information', 'Input files have been written in: \n' + pathfile )

        # Run model
        # Remark: here the subprocess library is used directly, overpassing
        # the flopy method: ml.run_model(), which seems to fail.
        if not self.chkOnlyInputFMP.isChecked():
            ml.exe_name = owhmdir
            QtGui.QMessageBox.information(None, 'Information', 'You are attempting to run: \n' + owhmdir +  '\n' + 'Click OK to start!'  )
            namefile = modelName + '.nam'
            proc = sub.Popen([owhmdir, namefile], stdin = sub.PIPE, stdout = sub.PIPE, stderr = sub.PIPE, cwd=pathfile)
            result = proc.communicate()[0]
            # Write results on LOG file.
            logfileName  = os.path.join(pathfile, modelName + '_MFOWHM_log.txt')
            logfile = open(logfileName, 'w')
            logfile.write(result)
            logfile.close()

            #txtInfo = 'Simulation done!! \n \nSee details on simulation run in log file: \n \n ' + logfileName + '\n \nYou can see the report of results by clicking Open Report'
            normal_msg = 'Normal'
            self.displayMessage(normal_msg,result)
            #QtGui.QMessageBox.information(None, 'Information', result  )
            #QtGui.QMessageBox.information(None, 'Information', txtInfo  )

#
    # # Run Crop Growth Module
    # def cropModule(self):
    #     import crop_dialog as cropDialog
    #     dlg = cropDialog.createCropDialog(self.iface)
    #     # show the dialog
    #     dlg.show()
    #     # Run the dialog event loop
    #     dlg.exec_()

#
    def writeUcodeInput(self):
        # Define starting time for tracking the process duration,
        #        and launch the progress dialog
        start_time = time.time()
        self.launchProgressBar()
        #
        modelName = self.cmbModelName.currentText()
        dirdict = self.getProgramLocation()
        # pop message if the value of the dictionary is empty
        if dirdict['UCODE'] == '':
            pop_message(self.tr('UCODE executable is needed!\nPlease enter the path of such executable in the prg_locations table!'), self.tr('warning'))
            return
        else:
            ucodedir = dirdict['UCODE']
        if dirdict['MF2005'] == '':
            pop_message(self.tr('MF2005 executable is needed!\nPlease enter the path of such executable in the prg_locations table!'), self.tr('warning'))
            return
        else:
            modflowdir = dirdict['MF2005']

        #Get basic info of MODFLOW model from modelname
        # path and stress periods
        (pathfile, nsp ) = getModelInfoByName(modelName)
        # MODFLOW model FloPy object
        namefile = modelName + '.nam'
        ml = Modflow.load(namefile, version='mf2005', model_ws=pathfile)
        # nrow, ncol, nlay
        (nrow, ncol, nlay, nper) = ml.get_nrow_ncol_nlay_nper()
        # Units
        (lenuni, itmuni) = getModelUnitsByName(modelName)
        #
        modelname = modelName
        modellengthunits= lenuni
        modelmassunits= 'kg' # This will be get from Transport Model, but so far this part is not running
        modeltimeunits= itmuni

        # ----- Get UCODE options from GUI
        verbose = self.cmbVerbose.currentText()
        if '0' in verbose:
            verbose = 0
        elif '1' in verbose:
            verbose = 1
        elif '2' in verbose:
            verbose = 2
        elif '3' in verbose:
            verbose = 3
        elif '4' in verbose:
            verbose = 4
        else:
            verbose = 5
        #
        sensitivities='no'
        if self.chkSensitivities.isChecked():
            sensitivities='yes'
        optimize='no'
        if self.chkOptimize.isChecked():
            optimize = 'yes'
        eigenvalues='no'
        if self.chkEigenValues.isChecked():
            eigenvalues='yes'
        startres='no'
        if self.chkStartRes.isChecked():
            startres='yes'
        intermedres='no'
        if self.chkIntermedRes.isChecked():
            intermedres='yes'
        finalres='no'
        if self.chkIntermedFinalRes.isChecked():
            finalres='yes'
        dataexchange='yes'
        tolpar = float(self.txtTolpar.text() ) # 0.01
        tolsosc= float(self.txtTolsocs.text()) # 0.001
        maxiter= int(self.txtMaxiter.text()) # 10
        xchange= float(self.txtMaxchange.text())  # 2.0
        senmethod = 1
        if 'Central' in self.cmbSenMethod.currentText():
            senmethod = 2
        trustregion='no'
        if self.chkTrustregion.isChecked():
            trustregion='yes'
        command= '"'+modflowdir + '  '  +  os.path.join(pathfile, modelName + '.nam') + '"' # this needs to be the command to run the actual model (name file included), so likely a batch file

        commandid='modflow'

        # ---- Get data for parameters, from both 3D and 2D parameters Table
        npar = 0
        modinfile = []
        templatefile = []
        parname = []
        mfpackage = []
        startvalue =[]
        adjustable = []
        transform= []
        mfpackuniq = []
        constrain=[]
        lowerconstraint=[]
        upperconstraint=[]
        perturbamt=[]

        #  Data from 3D parameters Table
        if self.chk3Dparams.isChecked():
            self.updateProgressBar(start_time, time.time(), 'Get data for LPF parameters')
            partable = getVectorLayerByName(self.cmbParameters.currentText())
            for f in partable.getFeatures():
                npar += 1
                parname.append(f['Name'])
                mfpackage.append(f['Package'])
                startvalue.append(float(f['StartValue']))
                adjustable.append(f['Adjustable'])
                transform.append(f['Transform'])
                constrain.append(str(f['Constr']).lower())
                lowerconstraint.append(f['Lower'])
                upperconstraint.append(f['Upper'])
                perturbamt.append(float(f['PerturbAmt']))

        # ---- Get data for parameters, from 2D time-variant parameters Table
        if self.chk2Dparams.isChecked():
            self.updateProgressBar(start_time, time.time(), 'Get data for RCH-EVT parameters')
            partable = getVectorLayerByName(self.cmb2DParameters.currentText())
            #
            for f in partable.getFeatures():
                name_temp = f['Name']
                if name_temp not in parname:
                    parname.append(name_temp)
                    mfpackage.append(f['Package'])
                    startvalue.append(float(f['StartValue']))
                    adjustable.append(f['Adjustable'])
                    transform.append(f['Transform'])
                    # For 2D Params: constr, lower, upper, perturbamt are NOT working so far (Dec. 2016), but needed by ucode_in.py
                    constrain.append('no')
                    lowerconstraint.append(0.0)
                    upperconstraint.append(0.0)
                    perturbamt.append(0.01)
                    npar += 1

        # --- set up names for template files based on parameterized packages
        mfpackuniq = np.unique(mfpackage) # Get the unique package names
        for packname in mfpackuniq:
            modinfile.extend([ modelName + '.' + packname.lower() ]) #swm: append package name
            templatefile.extend([modelName + '.' + packname.lower() +'.tpl'])    #swm: append package name +.tpl
        #
        # ----- Get data for Observations
        self.updateProgressBar(start_time, time.time(), 'Get data for Observations')
        hob_count = 0
        flwob_count = 0
        # obsname=[['h1','h21','h22','h31','h32','h33'],['Drob1_1','Drob1_2','Drob1_3','Drob2_1','Drob2_2','Drob2_3']]
        #	-------------------------------SWM: BEGIN CHANGES
        nobsf = 0  # initialize to 0
        nobs=[ ]
        obsname= [[ ]]
        obsvalue= [[ ]]
        statistic = [[ ]]
        statflag = [[ ]]
        groupname = [[ ]]
        modoutfile = [[ ]]
        instructionfile = [[ ]]
        if self.chkHob_ucode.isChecked():
            hoblayer = getVectorLayerByName(self.cmbHob.currentText())
            hob_count = hoblayer.featureCount()
            nobs.extend([hob_count])
            for f in hoblayer.getFeatures():
                obsname[nobsf].append(f['OBSNAM'])   # using nobsf as index
                obsvalue[nobsf].append(f['HOBS'])
                groupname[nobsf].append('head')
                statistic[nobsf].append(f['WEIGHT'])
                statflag[nobsf].append(f['STATISTICS'])

            modoutfile[nobsf]= modelname + '.obh'    # using nobsf as index and use the correct observation  output file name here (obh)
            instructionfile[nobsf] = modelname     + '.hob_ins'  # using nobsf as index use the correct instruction file name here (hob_ins)
            nobsf = nobsf + 1

        if self.chkRvob_ucode.isChecked():
            rvoblayer = getVectorLayerByName(self.cmbRvob_ucode.currentText())
            obsname.append([ ])  # append another observation type
            obsvalue.append([ ])
            statistic.append([ ])
            statflag.append([ ])
            groupname.append([ ])
            modoutfile.append([ ])
            instructionfile.append([ ])

            nqcfb_rvob = int(rvoblayer.featureCount()) # This is over estimated, but nqcfb is only used for max array size so this will work
            # Initialize the recarray representing the attributes table
            data_rvob = np.recarray(shape=(nqcfb_rvob), dtype =[('obsnam', object), ('flwobs', float),
                                    ('weight', float), ('statistics', object) ])
            for i, fh in enumerate(rvoblayer.getFeatures()):
                data_rvob['obsnam'][i] = fh['OBSNAM']
                data_rvob['flwobs'][i] = float(fh['FLOWOBS'])
                data_rvob['weight'][i] = float(fh['WEIGHT'])
                data_rvob['statistics'][i] = fh['STATISTICS']

# Get unique values based on obsnam and store indices and use for other variables
            obsnam_rvob   = getUniqueValuesArray(data_rvob['obsnam'])
            zz, rvob_idx   =  np.unique(data_rvob['obsnam'], return_index=True)
            flwobs_rvob   = data_rvob['flwobs'][rvob_idx]
            statistic_rvob   = data_rvob['weight'][rvob_idx]
            statflag_rvob   = data_rvob['statistics'][rvob_idx]
            flwob_count = len(obsnam_rvob)
            nobs.extend([flwob_count])

# Append unique values into arrays for ucode
            obsname[nobsf][:] = obsnam_rvob  # using nobsf as index
            obsvalue[nobsf][:] = flwobs_rvob
            groupname[nobsf][:] = [('riv_flow')]*flwob_count
            statistic[nobsf][:] = statistic_rvob
            statflag[nobsf][:] = statflag_rvob

            modoutfile[nobsf]=modelname + '.obr'    # using nobsf as index and use the correct observation output file name here (obr)
            instructionfile[nobsf] = modelname     + '.rvob_ins'  # using nobsf as index and use the correct instruction file name here (rvob_ins)
            nobsf = nobsf + 1

        if self.chkDrob_ucode.isChecked():
            droblayer = getVectorLayerByName(self.cmbDrob_ucode.currentText())
            obsname.append([ ])  # append another observation type
            obsvalue.append([ ])
            statistic.append([ ])
            statflag.append([ ])
            groupname.append([ ])
            modoutfile.append([ ])
            instructionfile.append([ ])

            nqcfb_drob = int(droblayer.featureCount()) # This is over estimated, but nqcfb is only used for max array size so this will work
            # Initialize the recarray representing the attributes table
            data_drob = np.recarray(shape=(nqcfb_drob), dtype =[('obsnam', object), ('flwobs', float),
                                    ('weight', float), ('statistics', object) ])
            for i, fh in enumerate(droblayer.getFeatures()):
                data_drob['obsnam'][i] = fh['OBSNAM']
                data_drob['flwobs'][i] = float(fh['FLOWOBS'])
                data_drob['weight'][i] = float(fh['WEIGHT'])
                data_drob['statistics'][i] = fh['STATISTICS']

# Get unique values based on obsnam and store indices and use for other variables
            obsnam_drob   = getUniqueValuesArray(data_drob['obsnam'])
            zz, drob_idx   =  np.unique(data_drob['obsnam'], return_index=True)
            flwobs_drob   = data_drob['flwobs'][drob_idx]
            statistic_drob   = data_drob['weight'][drob_idx]
            statflag_drob   = data_drob['statistics'][drob_idx]
            flwob_count = len(obsnam_drob)
            nobs.extend([flwob_count])

# Append unique values into arrays for ucode
            obsname[nobsf][:] = obsnam_drob  # using nobsf as index
            obsvalue[nobsf][:] = flwobs_drob
            groupname[nobsf][:] = [('drn_flow')]*flwob_count
            statistic[nobsf][:] = statistic_drob
            statflag[nobsf][:] = statflag_drob

            modoutfile[nobsf]=modelname + '.obd'    # using nobsf as index and use the correct observation output file name here (obd)
            instructionfile[nobsf] = modelname     + '.drob_ins'  # using nobsf as index and use the correct instruction file name here (drob_ins)
            nobsf = nobsf + 1

        if self.chkGbob_ucode.isChecked():
            gboblayer = getVectorLayerByName(self.cmbGbob_ucode.currentText())
            obsname.append([ ]) # append another observation type
            obsvalue.append([ ])
            statistic.append([ ])
            statflag.append([ ])
            groupname.append([ ])
            modoutfile.append([ ])
            instructionfile.append([ ])

            nqcfb_gbob = int(gboblayer.featureCount()) # This is over estimated, but nqcfb is only used for max array size so this will work
            # Initialize the recarray representing the attributes table
            data_gbob = np.recarray(shape=(nqcfb_gbob), dtype =[('obsnam', object), ('flwobs', float),
                                    ('weight', float), ('statistics', object) ])
            for i, fh in enumerate(gboblayer.getFeatures()):
                data_gbob['obsnam'][i] = fh['OBSNAM']
                data_gbob['flwobs'][i] = float(fh['FLOWOBS'])
                data_gbob['weight'][i] = float(fh['WEIGHT'])
                data_gbob['statistics'][i] = fh['STATISTICS']

# Get unique values based on obsnam and store indices and use for other variables
            obsnam_gbob   = getUniqueValuesArray(data_gbob['obsnam'])
            zz, gbob_idx   =  np.unique(data_gbob['obsnam'], return_index=True)
            flwobs_gbob   = data_ghbb['flwobs'][gbob_idx]
            statistic_gbob   = data_ghbb['weight'][gbob_idx]
            statflag_gbob   = data_ghbb['statistics'][gbob_idx]
            flwob_count = len(obsnam_gbob)
            nobs.extend([flwob_count])

# Append unique values into arrays for ucode
            obsname[nobsf][:] = obsnam_gbob  # using nobsf as index
            obsvalue[nobsf][:] = flwobs_gbob
            groupname[nobsf][:] = [('ghb_flow')]*flwob_count
            statistic[nobsf][:] = statistic_gbob
            statflag[nobsf][:] = statflag_gbob

            modoutfile[nobsf]=modelname + '.obg'    # using nobsf as index and use the correct observation output file name here (obg)
            instructionfile[nobsf] = modelname     + '.gbob_ins'  # using nobsf as index and use the correct instruction file name here (gbob_ins)
            nobsf = nobsf + 1

        if self.chkChob_ucode.isChecked():
            choblayer = getVectorLayerByName(self.cmbChob_ucode.currentText())
            obsname.append([ ]) # append another observation type
            obsvalue.append([ ])
            statistic.append([ ])
            statflag.append([ ])
            groupname.append([ ])
            modoutfile.append([ ])
            instructionfile.append([ ])

            nqcfb_chob = int(choblayer.featureCount()) # This is over estimated, but nqcfb is only used for max array size so this will work
            # Initialize the recarray representing the attributes table
            data_chob = np.recarray(shape=(nqcfb_chob), dtype =[('obsnam', object), ('flwobs', float),
                                    ('weight', float), ('statistics', object) ])
            for i, fh in enumerate(choblayer.getFeatures()):
                data_chob['obsnam'][i] = fh['OBSNAM']
                data_chob['flwobs'][i] = float(fh['FLOWOBS'])
                data_chob['weight'][i] = float(fh['WEIGHT'])
                data_chob['statistics'][i] = fh['STATISTICS']

# Get unique values based on obsnam and store indices and use for other variables
            obsnam_chob   = getUniqueValuesArray(data_chob['obsnam'])
            zz, chob_idx   =  np.unique(data_chob['obsnam'], return_index=True)
            flwobs_chob   = data_chob['flwobs'][chob_idx]
            statistic_chob   = data_chob['weight'][chob_idx]
            statflag_chob   = data_chob['statistics'][chob_idx]
            flwob_count = len(obsnam_chob)
            nobs.extend([flwob_count])

# Append unique values into arrays for ucode
            obsname[nobsf][:] = obsnam_chob  # using nobsf as index
            obsvalue[nobsf][:] = flwobs_chob
            groupname[nobsf][:] = [('chd_flow')]*flwob_count
            statistic[nobsf][:] = statistic_chob
            statflag[nobsf][:] = statflag_chob

            modoutfile[nobsf]=modelname + '.obc'    # using nobsf as index and use the correct observation output file name here (obc)
            instructionfile[nobsf] = modelname     + '.chob_ins'  # using nobsf as index and use the correct instruction file name here (chob_ins)
            nobsf = nobsf + 1

#        QtGui.QMessageBox.information(None, 'Information', 'nobs= \n' + str(max(nobs)) )
#        QtGui.QMessageBox.information(None, 'Information', 'nobsf= \n' + str(nobsf) )

#      ----------------------------SWM: END CHANGES
        self.updateProgressBar(start_time, time.time(), 'Writing UCODE input file...')

        # Actually, tolpar, senmethod and maxchnge could be set for each parameter: set a list of a constant value, so far!
        tolpar = [tolpar for i in range(npar)]
        xchange = [xchange for i in range(npar)]
        senmethod = [senmethod for i in range(npar)]

        # --- Pass data to Ucode Writer Object and run WRITING method
##        old: ucodewriter= UcodeWriter(ml,verbose=verbose,modelname=modelname,modellengthunits=modellengthunits,modelmassunits=modelmassunits,modeltimeunits=modeltimeunits, sensitivities= sensitivities, optimize= optimize, eigenvalues= eigenvalues,
##                              tolpar=tolpar,tolsosc=tolsosc,maxiter=maxiter, maxchange=xchange, command=command, commandid=commandid,npar=npar,parname=parname,
##                              startvalue=startvalue,adjustable=adjustable,transform=transform,nobsf=nobsf,nobs=nobs,obsname=obsname,obsvalue=obsvalue,statistic=statistic,
##                              statflag=statflag,groupname=groupname,modinfile=modinfile,templatefile=templatefile,modoutfile=modoutfile,instructionfile=instructionfile )

        ucodewriter= UcodeWriter(ml,verbose=verbose,modelname=modelname,modellengthunits=modellengthunits,modelmassunits=modelmassunits,modeltimeunits=modeltimeunits, sensitivities= sensitivities, optimize= optimize, eigenvalues= eigenvalues,
                              tolpar=tolpar,senmethod=senmethod,tolsosc=tolsosc,maxiter=maxiter, maxchange=xchange, command=command, commandid=commandid,npar=npar,parname=parname,
                              startvalue=startvalue,adjustable=adjustable,transform=transform,nobsf=nobsf,nobs=nobs,obsname=obsname,obsvalue=obsvalue,statistic=statistic,
                              statflag=statflag,groupname=groupname,modinfile=modinfile,templatefile=templatefile,modoutfile=modoutfile,instructionfile=instructionfile,
                              constrain=constrain, lowerconstraint=lowerconstraint, upperconstraint=upperconstraint, perturbamt = perturbamt )

        ucodewriter.write_file()

        # Sleep for a while ... so the User knows from the ProgressBar that everything was done!
        self.closeProgressBar()
        time.sleep(3)

        QtGui.QMessageBox.information(None, 'Information', 'UCODE Input files has been written in: \n' + pathfile )


    def writeTemplate(self):
        # Define starting time for tracking the process duration,
        #        and launch the progress dialog
        start_time = time.time()
        self.launchProgressBar()
        #
        modelName = self.cmbModelName.currentText()
        #Get basic info of MODFLOW model from modelname
        # path and stress periods
        (pathfile, nsp ) = getModelInfoByName(modelName)
        # MODFLOW model FloPy object
        namefile = modelName + '.nam'
        ml = Modflow.load(namefile, version='mf2005', model_ws=pathfile)
        # nrow, ncol, nlay, nper
        (nrow, ncol, nlay, nper) = ml.get_nrow_ncol_nlay_nper()

        # Get info on parameters from GUI, parameters Table
        # Initialize the list of Parameters objects, and npar (number of parameters)
        p = []
        npar = 0

        # 3D parameters (LPF)
        if self.chk3Dparams.isChecked():
            self.updateProgressBar(start_time, time.time(), 'Getting data for LPF package')
            partable = getVectorLayerByName(self.cmbParameters.currentText())
            npar = npar + partable.featureCount()
            parname = []
            partype = []
            startvalue =[]
            adjustable = []
            transform= []
            layers = []
            # These parameters and "transform" have not affect yet, but may in the future
            lbound = 0.001
            ubound = 1000.

            for f in partable.getFeatures():
                parname = f['Name']
                partype = f['Type']
                mfpackage = 'LPF' # iac: so far only LPF is supported (June 2016)
                startvalue = float(f['StartValue'])
                transform = f['Transform']
                layer = int(f['layer'])
                if 'Yes' in f['Constr']:
                    ubound = float(f['Upper'])
                    lbound = float(f['Lower'])

                # Check if Zone Option is active
                if 'Yes' in f['UseZones']:
                    # Create the zone array from GUI
                    zonearray = np.ones((nlay, nrow, ncol), dtype=int)
                    for l in range(0,nlay):
                        ztemp = freewat.createGrid_utils.get_param_array(getVectorLayerByName(f['ZoneLayer']), fieldName = 'zone_lay_%s'%str(l + 1))
                        zonearray[l, : ,: ] = ztemp
                    # create the span dictionary, according to Zone Id
                    span = {}
                    span['idx'] = np.where(zonearray == int(f['ZoneId']))
                else:
                    idx = np.empty((nlay, nrow, ncol), dtype=np.bool)
                    idx[:,:,:] = False
                    idx[layer - 1] = True
                    span = {'idx': idx}
                # The span variable defines how the parameter spans the package
                p.append(Params(mfpackage, partype, parname, startvalue, lbound, ubound, span))

        # 2D-time varying parameters (RCH, EVT)
        if self.chk2Dparams.isChecked():
            self.updateProgressBar(start_time, time.time(), 'Getting data for RCH/LPF package')
            partable = getVectorLayerByName(self.cmb2DParameters.currentText())
            parname = []
            partype = []
            mfpackage = []
            startvalue =[]
            adjustable = []
            transform= []
            kperlist = [] # List of stress period list, for each parameter; e.g. [[1],[2,4],[1,3,4]]
            idxlist = [] # List of idx value, for each parameter (None if multiplier, zone array if Zones are used)

            # These parameters and "transform" have not affect yet, but may in the future
            lbound = 0.001
            ubound = 1000.
            #
            ip = 0
            for f in partable.getFeatures():
                name_temp = f['Name']
                sp_id = f['StressPeriod']
                if name_temp not in parname:
                    parname.append(name_temp)
                    kperlist.append([sp_id - 1])   #swm: 0 indexing in FloPy
                    partype.append(f['Type'])
                    mfpackage.append(f['Package'])
                    startvalue.append(float(f['StartValue']))
                    transform.append(f['Transform'].lower())
                    if 'Yes' in f['Constr']:
                        ubound = float(f['Upper'])
                        lbound = float(f['Lower'])
                    # Check if Zone Option is active
                    if 'Yes' in f['UseZones']:
                        # Create the zone array from GUI
                        zonearray = np.ones((nrow, ncol), dtype=int)
                        ztemp = freewat.createGrid_utils.get_param_array(getVectorLayerByName(f['ZoneLayer']), fieldName = 'zone_lay_1')
                        zonearray[ : ,: ] = ztemp
                        # create the idx key for span dictionary, according to Zone Id
                        idxtemp = np.empty((nrow, ncol), dtype=np.bool)
                        idxtemp[:,:] = False
                        for i in range(nrow):
                            for j in range(ncol):
                                if zonearray[i,j] == int(f['ZoneId']):
                                    idxtemp[i,j] = True
                        #idx = np.where(zonearray == int(f['ZoneId']))
                        idxlist.append(idxtemp)
                    else:
                        # If Zone Option is NO, use multiplier!
                        idxlist.append(None)

                    ip += 1

                else:
                    kperlist[ip-1].append(sp_id - 1)   #swm: 0 indexing in FloPy


            # Define Parameters objects and save in list p
            for jp, prm in enumerate(parname):
                npar += 1
                span = {'kpers': kperlist[jp], 'idx': idxlist[jp]}
                p.append(Params(mfpackage[jp], partype[jp], parname[jp], startvalue[jp], lbound, ubound, span))
            #

        # Write template files (here p is the list of Params objects, both 3D and 2D time-variant)
        self.updateProgressBar(start_time, time.time(), 'Writing Template File ... ')
        #
        tw = TemplateWriter(ml, p)
        tw.write_template()
        #
        # Sleep for a while ... so the User knows from the ProgressBar that everything was done!
        self.closeProgressBar()
        time.sleep(3)
        #
        QtGui.QMessageBox.information(None, 'Information', 'UCODE Template files have been written in: \n' + pathfile )


    def runUcode(self):
        # Define starting time for tracking the process duration,
        #        and launch the progress dialog
        start_time = time.time()
        self.launchProgressBar()
        #
        modelName = self.cmbModelName.currentText()
        # pop message if the value of the dictionary is empty
        dirdict = self.getProgramLocation()
        if dirdict['UCODE'] == '':
            pop_message(self.tr('UCODE executable is needed!\nPlease enter the path of such executable in the prg_locations table!'), self.tr('warning'))
            return
        else:
            ucodedir = dirdict['UCODE']
        if dirdict['MF2005'] == '':
            pop_message(self.tr('MODFLOW executable is needed! Please enter the path of such executable in the prg_locations table!'), self.tr('warning'))
            return
        else:
            modflowdir = dirdict['MF2005']

        # path and stress periods
        (pathfile, nsp ) = getModelInfoByName(modelName)

        inputfile = modelName + '_ucode.in' # self.model.model_ws + '\\' +
        ucodeoutput =  modelName # '  \\UCODE_out\\' + modelName
        #ucodecommand = inputfile + ' ' + modelName

        messageRun = 'You are attempting to run: \n' + ucodedir + '\n' + 'and output is saved under the folder \UCODE_out' + '\n' + 'Click OK to start!'
        QtGui.QMessageBox.information(None, 'Information', messageRun  )
        #
        self.updateProgressBar(start_time, time.time(), 'UCODE is running ... ')
        #
        proc = sub.Popen([ucodedir, inputfile, ucodeoutput], stdin = sub.PIPE, stdout = sub.PIPE, stderr = sub.PIPE, cwd=pathfile)
        #
        result = proc.communicate()[0]
        #
        self.updateProgressBar(start_time, time.time(), 'Writing simulation results ... ')
        # Write results on LOG file.
        logfileName  = os.path.join(pathfile, modelName + '_UCODE_log.txt')
        logfile = open(logfileName, 'w')
        logfile.write(result)
        logfile.close()

        # Sleep for a while ... so the User knows from the ProgressBar that everything was done!
        self.updateProgressBar(start_time, time.time(), 'Simulation concluded ! ')
        time.sleep(3)
        self.closeProgressBar()

        txtInfo = 'To check if simulation was successful \n \nsee details on log file: \n \n ' + logfileName + '\n \nYou can see the report of results by clicking Open Report'
        #self.displayMessage(logfileName,result)
        QtGui.QMessageBox.information(None, 'Information', txtInfo  )

##        QtGui.QMessageBox.information(None, 'Information', result  )

    def openReportMF(self):
        # Open the list file of the MODFLOW model selected
        import webbrowser
        #
        modelName = self.cmbModelName.currentText()
        # Retrieve numebr of stress periods from FlowModel
        (pathFile, nsp ) = getModelInfoByName(modelName)
        report = os.path.join(pathFile, modelName + '.list')
        # report = pathFile + "\\" + modelName + '.list'
        webbrowser.open(report)
#
    def openReportMT(self):
        import webbrowser
        # Open the list file of the MT3DMS or SEAWAT model selected
        modelName = self.cmbModelName.currentText()
        transportName = self.cmbTransportName.currentText()
        # Retrieve numebr of stress periods from FlowModel
        (pathFile, nsp ) = getModelInfoByName(modelName)
        report = os.path.join(pathFile, transportName + '.list')
        # report = pathFile + "\\" + transportName + '.list'
        webbrowser.open(report)

#
    def openReportUCODE(self):
        import webbrowser
        # Open the UCODE output file for the model selected
        modelName = self.cmbModelName.currentText()
        # Retrieve model path
        (pathFile, nsp ) = getModelInfoByName(modelName)
        report = os.path.join(pathFile, modelName + '.#uout')
        # report = pathFile + "\\" + modelName + '.#uout'
        webbrowser.open(report)
