# -*- 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, Mt3dSft, Mt3dLkt, 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
        chd_layers_flow = []
        for nametemp in layerNameList:
            if nametemp[-4:] == '_chd':
                chd_layers_flow.append(nametemp)
        self.cmbCHD.addItems(chd_layers_flow)

        chd_layers = []
        for nametemp in layerNameList:
            if '_chd' in nametemp:
                chd_layers.append(nametemp)
        self.cmbSSM_CHD.addItems(chd_layers)
        #
        # wel
        wel_layers_flow = []
        for nametemp in layerNameList:
            if nametemp[-5:] == '_well':
                wel_layers_flow.append(nametemp)
        self.cmbWEL.addItems(wel_layers_flow)

        wel_layers = []
        for nametemp in layerNameList:
            if '_wel' in nametemp:
                wel_layers.append(nametemp)
        self.cmbSSM_WEL.addItems(wel_layers)
        # mnw
        mnw_layers_flow = []
        for nametemp in layerNameList:
            if nametemp[-4:] == '_mnw':
                mnw_layers_flow.append(nametemp)
        self.cmbMNW.addItems(mnw_layers_flow)
        # riv
        riv_layers_flow = []
        for nametemp in layerNameList:
            if nametemp[-4:] == '_riv':
                riv_layers_flow.append(nametemp)
        self.cmbRIV.addItems(riv_layers_flow)

        riv_layers = []
        for nametemp in layerNameList:
            if '_riv' in nametemp:
                riv_layers.append(nametemp)
        self.cmbSSM_RIV.addItems(riv_layers)
        #lak
        lak_layers = []
        for nametemp in layerNameList:
            if 'lak_layer_group' in nametemp:
                lak_layers.append(nametemp)
        self.cmbLAK.addItems(lak_layers)
        # ghb
        ghb_layers_flow = []
        for nametemp in layerNameList:
            if nametemp[-4:] == '_ghb':
                ghb_layers_flow.append(nametemp)
        self.cmbGHB.addItems(ghb_layers_flow)

        ghb_layers = []
        for nametemp in layerNameList:
            if '_ghb' in nametemp:
                ghb_layers.append(nametemp)
        self.cmbSSM_GHB.addItems(ghb_layers)
        # sfr
        sfr_layers_flow = []
        for nametemp in layerNameList:
            if nametemp[-4:] == '_sfr':
                sfr_layers_flow.append(nametemp)
        self.cmbSFR.addItems(sfr_layers_flow)

        sfr_tables = []
        for nametemp in layerNameList:
            if nametemp[-10:] == '_sfr_table':
                sfr_tables.append(nametemp)
        self.cmbSfrTable.addItems(sfr_tables)
        # rch
        rch_layers_flow = []
        for nametemp in layerNameList:
            if nametemp[-4:] == '_rch':
                rch_layers_flow.append(nametemp)
        self.cmbRCH.addItems(rch_layers_flow)
        #
        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_flow = []
        for nametemp in layerNameList:
            if nametemp[-4:] == '_drn':
                drn_layers_flow.append(nametemp)
        self.cmbDRN.addItems(drn_layers_flow)
        #
        # evt
        evt_layers_flow = []
        for nametemp in layerNameList:
            if nametemp[-4:] == '_evt':
                evt_layers_flow.append(nametemp)
        self.cmbEVT.addItems(evt_layers_flow)
        #
        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_flow = []
        surf_layers_flow = []
        for nametemp in layerNameList:
            if nametemp[-4:] == '_uzf':
                uzf_layers_flow.append(nametemp)
            if nametemp[-4:] == '_sml':
                surf_layers_flow.append(nametemp)
        self.cmbUZF.addItems(uzf_layers_flow)
        self.cmbSurf.addItems(surf_layers_flow)
        # swi2
        swi_layers_flow = []
        for nametemp in layerNameList:
            if nametemp[-4:] == '_swi':
                swi_layers_flow.append(nametemp)
        self.cmbSWI.addItems(swi_layers_flow)

        # Observations
        # hob
        hob_layers = []
        for nametemp in layerNameList:
            if nametemp[-4:] == '_hob':
                hob_layers.append(nametemp)
        self.cmbHOB.addItems(hob_layers)

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

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

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

        # chob
        chob_layers = []
        for nametemp in layerNameList:
            if nametemp[-5:] == '_chob':
                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)

        # sft
        #sft_layers is a list which will be filled with all the _sft layers available in the Legend
        sft_layers = []
        #sft_tables is a list which will be filled with all the sft tables available in the Legend.
        #The first element of the list sft_tables is the string 'No stream boundary conditions'
        sft_tables=['No stream boundary conditions']
        #nametemp iterates over the elements of the list layerNameList.
        #The list layerNameList is filled with the names of all the vector layers available
        #in the Legend (at the beginning of this function, layerNameList is created using the
        #getVectorLayerNames function of the freewat_utils.py module; such function returns a list)
        for nametemp in layerNameList:
            #The elements of the list layerNameList are strings.
            #If the last 4 characters of a string correspond to '_sft'...
            if nametemp[-4:] == '_sft':
                #...then the list sft_layers is updated
                sft_layers.append(nametemp)
            #If the elements of the list layerNameList contain '_sft_table_'...
            if '_sft_table_' in nametemp:
                #...then the list sft_tables is updated
                sft_tables.append(nametemp)
        # fill the cmbSFT combobox with only the _sft layers
        self.cmbSFT.addItems(sft_layers)
        # fill the cmbSFT_bc combobox with only the sft tables
        self.cmbSFT_bc.addItems(sft_tables)

        #The combo box cmbISFSOLV is updated with the item 'Finite-Difference'
        self.cmbISFSOLV.addItems(['Finite-Difference'])

        #The combo box cmbWIMP is updated with the items 'Crank-Nicolson scheme', 'Explicit scheme', 'Fully implicit scheme'
        self.wimp_list = ['Crank-Nicolson scheme', 'Explicit scheme', 'Fully implicit scheme']
        self.cmbWIMP.addItems(self.wimp_list)

        #The combo box cmbWUPS is updated with the items 'Central-in-space scheme', 'Upstream scheme'
        self.wups_list = ['Central-in-space scheme', 'Upstream scheme']
        self.cmbWUPS.addItems(self.wups_list)

        #lkt
        #lake_start_conc is a list which will be filled with the lake_starting_conc_ table available in the Legend
        lake_start_conc=[]
        #lake_bc is a list which will be filled with the lake_table_bc_ table available in the Legend
        lake_bc=[]
        #nametemp iterates over the elements of the list layerNameList.
        #The list layerNameList is filled with the names of all the vector layers available
        #in the Legend (at the beginning of this function, layerNameList is created using the
        #getVectorLayerNames function of the freewat_utils.py module; such function returns a list)
        for nametemp in layerNameList:
            #If the elements of the list layerNameList contain 'lake_starting_conc_'...
            if 'lake_starting_conc_' in nametemp:
                #...then the list lake_start_conc is updated
                lake_start_conc.append(nametemp)
            #If the elements of the list layerNameList contain 'lake_table_bc_'...
            if 'lake_table_bc_' in nametemp:
                #...then the list lake_bc is updated
                lake_bc.append(nametemp)
        # fill the cmbLKT_start_conc combobox with only the lake_starting_conc_ tables
        self.cmbLKT_start_conc.addItems(lake_start_conc)
        # fill the cmbLKT_bc combobox with only the lake_table_bc_ tables
        self.cmbLKT_bc.addItems(lake_bc)


        #
        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_flow)

        # 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 = str(self.cmbMNW.currentText()) + '_table'
            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())
        #
        if self.boxSWI.isChecked():
            swilayer = getVectorLayerByName(self.cmbSWI.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 = []

            for f in mnwlayer.getFeatures():
                wellidlist.append(f['WELLID'])

            mnwmax = len(wellidlist)

            itmp = np.zeros(nper, dtype='int32')
            itmp[:] = mnwmax

            ##########Build node_data recarray##########

            #The number of rows in the table mnwtable is counted (this number is stored in the variable rows_mnwtable)
            rows_mnwtable=0
            for f in mnwtable.getFeatures():
                rows_mnwtable=rows_mnwtable+1

            #Initialize node_data recarray and its columns
            node_data = np.recarray((rows_mnwtable,), dtype=[('k', '<i8'), ('i', '<i8'), ('j', '<i8'), ('wellid', 'O'), ('losstype', 'O'), ('pumploc', '<i8'), ('qlimit', '<i8'), ('ppflag', '<i8'), ('pumpcap', '<i8'), ('rw', '<f4'), ('B', '<f4'), ('C', '<f4'), ('P', '<f4')])

            for i, f in enumerate(mnwtable.getFeatures()):
                    node_data['k'][i]  =  f['layer'] - 1       #get the layer number (decreased by 1)
                    for row in mnwlayer.getFeatures():
                        if row['WELLID']==f['well_id']:
                            node_data['i'][i] = row['ROW'] - 1 #get the row number (decreased by 1)
                            node_data['j'][i] = row['COL'] - 1 #get the column number (decreased by 1)
                    node_data['wellid'][i] = f['well_id']      #get the well ID
                    node_data['losstype'][i] = 'general'
                    node_data['pumploc'][i] = 0
                    node_data['qlimit'][i] = 0
                    node_data['ppflag'][i] = 0
                    node_data['pumpcap'][i] = 0
                    node_data['rw'][i] = f['Rw']               #get the well radius
                    node_data['B'][i] = f['B']                 #get the B parameter for the well
                    node_data['C'][i] = f['C']                 #get the C parameter for the well
                    node_data['P'][i] = f['P']                 #get the P parameter for the well

            ##########Build stress_period_data dictionary##########

            # Initialize stress_period_data dictionary
            stress_period_data = {}

            for j in range(nper):
                # Initialize rec_array for current stress period
                mnw_data = np.recarray((mnwmax,),  dtype= [('wellid', 'O'), ('qdes', '<f4')])
                for i, f in enumerate(mnwlayer.getFeatures()):
                    # get data only for current stress period
                    mnw_data['wellid'][i] = f['WELLID']        #get the well ID
                    mnw_data['qdes'][i] = f['Qw_' + str(j+1)]    #get the well pumping rate

                stress_period_data[j] = mnw_data

        # ----
        # 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
        # ----
        # Retrieve data for SWI:
        if self.boxSWI.isChecked():
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for SWI2')
            #
            zeta_array = np.zeros((nlay, nrow, ncol))
            iso = np.zeros((nlay, nrow, ncol), dtype=np.int)
            ssz = np.zeros((nlay, nrow, ncol))

            # Select IF swilayer is in DB (Spatialite) or not:
            if checkIfTableExists(self.dbpath,swilayer.name()) :
                # This means that SWI MDO is in DB... so Use the Fast Method
                for i in range(1,nlay+1):
                    z_record  = getTableArrayData(self.dbpath, swilayer.name() , col = 'zeta_%s'%(i) , nrow = nrow, ncol = ncol)
                    zeta_array[i-1,:,:] = z_record
                    ssz_record = getTableArrayData(self.dbpath, swilayer.name() , col = 'ssz_%s'%(i), nrow = nrow, ncol = ncol)
                    ssz[i-1,:,:] = ssz_record
                    iso_record = getTableArrayData(self.dbpath, swilayer.name() , col = 'isource_%s'%(i), nrow = nrow, ncol = ncol)
                    iso[i-1,:,:] = iso_record
            else:
                # This means that SWI MDO is NOT in DB... so Use the Old Slower Method
                for i in range(1,nlay+1):
                    for f in swilayer.getFeatures():
                        nr = f['row'] - 1
                        nc = f['col'] - 1
                        zeta_array[i-1,nr,nc] = f['zeta_%s'%(i)]
                        ssz[i-1,nr,nc] = f['ssz_%s'%(i)]
                        iso[i-1,nr,nc] = f['isource_%s'%(i)]
            # z list should be long as NSRF, which here is 1 for default!
            z_list = [zeta_array]

            # Retrieve SWI data from GUI:
            fresh_density = float(self.spinFreshDens.value())
            salt_density = float(self.spinSaltDens.value())
            nu = [fresh_density, salt_density]
            toeslope = float(self.spinToeSlope.value())
            tipslope = float(self.spinTipSlope.value())
        # ----
        ### ---------- 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, ipakcb=unit_cbcfile, stress_period_data = layer_row_column_Q )

        # MNW2 (MNW2 Package)
        if self.ckMnw.isChecked():
            mnw = ModflowMnw2(ml, mnwmax= mnwmax, ipakcb=unit_cbcfile, mnwprnt= 0, node_data= node_data, mnw= None,
                              stress_period_data= stress_period_data, itmp= itmp, extension= 'mnw2', unitnumber= 34, gwt= False)

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

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

        ### ---------- Add lak to model ----------###
        if self.ckLAK.isChecked():
            ModflowLak(ml, nlakes=lak_nlakes, ipakcb=unit_cbcfile, 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, ipakcb=unit_cbcfile, stress_period_data = drn_dict )

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

        # EVT (EVAPORATION package)
        if self.ckEvt.isChecked():
            # default value for ievt = 1
            evt = ModflowEvt(ml, nevtop = evt_option, ipakcb=unit_cbcfile, 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=unit_cbcfile, 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, ipakcb=unit_cbcfile, 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, ipakcb=unit_cbcfile, istcb2=0, 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, ipakcb=unit_cbcfile, istcb2=0, 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 = int(self.txtOuter.text())
        iter1 = int(self.txtInner.text())
        hclose = float(self.txtHclose.text())
        rclose = float(self.txtRclose.text())
        relax = float(self.txtRelax.text())

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

        iprpcg = int(self.txtIprpcg.text())
        damp = float(self.txtDamp.text())
        # Retrieve also MUTPCG from UI
        mutpcgText = str(self.cmbMUTPCG.currentText()).lower()
        if 'maximum' in mutpcgText:
            mutpcg = 0
        elif 'iterations' in mutpcgText:
            mutpcg = 1
        elif 'printing' in mutpcgText:
            mutpcg = 2
        else:
            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)

        #
        if self.boxSWI.isChecked():
            # SWI2 package (here also PCG info before loaded are used and passed to the package
            """ LIMITATIONS:
            NSRF = 1 (always 2 zones only: fresh water and salt water)
            ISTRAT = 1 (density is constant between surfaces.)
            NO "OPTIONS" (not ADAPTIVE method)
            NSOLVER = 2 (PCG)
            NBPOL = 2
            DAMPT = 1.0
            """
            solver2parameters ={'mxiter': mxiter, 'iter1': iter1, 'npcond': npcond,
                                'zclose': hclose, 'rclose': rclose, 'relax': relax,
                                'nbpol': 2, 'damp': damp, 'dampt': 1.0}

            swi = ModflowSwi2(ml, nsrf=1, istrat=1, toeslope=toeslope, tipslope=tipslope, nu=nu,
                                    zeta=z_list, ssz=ssz, isource=iso, nsolver=2,
                                    iswizt=55, solver2params=solver2parameters)

        # 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 SFT (Stream Flow Transport)
        sftActive = False
        #If the check box boxSFT is checked...
        if self.boxSFT.isChecked():
            #Update the progress bar while the .sft file is being written
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for SFT')

            # Check that SFR is present in flow model:
            #package_list is a list with all the packages activated in the flow model ml
            #(ml is a Modflow object created recalling above the class Modflow
            #present in the module mf.py within flopy/flopy/modflow/;
            #the method get_package_list reads the .nam file; the flow model
            #taken into account is the one whose name is read from the combo box
            #cmbModelName at the top of the Rum Model window).
            #ATTENTION: this check is necessary because to run SFT, the SFR package is needed!!!
            package_list = ml.get_package_list()
            #If the SFR package is not in the list package_list...
            if 'SFR' not in package_list:
                # message error:
                errorMsg = '''
                You are attempting to simulate solute transport in '
                surface streams, but you didn't apply SFR package
                in the flow model !!!
                You need to run MODFLOW once again
                activating SFR.
                Please, do it, and come back here later!!  '''

                #The message box which appears is entitled 'No SFR package found '
                QtGui.QMessageBox.warning(None, 'No SFR package found ', errorMsg )
                #If this error occurs, nothing is done
                return

            #Check if MODFLOW-NWT executable is present.
            #ATTENTION: MODFLOW-NWT is needed to run MT3D-USGS afterwards!!!
            if mfnwtdir == '' :
                # message error:
                errorMsg = '''
                You are attempting to simulate solute transport in
                surface streams, 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  '''

                #The message box which appears is entitled 'check MF-NEWTON Executable in Program Locations '
                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 !!

            #Check if MT3D-USGS executable is present.
            #ATTENTION: MT3D-USGS is needed to run the SFT package!!!
            if mtusgsdir == '' :
                # message error:
                errorMsg = ''' You are attempting to simulate solute transport in '
                               surface streams, but it seems that you didn't enter
                               MT3D-USGS executable in Program Locations!
                               Please, do it, and come back here later!!  '''

                #The message box which appears is entitled 'No MT3D-USGS executable in Program Locations '
                QtGui.QMessageBox.warning(None, 'No MT3D-USGS executable in Program Locations ', errorMsg )
                #If this error occurs, nothing is done
                return

            # If yes, proceed...
            #If the SFR package is in the list package_list...
            sftActive = True

            #isfsolv_opt is the content of the combo box cmbISFSOLV
            isfsolv_opt = self.cmbISFSOLV.currentText()
            #If the option 'Finite-Difference' is chosen...
            if isfsolv_opt == 'Finite-Difference' :
                #...then isfsolv_param = 1
                isfsolv_param = 1

            #wimp_opt is the content of the combo box cmbWIMP
            wimp_opt = self.cmbWIMP.currentText()
            #If the option 'Crank-Nicolson scheme' is chosen...
            if wimp_opt == 'Crank-Nicolson scheme' :
                #...then wimp_param = 0.5
                wimp_param = 0.5
            #If the option 'Explicit scheme' is chosen...
            elif wimp_opt == 'Explicit scheme' :
                #...then wimp_param = 0.0
                wimp_param = 0.0
            #If the option 'Fully implicit scheme' is chosen...
            elif wimp_opt == 'Fully implicit scheme' :
                #...then wimp_param = 1.0
                wimp_param = 1.0

            #wups_opt is the content of the combo box cmbWUPS
            wups_opt = self.cmbWUPS.currentText()
            #If the option 'Central-in-space scheme' is chosen...
            if wups_opt == 'Central-in-space scheme' :
                #...then wups_param = 0.0
                wups_param = 0.0
            #If the option 'Upstream scheme' is chosen...
            elif wups_opt == 'Upstream scheme' :
                #...then wups_param = 1.0
                wups_param = 1.0

            #cclosesf_param is the content of the line edit txtCCLOSESF
            #translated in a float number
            cclosesf_param = float(self.txtCCLOSESF.text())

            #mxitersf_param is the content of the line edit txtMXITERSF
            #translated in an integer number
            mxitersf_param = int(self.txtMXITERSF.text())

            #crntsf_param is the content of the line edit txtCRNTSF
            #translated in a float number
            crntsf_param = float(self.txtCRNTSF.text())

            # Get info from MDO
            #sftlayer is the _sft MDO whose name is read from the combo box cmbSFT
            sftlayer = getVectorLayerByName(self.cmbSFT.currentText())

            #sfttable is the sft table whose name is read from the combo box cmbSFT_bc
            sfttable = getVectorLayerByName(self.cmbSFT_bc.currentText())

            #modelName is the name of the flow model read from the combo box cmbModelName
            modelName = self.cmbModelName.currentText()
            #The tuple (pathfile, nsp ) is made of two objects:
            # -pathfile is the path of the working directiry for a specific flow model
            # -nsp is the number of stress period for a specific flow model
            #(getModelInfoByName is a function of the freewat_utils.py module; such function
            #recalls the getModelsInfoLists function to retrieve the names of all the flow
            #models available in the legend and the corresponding paths, then it retrieves
            #the path of the specific flow model selected in the cmbFlowName combo box and
            #iterates over the rows of the timetable_modelName table to get the number of
            #stress periods; it returns a tuple)
            (pathfile, nsp ) = getModelInfoByName(modelName)

            #If stream boundary conditions exist...
            if self.cmbSFT_bc.currentText() != 'No stream boundary conditions':

                reach_ID=[]
                coldsf=[]
                dispsf=[]
                #f iterates over the rows of the Attribute Table of the sftlayer layer
                for f in sftlayer.getFeatures():
                    #The list reach_ID is updated with values read from the field 'reach_ID'
                    reach_ID.append(f['reach_ID'])
                    #The list coldsf is updated with values read from the field 'coldsf_spec_1'
                    coldsf.append(f['coldsf_spec_1'])
                    #The list dispsf is updated with values read from the field 'dispsf_spec_1'
                    dispsf.append(f['dispsf_spec_1'])
                #The zip method allows to get a tuple, whose elements are lists of three
                #elements each ([reach_ID,coldsf,dispsf],[reach_ID,coldsf,dispsf],...,[reach_ID,coldsf,dispsf]).
                #The sorted method allows then to sort the tuple. The sorting is performed
                #based on the first elements of each list (i.e., based on the values of the
                #'reach_ID' field, which are unique values)
                a=sorted(zip(reach_ID,coldsf,dispsf))

                coldsf_params=[]
                dispsf_params=[]
                #j iterates over the length of the tuple a
                for j in range(len(a)):
                    #The list coldsf_params is updated with the second elements of each list
                    #in the tuple a => the list coldsf_params contains coldsf values of the
                    #sft cells ordered from reach_ID=1 to reach_ID=n
                    coldsf_params.append(a[j][1])
                    #The list dispsf_params is updated with the third elements of each list
                    #in the tuple a => the list dispsf_params contains dispsf values of the
                    #sft cells ordered from reach_ID=1 to reach_ID=n
                    dispsf_params.append(a[j][2])

                bc_params={}
                bc=[]
                bc_list=[]
                #n iterates over the number of stress periods
                for n in range(nsp):
                    #f iterates over the rows of the sft table
                    for f in sfttable.getFeatures():
                        #The list bc_list is updated with values read from the field 'reach_ID'
                        #ATTENTION: f['reach_ID'] is decreased by 1 because in the file .sft
                        #it is not well written, even if the dictionary bc_params is correct!!!
                        bc_list.append(f['reach_ID']-1)
                        #The list bc_list is updated with values read from the field 'bc_type'
                        #If the content of the field 'bc_type' is 'headwater'...
                        if f['bc_type']=='headwater':
                            #...then the list bc_list is updated with 0
                            bc_list.append(0)
                        #If the content of the field 'bc_type' is 'precipitation'...
                        elif f['bc_type']=='precipitation':
                            #...then the list bc_list is updated with 1
                            bc_list.append(1)
                        #If the content of the field 'bc_type' is 'runoff'...
                        elif f['bc_type']=='runoff':
                            #...then the list bc_list is updated with 2
                            bc_list.append(2)
                        #If the content of the field 'bc_type' is 'constant concentration'...
                        elif f['bc_type']=='constant concentration':
                            #...then the list bc_list is updated with 3
                            bc_list.append(3)
                        #If the content of the field 'bc_type' is 'pumping'...
                        elif f['bc_type']=='pumping':
                            #...then the list bc_list is updated with 4
                            bc_list.append(4)
                        #If the content of the field 'bc_type' is 'evaporation'...
                        elif f['bc_type']=='evaporation':
                            #...then the list bc_list is updated with 5
                            bc_list.append(5)
                        #If a different string is read...
                        else:
                            #...then a warning message appears and nothing is done
                            pop_message(self.tr('Invalid bc_type in the sft table!\nPlease, check the spelling of bc_type'), self.tr('warning'))
                            return
                        #The list bc_list is updated with values read from the field
                        #'conc_sp_n_spec_1' (where n is the n-th stress period)
                        #(the list bc_list is something like [reach,bc_type,conc]
                        bc_list.append(f['conc_sp_' + str(n+1) + '_spec_1'])
                        #The list bc is updated with the list bc_list
                        #(the list bc is something like [[reach1,bc_type1,conc1],[reach2,bc_type2,conc2],...])
                        bc.append(bc_list)
                        #The list bc_list is emptied
                        bc_list=[]
                    #Update the dictionary bc_params (the key is the number of the current
                    #stress period; as required by the flopy, the argument is the list bc)
                    bc_params[n]=bc
                    #The list bc is emptied
                    bc=[]

                #If there are more than one species...
                if nspec > 1:

                    reach_ID=[]
                    coldsf=[]
                    dispsf=[]

                    sft_multiComponent={}

                    #i iterates over the number of species, starting from species #2
                    for i in range(2,nspec+1):
                        #f iterates over the rows of the Attribute Table of the sftlayer layer
                        for f in sftlayer.getFeatures():
                            #The list reach_ID is updated with values read from the field 'reach_ID'
                            reach_ID.append(f['reach_ID'])
                            #The list coldsf is updated with values read from the field 'coldsf_spec_i'
                            coldsf.append(f['coldsf_spec_' + str(i)])
                            #The list dispsf is updated with values read from the field 'dispsf_spec_i'
                            dispsf.append(f['dispsf_spec_' + str(i)])
                        #The zip method allows to get a tuple, whose elements are lists of three
                        #elements each ([reach_ID,coldsf,dispsf],[reach_ID,coldsf,dispsf],...,[reach_ID,coldsf,dispsf]).
                        #The sorted method allows then to sort the tuple. The sorting is performed
                        #based on the first elements of each list (i.e., based on the values of the
                        #'reach_ID' field, which are unique values)
                        a=sorted(zip(reach_ID,coldsf,dispsf))

                        coldsf_new=[]
                        dispsf_new=[]
                        #j iterates over the length of the tuple a
                        for j in range(len(a)):
                            #The list coldsf_new is updated with the second elements of each list
                            #in the tuple a => the list coldsf_new contains coldsf values of the
                            #sft cells ordered from reach_ID=1 to reach_ID=n
                            coldsf_new.append(a[j][1])
                            #The list dispsf_new is updated with the second elements of each list
                            #in the tuple a => the list dispsf_new contains coldsf values of the
                            #sft cells ordered from reach_ID=1 to reach_ID=n
                            dispsf_new.append(a[j][2])

                        #The dictionary sft_multiComponent is updated with pairs (key:value).
                        #The keys are strings coldsfi and dispsfi (where i is the number of the i-th species).
                        #The values are the lists coldsf_new and dispsf_new (i.e., the coldsf and dispsf values
                        #for the i-th species, sorted by reach_ID).
                        #ATTENTION: the elements of the disctinary sft_multiComponent might not be sorted
                        #by keys! What is important is the name assigned to each key (such name is the same
                        #as that indicated in the module mtsft.py of the flopy)
                        sft_multiComponent['coldsf' + str(i)]=coldsf_new
                        sft_multiComponent['dispsf' + str(i)]=dispsf_new

                bc_multiComponent={}
                bc=[]
                bc_list=[]
                #n iterates over the number of stress periods
                for n in range(nsp):
                    #f iterates over the rows of the sft table
                    for f in sfttable.getFeatures():
                        #The list bc_list is updated with values read from the field 'reach_ID'
                        #ATTENTION: f['reach_ID'] is decreased by 1 because in the file .sft
                        #it is not well written, even if the dictionary bc_multiComponent is correct!!!
                        bc_list.append(f['reach_ID']-1)
                        #The list bc_list is updated with values read from the field 'bc_type'
                        #If the content of the field 'bc_type' is 'headwater'...
                        if f['bc_type']=='headwater':
                            #...then the list bc_list is updated with 0
                            bc_list.append(0)
                        #If the content of the field 'bc_type' is 'precipitation'...
                        elif f['bc_type']=='precipitation':
                            #...then the list bc_list is updated with 1
                            bc_list.append(1)
                        #If the content of the field 'bc_type' is 'runoff'...
                        elif f['bc_type']=='runoff':
                            #...then the list bc_list is updated with 2
                            bc_list.append(2)
                        #If the content of the field 'bc_type' is 'constant concentration'...
                        elif f['bc_type']=='constant concentration':
                            #...then the list bc_list is updated with 3
                            bc_list.append(3)
                        #If the content of the field 'bc_type' is 'pumping'...
                        elif f['bc_type']=='pumping':
                            #...then the list bc_list is updated with 4
                            bc_list.append(4)
                        #If the content of the field 'bc_type' is 'evaporation'...
                        elif f['bc_type']=='evaporation':
                            #...then the list bc_list is updated with 5
                            bc_list.append(5)
                        #If a different string is read...
                        else:
                            #...then a warning message appears and nothing is done
                            pop_message(self.tr('Invalid bc_type in the sft table!\nPlease, check the spelling of bc_type'), self.tr('warning'))
                            return
                        for i in range(1,nspec+1):
                            #The list bc_list is updated with values read from the field
                            #'conc_sp_n_spec_m' (where n is the n-th stress period and
                            #m is the m.th species)
                            #(the list bc_list is something like [reach,bc_type,conc_species1,conc_species2,...]
                            bc_list.append(f['conc_sp_' + str(n+1) + '_spec_' + str(i)])
                        #The list bc is updated with the list bc_list
                        #(the list bc is something like [[reach1,bc_type1,conc1_species1,conc1_species2,...],[reach2,bc_type2,conc2_species1,conc2_species2,...],...])
                        bc.append(bc_list)
                        #The list bc_list is emptied
                        bc_list=[]
                    #Update the dictionary bc_multiComponent (the key is the number of the current
                    #stress period; as required by the flopy, the argument is the list bc)
                    bc_multiComponent[n]=bc
                    #The list bc is emptied
                    bc=[]

                # build the SFT package:
                #(mt is an Mt3d object created recalling above the class Mt3dmsUSGS
                #present in the module mt.py within flopy/flopy/mt3d/.
                #ATTENTION: in the  module mt.py within flopy/flopy/mt3d/ the class
                #is actually called Mt3dms, while in the module flopyaddon/mtusgs.py
                #it is called Mt3dmsUSGS).
                #If there are more than one species...
                if nspec > 1:
                    #...coldsf and dispsf for the first species are read from the fields 'coldsf_spec_1'
                    #and 'dispsf_spec_1' (see the lists coldsf_params and dispsf_params above);
                    #for species #2 through #n, coldsf and dispsf are read from the dictionary sft_multiComponent.
                    #For all species, sf_stress_period_data is read from the dictionary bc_multiComponent
                    sft = Mt3dSft(mt, nsfinit=max(reach_ID), mxsfbc=max(reach_ID), icbcsf=0, ioutobs=0,
                    ietsfr=0, isfsolv=isfsolv_param, wimp=wimp_param, wups=wups_param, cclosesf=cclosesf_param,
                    mxitersf=mxitersf_param, crntsf=crntsf_param, iprtxmd=0, coldsf=coldsf_params, dispsf=dispsf_params,
                    nobssf=0, obs_sf=None, sf_stress_period_data=bc_multiComponent,
                    unitnumber=19, **sft_multiComponent)
                #If there is just one species...
                else:
                    #...the dictionaries sft_multiComponent and bc_multiComponent are no longer needed
                    sft = Mt3dSft(mt, nsfinit=max(reach_ID), mxsfbc=max(reach_ID), icbcsf=0, ioutobs=0,
                    ietsfr=0, isfsolv=isfsolv_param, wimp=wimp_param, wups=wups_param, cclosesf=cclosesf_param,
                    mxitersf=mxitersf_param, crntsf=crntsf_param, iprtxmd=0, coldsf=coldsf_params, dispsf=dispsf_params,
                    nobssf=0, obs_sf=None, sf_stress_period_data=bc_params,
                    unitnumber=19)

            #If no stream boundary conditions exist...
            else:

                reach_ID=[]
                coldsf=[]
                dispsf=[]
                #f iterates over the rows of the Attribute Table of the sftlayer layer
                for f in sftlayer.getFeatures():
                    #The list reach_ID is updated with values read from the field 'reach_ID'
                    reach_ID.append(f['reach_ID'])
                    #The list coldsf is updated with values read from the field 'coldsf_spec_1'
                    coldsf.append(f['coldsf_spec_1'])
                    #The list dispsf is updated with values read from the field 'dispsf_spec_1'
                    dispsf.append(f['dispsf_spec_1'])
                #The zip method allows to get a tuple, whose elements are lists of three
                #elements each ([reach_ID,coldsf,dispsf],[reach_ID,coldsf,dispsf],...,[reach_ID,coldsf,dispsf]).
                #The sorted method allows then to sort the tuple. The sorting is performed
                #based on the first elements of each list (i.e., based on the values of the
                #'reach_ID' field, which are unique values)
                a=sorted(zip(reach_ID,coldsf,dispsf))

                coldsf_params=[]
                dispsf_params=[]
                #j iterates over the length of the tuple a
                for j in range(len(a)):
                    #The list coldsf_params is updated with the second elements of each list
                    #in the tuple a => the list coldsf_params contains coldsf values of the
                    #sft cells ordered from reach_ID=1 to reach_ID=n
                    coldsf_params.append(a[j][1])
                    #The list dispsf_params is updated with the third elements of each list
                    #in the tuple a => the list dispsf_params contains dispsf values of the
                    #sft cells ordered from reach_ID=1 to reach_ID=n
                    dispsf_params.append(a[j][2])

                #If there are more than one species...
                if nspec > 1:

                    reach_ID=[]
                    coldsf=[]
                    dispsf=[]

                    sft_multiComponent={}

                    #i iterates over the number of species, starting from species #2
                    for i in range(2,nspec+1):
                        #f iterates over the rows of the Attribute Table of the sftlayer layer
                        for f in sftlayer.getFeatures():
                            #The list reach_ID is updated with values read from the field 'reach_ID'
                            reach_ID.append(f['reach_ID'])
                            #The list coldsf is updated with values read from the field 'coldsf_spec_i'
                            coldsf.append(f['coldsf_spec_' + str(i)])
                            #The list dispsf is updated with values read from the field 'dispsf_spec_i'
                            dispsf.append(f['dispsf_spec_' + str(i)])
                        #The zip method allows to get a tuple, whose elements are lists of three
                        #elements each ([reach_ID,coldsf,dispsf],[reach_ID,coldsf,dispsf],...,[reach_ID,coldsf,dispsf]).
                        #The sorted method allows then to sort the tuple. The sorting is performed
                        #based on the first elements of each list (i.e., based on the values of the
                        #'reach_ID' field, which are unique values)
                        a=sorted(zip(reach_ID,coldsf,dispsf))

                        coldsf_new=[]
                        dispsf_new=[]
                        #j iterates over the length of the tuple a
                        for j in range(len(a)):
                            #The list coldsf_new is updated with the second elements of each list
                            #in the tuple a => the list coldsf_new contains coldsf values of the
                            #sft cells ordered from reach_ID=1 to reach_ID=n
                            coldsf_new.append(a[j][1])
                            #The list dispsf_new is updated with the second elements of each list
                            #in the tuple a => the list dispsf_new contains coldsf values of the
                            #sft cells ordered from reach_ID=1 to reach_ID=n
                            dispsf_new.append(a[j][2])

                        #The dictionary sft_multiComponent is updated with pairs (key:value).
                        #The keys are strings coldsfi and dispsfi (where i is the number of the i-th species).
                        #The values are the lists coldsf_new and dispsf_new (i.e., the coldsf and dispsf values
                        #for the i-th species, sorted by reach_ID).
                        #ATTENTION: the elements of the disctinary sft_multiComponent might not be sorted
                        #by keys! What is important is the name assigned to each key (such name is the same
                        #as that indicated in the module mtsft.py of the flopy)
                        sft_multiComponent['coldsf' + str(i)]=coldsf_new
                        sft_multiComponent['dispsf' + str(i)]=dispsf_new

                # build the SFT package:
                #(mt is an Mt3d object created recalling above the class Mt3dmsUSGS
                #present in the module mt.py within flopy/flopy/mt3d/.
                #ATTENTION: in the  module mt.py within flopy/flopy/mt3d/ the class
                #is actually called Mt3dms, while in the module flopyaddon/mtusgs.py
                #it is called Mt3dmsUSGS).
                #If there are more than one species...
                if nspec > 1:
                    #...coldsf and dispsf for the first species are read from the fields 'coldsf_spec_1'
                    #and 'dispsf_spec_1' (see the lists coldsf_params and dispsf_params above);
                    #for species #2 through #n, coldsf and dispsf are read from the dictionary sft_multiComponent.
                    #For all species, sf_stress_period_data is read from the dictionary bc_multiComponent
                    sft = Mt3dSft(mt, nsfinit=max(reach_ID), mxsfbc=max(reach_ID), icbcsf=0, ioutobs=0,
                    ietsfr=0, isfsolv=isfsolv_param, wimp=wimp_param, wups=wups_param, cclosesf=cclosesf_param,
                    mxitersf=mxitersf_param, crntsf=crntsf_param, iprtxmd=0, coldsf=coldsf_params, dispsf=dispsf_params,
                    nobssf=0, obs_sf=None, sf_stress_period_data=None,
                    unitnumber=19, **sft_multiComponent)
                #If there is just one species...
                else:
                    #...the dictionaries sft_multiComponent and bc_multiComponent are no longer needed
                    sft = Mt3dSft(mt, nsfinit=max(reach_ID), mxsfbc=max(reach_ID), icbcsf=0, ioutobs=0,
                    ietsfr=0, isfsolv=isfsolv_param, wimp=wimp_param, wups=wups_param, cclosesf=cclosesf_param,
                    mxitersf=mxitersf_param, crntsf=crntsf_param, iprtxmd=0, coldsf=coldsf_params, dispsf=dispsf_params,
                    nobssf=0, obs_sf=None, sf_stress_period_data=None,
                    unitnumber=19)

        # Add LKT (Lake Transport)
        lktActive = False
        #If the check box chkLkt is checked...
        if self.chkLkt.isChecked():
            #Update the progress bar while the .lkt file is being written
            self.updateProgressBar(start_time, time.time(), string_info = 'Getting data for LKT')

            # Check that LAK is present in flow model:
            #package_list is a list with all the packages activated in the flow model ml
            #(ml is a Modflow object created recalling above the class Modflow
            #present in the module mf.py within flopy/flopy/modflow/;
            #the method get_package_list reads the .nam file; the flow model
            #taken into account is the one whose name is read from the combo box
            #cmbModelName at the top of the Rum Model window).
            #ATTENTION: this check is necessary because to run LKT, the LAK package is needed!!!
            package_list = ml.get_package_list()
            #If the LAK package is not in the list package_list...
            if 'LAK' not in package_list:
                # message error:
                errorMsg = '''
                You are attempting to simulate solute transport in '
                lakes, but you didn't apply LAK package
                in the flow model !!!
                You need to run MODFLOW once again
                activating LAK.
                Please, do it, and come back here later!!  '''

                #The message box which appears is entitled 'No LAK package found '
                QtGui.QMessageBox.warning(None, 'No LAK package found ', errorMsg )
                #If this error occurs, nothing is done
                return

            #Check if MODFLOW-NWT executable is present.
            #ATTENTION: MODFLOW-NWT is needed to run MT3D-USGS afterwards!!!
            if mfnwtdir == '' :
                # message error:
                errorMsg = '''
                You are attempting to simulate solute transport in
                lakes, 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  '''

                #The message box which appears is entitled 'check MF-NEWTON Executable in Program Locations '
                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 !!

            #Check if MT3D-USGS executable is present.
            #ATTENTION: MT3D-USGS is needed to run the LKT package!!!
            if mtusgsdir == '' :
                # message error:
                errorMsg = ''' You are attempting to simulate solute transport in '
                               lakes, but it seems that you didn't enter
                               MT3D-USGS executable in Program Locations!
                               Please, do it, and come back here later!!  '''

                #The message box which appears is entitled 'No MT3D-USGS executable in Program Locations '
                QtGui.QMessageBox.warning(None, 'No MT3D-USGS executable in Program Locations ', errorMsg )
                #If this error occurs, nothing is done
                return

            # If yes, proceed...
            #If the LAK package is in the list package_list...
            lktActive = True

            #modelName is the name of the flow model read from the combo box cmbModelName
            modelName = self.cmbModelName.currentText()
            #The tuple (pathfile, nsp ) is made of two objects:
            # -pathfile is the path of the working directiry for a specific flow model
            # -nsp is the number of stress period for a specific flow model
            #(getModelInfoByName is a function of the freewat_utils.py module; such function
            #recalls the getModelsInfoLists function to retrieve the names of all the flow
            #models available in the legend and the corresponding paths, then it retrieves
            #the path of the specific flow model selected in the cmbFlowName combo box and
            #iterates over the rows of the timetable_modelName table to get the number of
            #stress periods; it returns a tuple)
            (pathfile, nsp ) = getModelInfoByName(modelName)

            # Get info from LKT Tables
            #lkt_start_conc is the lkt table whose name is read from the combo box cmbLKT_start_conc
            lkt_start_conc = getVectorLayerByName(self.cmbLKT_start_conc.currentText())

            #nlakes is the number of lakes implemented in the flow model
            nlakes=0
            #coldlak_params is a list containing the starting
            #concentration of species n.1 for each lake
            coldlak_params=[]

            #f iterates over the rows of the table 'lake_starting_conc_'
            for f in lkt_start_conc.getFeatures():
                #nlakes is increased by 1
                nlakes+=1
                ##The list coldlak_params is updated with values read from the field 'coldlak_spec_1'
                coldlak_params.append(f['coldlak_spec_1'])

                #If there are more than one species...
                if nspec > 1:
                    coldlak=[]
                    lkt_multiComponent={}

                    #i iterates over the number of species, starting from species #2
                    for i in range(2,nspec+1):
                        #f iterates over the rows of the table 'lake_starting_conc_'
                        for f in lkt_start_conc.getFeatures():
                            #The list coldlak is updated with values read from the field 'coldlak_spec_i'
                            coldlak.append(f['coldlak_spec_' + str(i)])

                        #The dictionary lkt_multiComponent is updated with pairs (key:value).
                        #The keys are strings coldlaki (where i is the number of the i-th species).
                        #The values are the elements of the list coldlak (i.e., the coldlak values
                        #for the i-th species).
                        #ATTENTION: the elements of the dictionary lkt_multiComponent might not be sorted
                        #by keys! What is important is the name assigned to each key (such name is the same
                        #as that indicated in the module mtlkt.py of the flopy)
                        lkt_multiComponent['coldlak' + str(i)]=coldlak
                        #The list coldlak is emptied
                        coldlak=[]

            #lkt_bc is the lkt table whose name is read from the combo box cmbLKT_bc
            lkt_bc = getVectorLayerByName(self.cmbLKT_bc.currentText())

            bc_params={}
            bc=[]
            bc_list=[]

            #n iterates over the number of stress periods
            for n in range(nsp):
                #f iterates over the rows of the table 'lake_table_bc_'
                for f in lkt_bc.getFeatures():

                    #The list bc_list is updated with values read from the field 'lake_ID'
                    #ATTENTION: f['lake_ID'] is decreased by 1 because in the file .lkt
                    #it is not well written, even if the dictionary bc_params is correct!!!
                    bc.append(f['lake_ID']-1)
                    #The list bc_list is updated with values read from the field 'bc_type'
                    #If the content of the field 'bc_type' is 'precipitation'...
                    if f['bc_type']=='precipitation':
                        #...then the list bc_list is updated with 1
                        bc.append(1)
                    #If the content of the field 'bc_type' is 'runoff'...
                    elif f['bc_type']=='runoff':
                        #...then the list bc_list is updated with 2
                        bc.append(2)
                    #If the content of the field 'bc_type' is 'pumping'...
                    elif f['bc_type']=='pumping':
                        #...then the list bc_list is updated with 3
                        bc.append(3)
                    #If the content of the field 'bc_type' is 'evaporation'...
                    elif f['bc_type']=='evaporation':
                        #...then the list bc_list is updated with 4
                        bc.append(4)
                    #If a different string is read...
                    else:
                        #...then a warning message appears and nothing is done
                        pop_message(self.tr('Invalid bc_type in the lake_table_bc table!\nPlease, check the spelling of bc_type'), self.tr('warning'))
                        return
                    #The list bc_list is updated with values read from the field
                    #'conc_sp_n_spec_1' (where n is the n-th stress period)
                    #(the list bc_list is something like [lake,bc_type,conc]
                    bc.append(f['conc_sp_' + str(n+1) + '_spec_1'])

                    #If there are more than one species...
                    if nspec > 1:

                        #i iterates over the number of species, starting from species #2
                        for i in range(2,nspec+1):
                            #The list bc_list is updated with values read from the field
                            #'conc_sp_n_spec_m' (where n is the n-th stress period and
                            #m is the m.th species)
                            #(the list bc_list is something like [lake,bc_type,conc_species1,conc_species2,...]
                            bc.append(f['conc_sp_' + str(n+1) + '_spec_' + str(i)])
                    bc_list.append(bc)
                    bc=[]

                #The dictionary bc_params is updated (the key is the number of the current
                #stress period; as required by the flopy, the argument is the list bc_list)
                bc_params[n]=bc_list
                #The list bc_list is emptied
                bc_list=[]

            #####Check if concentration values are all null at certain stress periods#####
            #Extract values of the keys of the dictionary bc_params
            #(the elements of the list list_values are lists of lists)
            list_values=bc_params.values()
            summ=[]
            #i iterates over all the possible boundary conditions
            #(i is the second index of all the elements of the list list_values)
            for i in range(nlakes*4):
                count=0.0
                #n iterates over the number of stress periods
                #(n is the first index of all the elements of the list list_values)
                for n in range(nsp):
                    #The float count is updated
                    count+=sum(list_values[n][i][2:])
                #The list summ is updated (at the end, the list summ will contain
                #nlakes*4 float elements. Such floats will be the sum of
                #concentration values at each lake, for each bc type and for each sp)
                summ.append(count)

            #The list new_list_values will be filled with lists only referring to lakes,
            #bc type and sp where at least one non null concentration value exists
            new_list_values=[]
            #The list new_list_values is initialized with as many elements as the list
            #list_values has. Such elements are (at the beginning) lists of null lists
            for i in range(len(list_values)):
                new_list_values.append([])

            bc_Active=False
            #i iterates over the indexes of the list summ
            #(i is the second index of all the elements of the list list_values)
            #count iterates over the elements of the list summ
            for i, count in enumerate(summ):
                #n iterates over the number of stress periods
                #(n is the first index of all the elements of the list list_values)
                #(n is also the first index of all the elements of the list new_list_values)
                for n in range(nsp):
                    #If count is not null (i.e., if at least one non null
                    #concentration value exists for a certain lake, for a
                    #certain bc type and for a certain sp)...
                    if count != 0.0:
                        #...then at least one bc exists
                        bc_Active=True
                        #...and the list new_list_values is updated
                        new_list_values[n].append(list_values[n][i])

            #n iterates over the number of stress periods
            #(n is the index of all the keys of the dictionary bc_params)
            for n in range(nsp):
                #The dictionary bc_params is updated with the elements of the list new_list_values
                bc_params[n]=new_list_values[n]

            # build the LKT package:
            #(mt is an Mt3d object created recalling above the class Mt3dmsUSGS
            #present in the module mt.py within flopy/flopy/mt3d/.
            #ATTENTION: in the  module mt.py within flopy/flopy/mt3d/ the class
            #is actually called Mt3dms, while in the module flopyaddon/mtusgs.py
            #it is called Mt3dmsUSGS).
            #If there are more than one species and there are some boundary conditions...
            if (nspec > 1 and bc_Active == True):
                #...coldlak for the first species is read from the field 'coldlak_spec_1'
                #(see the list coldlak_params above);
                #for species #2 through #n, coldlak is read from the dictionary lkt_multiComponent
                lkt = Mt3dLkt(mt, nlkinit=nlakes, mxlkbc=nlakes*4, icbclk=None, ietlak=0,
                coldlak=coldlak_params, lk_stress_period_data=bc_params,
                unitnumber=18, **lkt_multiComponent)
            #If there are more than one species but no boundary conditions exist...
            elif (nspec > 1 and bc_Active == False):
                #...coldlak for the first species is read from the field 'coldlak_spec_1'
                #(see the list coldlak_params above);
                #for species #2 through #n, coldlak is read from the dictionary lkt_multiComponent
                lkt = Mt3dLkt(mt, nlkinit=nlakes, mxlkbc=nlakes*4, icbclk=None, ietlak=0,
                coldlak=coldlak_params, lk_stress_period_data=None,
                unitnumber=18, **lkt_multiComponent)
            #If there is just one species and there are some boundary conditions...
            elif (nspec == 1 and bc_Active == True):
                #...the dictionary lkt_multiComponent is no longer needed
                lkt = Mt3dLkt(mt, nlkinit=nlakes, mxlkbc=nlakes*4, icbclk=None, ietlak=0,
                coldlak=coldlak_params, lk_stress_period_data=bc_params,
                unitnumber=18)
            #If there is just one species and no boundary conditions exist...
            elif (nspec == 1 and bc_Active == False):
                #...the dictionary lkt_multiComponent is no longer needed
                lkt = Mt3dLkt(mt, nlkinit=nlakes, mxlkbc=nlakes*4, icbclk=None, ietlak=0,
                coldlak=coldlak_params, lk_stress_period_data=None,
                unitnumber=18)

        # 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

            #If the SFT or the LKT packages are going to be run...
            if sftActive or lktActive :
                #...if the MT3D-USGS executable is not present...
                if mtusgsdir == '' :
                    # pop message if the value of the dictionary is empty
                    #...then a warning message appears
                    pop_message(self.tr('MT3D-USGS executable is needed!\nPlease enter the path of such executable in the prg_locations table!'), self.tr('warning'))
                    #...and nothing is done
                    return
                #...otherwise, if the MT3D-USGS executable is present...
                else:
                    #...such executable is retrieved from the prg_locations table
                    mt3dexe = mtusgsdir

            #If the UZT, the SFT ald the LKT packages are not going to be run...
            if not uztActive and not sftActive and not lktActive :
                #...if the MT3DMS executable is not present...
                if mt3ddir == '' :
                    # pop message if the value of the dictionary is empty
                    #...then a warning message appears
                    pop_message(self.tr('MT3DMS executable is needed!\nPlease enter the path of such executable in the prg_locations table!'), self.tr('warning'))
                    #...and nothing is done
                    return
                #...otherwise, if the MT3DMS executable is present...
                else:
                    #...such executable is retrieved from the prg_locations table
                    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)
