# gSSURGO_ValuTable.py
#
# Steve Peaslee, National Soil Survey Center
# 2014-09-27

# Adapted from HorizonData_Mapping
# Calculates root zone depth, root zone available water supply and a set of
# available water supply values at various defined depth ranges.
#
# The root zone available water supply is calculated down to the root zone for
# each component and then a weighted average is calculated for the map unit.
#
# The rest of the AWS columns are calculated for all components at a variety of depths.
# Output is to a geodatabase table (Valu2)
#
# To do list for FY2016....
# In order for this script to meet the 2015 standard for the Valu1 table it needs to:
#   skip components where compkind is null
#   stop horizon calculations for all data where there is a discrepancy in horizon depths
#   null map unit AWS, SOC and PWSL values where sum of all components > 100
#   null Rootzone Depth, Rootzone AWS and NCCPI* where sum of major components > 100
#   perform independent validation of SOC results
#   compare metadata dates to make sure they are properly coded and inserted
#   look at array processing for calculations based upon depth ranges

## ===================================================================================
class MyError(Exception):
    pass

## ===================================================================================
def errorMsg():
    try:
        tb = sys.exc_info()[2]
        tbinfo = traceback.format_tb(tb)[0]
        theMsg = tbinfo + " \n" + str(sys.exc_type)+ ": " + str(sys.exc_value) + " \n"
        PrintMsg(theMsg, 2)

    except:
        PrintMsg("Unhandled error in errorMsg method", 2)
        pass

## ===================================================================================
def PrintMsg(msg, severity=0):
    # Adds tool message to the geoprocessor
    #
    #Split the message on \n first, so that if it's multiple lines, a GPMessage will be added for each line
    try:
        for string in msg.split('\n'):
            #Add a geoprocessing message (in case this is run as a tool)
            if severity == 0:
                arcpy.AddMessage(string)

            elif severity == 1:
                arcpy.AddWarning(string)

            elif severity == 2:
                arcpy.AddMessage("    ")
                arcpy.AddError(string)

    except:
        pass

## ===================================================================================
def Number_Format(num, places=0, bCommas=True):
    try:
    # Format a number according to locality and given places
        #locale.setlocale(locale.LC_ALL, "")
        locale.getlocale()

        if bCommas:
            theNumber = locale.format("%.*f", (places, num), True)

        else:
            theNumber = locale.format("%.*f", (places, num), False)
        return theNumber

    except:
        PrintMsg("Unhandled exception in Number_Format function (" + str(num) + ")", 2)
        return False

## ===================================================================================
def GetLastDate(inputDB):
    # Get the most recent date 'YYYYMMDD' from SACATALOG.SAVEREST and use it to populate metadata
    #
    try:
        tbl = os.path.join(inputDB, "SACATALOG")
        today = ""

        with arcpy.da.SearchCursor(tbl, ['SAVEREST'], "", "", [None, "ORDER_BY SAVEREST DESC"]) as cur:
            for rec in cur:
                #lastDate = rec[0].split(" ")[0].replace("-", "")
                lastDate = rec[0].strftime('%Y%m%d')
                break

        return lastDate


    except MyError, e:
        # Example: raise MyError("this is an error message")
        PrintMsg(str(e) + " \n", 2)
        return ""

    except:
        errorMsg()
        return ""

## ===================================================================================
def CreateQueryTables(inputDB, outputDB, maxD):
    #
    try:

        # Problems with query table not getting all records where join table record is null
        # FoQueryTable_CRr instance, I appear to be missing a component record when there is no chorizon data,
        # such as for 'Miscellaneous area'
        env.workspace = inputDB
        queryMU = "MU"
        queryCO = "CO"
        queryHZ = "HZ"
        queryCR = "CR"
        queryCT = "CT"
        queryTemp = "Tmp"  # empty table uses as template

        # Create an empty table using the original queryHZ as a template

        # Mapunit query
        PrintMsg(" \nReading MAPUNIT table...", 0)
        whereClause = ""
        fldMu = [["mukey", "mukey"], ["musym", "musym"], ["muname", "muname"]]
        #arcpy.MakeQueryTable_management(["mapunit"], queryMU, "USE_KEY_FIELDS", "#", fldMu, whereClause)
        fldMu2 = list()
        dMu = dict()
        #sqlClause = (None, "ORDER BY cokey, resdept_r ASC")
        sqlClause = (None, "ORDER BY mukey")

        for fld in fldMu:
            fldMu2.append(fld[0])

        muList = list()

        #with arcpy.da.SearchCursor(queryMU, fldMu2, "", "", "", sqlClause) as mcur:
        muTbl = os.path.join(inputDB, "mapunit")
        tmukey = 515455

        with arcpy.da.SearchCursor(muTbl, fldMu2, "", "", "", sqlClause) as mcur:
            for mrec in mcur:
                rec = list(mrec)
                mukey = int(rec[0])
                rec.pop(0)
                dMu[mukey] = rec
                muList.append(mukey)

                #if mukey == tmukey:
                #    PrintMsg(" \nFound mukey " + str(tmukey) + " in mapunit table", 1)

        muList.sort()

        # Component query
        PrintMsg(" \nReading COMPONENT table...", 0)
        fldCo = [["mukey", "mukey"], ["cokey", "cokey"], ["comppct_r", "comppct_r"], ["majcompflag", "majcompflag"], \
        ["compname", "compname"], ["compkind", "compkind"], ["taxorder", "taxorder"], ["taxsubgrp", "taxsubgrp"], \
        ["localphase", "localphase"], ["otherph", "otherph"], ["hydricrating", "hydricrating"], ["drainagecl", "drainagecl"]]
        fldCo2 = list()
        dCo = dict()
        whereClause = "comppct_r is not NULL"
        sqlClause = (None, "ORDER BY cokey, comppct_r DESC")
        coTbl = os.path.join(inputDB, "component")

        for fld in fldCo:
            fldCo2.append(fld[0])

        with arcpy.da.SearchCursor(coTbl, fldCo2, whereClause, "", "", sqlClause) as ccur:
            for crec in ccur:
                rec = list(crec)
                mukey = int(rec.pop(0))  # get rid of mukey from component record

                try:
                    # Add next component record to list
                    dCo[mukey].append(rec)

                except:
                    # initialize list of records
                    dCo[mukey] = [rec]

        # HORIZON TABLE
        PrintMsg(" \nReading HORIZON table...", 0)
        fldHz = [["cokey", "cokey"], ["chkey", "chkey"], ["hzname", "hzname"], ["desgnmaster", "desgnmaster"], \
        ["hzdept_r", "hzdept_r"], ["hzdepb_r", "hzdepb_r"], ["sandtotal_r", "sandtotal_r"], \
        ["silttotal_r", "silttotal_r"], ["claytotal_r", "claytotal_r"], ["om_r", "om_r"], \
        ["dbthirdbar_r", "dbthirdbar_r"], ["ec_r", "ec_r"], ["ph1to1h2o_r", "ph1to1h2o_r"], \
        ["awc_r", "awc_r"]]
        fldHz2 = list()
        dHz = dict()
        whereClause = "hzdept_r is not NULL and hzdepb_r is not NULL"
        sqlClause = (None, "ORDER BY chkey, hzdept_r ASC")

        for fld in fldHz:
            fldHz2.append(fld[0])

        hzTbl = os.path.join(inputDB, "chorizon")

        with arcpy.da.SearchCursor(hzTbl, fldHz2, whereClause, "", "", sqlClause) as hcur:
            for hrec in hcur:
                rec = list(hrec)
                cokey = int(rec.pop(0))

                try:
                    # Add next horizon record to list
                    dHz[cokey].append(rec)

                except:
                    # initialize list of horizon records
                    dHz[cokey] = [rec]

            # HORIZON TEXTURE QUERY
            #
            PrintMsg(" \nMaking query table (CT) for texture information", 0)
            inputTbls = list()
            tbls = ["chtexturegrp", "chtexture"]
            for tbl in tbls:
                inputTbls.append(os.path.join(inputDB, tbl))

            txList1 = [["chtexturegrp.chkey", "chkey"], ["chtexturegrp.texture", "texture"], ["chtexture.lieutex", "lieutex"]]
            whereClause = "chtexturegrp.chtgkey = chtexture.chtgkey and chtexturegrp.rvindicator = 'Yes'"
            arcpy.MakeQueryTable_management(inputTbls, queryCT, "USE_KEY_FIELDS", "#", txList1, whereClause)

            # Read texture query into dictionary
            txList2 = ["chtexturegrp.chkey", "chtexturegrp.texture", "chtexture.lieutex"]
            dTexture = dict()
            ctCnt = int(arcpy.GetCount_management(queryCT).getOutput(0))

            arcpy.SetProgressor ("step", "Getting horizon texture information for QueryTable_HZ...", 0, ctCnt, 1)

            # Have been running out of memory here if other applications are running. 5.66GB

            with arcpy.da.SearchCursor(queryCT, txList2) as cur:
                for rec in cur:
                    dTexture[int(rec[0])] = [rec[1], rec[2]]
                    arcpy.SetProgressorPosition()

            arcpy.Delete_management(queryCT)

        # COMPONENT RESTRICTIONS which will be saved to a gdb table
        #
        crTbl = os.path.join(inputDB, "corestrictions")
        PrintMsg(" \nWriting component restrictions to " + os.path.join(outputDB, "QueryTable_CR"), 0)
        fldCr = [["cokey", "cokey"], \
        ["reskind", "reskind"],\
        ["reshard", "reshard"],\
        ["resdept_r", "resdept_r"]]
        whereClause = "OBJECTID = 1"
        arcpy.MakeQueryTable_management(crTbl, queryCR, "USE_KEY_FIELDS", "#", fldCr, whereClause)
        arcpy.CreateTable_management(outputDB, "QueryTable_CR", queryCR)
        arcpy.Delete_management(queryCR)

        fldCr2 = list()
        dCr = dict()
        sqlClause = (None, "ORDER BY cokey, resdept_r ASC")

        for fld in fldCr:
            fldCr2.append(fld[0])

        whereClause = "resdept_r is not NULL and resdept_r < " + str(maxD)
        outputCR = os.path.join(outputDB, "QueryTable_CR")
        crList = list() # list of components with a restriction
        crCnt = int(arcpy.GetCount_management(crTbl).getOutput(0))
        arcpy.SetProgressor ("step", "Saving component restriction information...", 0, crCnt, 1)

        with arcpy.da.SearchCursor(crTbl, fldCr2, whereClause, "", "", sqlClause) as crcur:

            for crrec in crcur:
                rec = list(crrec)
                cokey = int(rec[0])
                arcpy.SetProgressorPosition()

                if not cokey in crList:
                    # Only save the highest level restriction above 150cm
                    crList.append(cokey)
                    dCr[cokey] = rec

        PrintMsg(" \nCreating QueryTable_HZ in " + outputDB, 0)
        fldCo.pop(0)
        fldCo2.pop(0)
        fldHz.pop(0)
        fldHz2.pop(0)

        # Create list of fields for query table
        fldAll = list()
        # Create list of fields for output cursor
        fldAll2 = list()

        for fld in fldMu:
            fldAll.append(["mapunit." + fld[0], fld[1]])
            fldAll2.append(fld[1])

        for fld in fldCo:
            fldAll.append(["component." + fld[0], fld[1]])
            fldAll2.append(fld[1])

        for fld in fldHz:
            fldAll.append(["chorizon." + fld[0], fld[1]])
            fldAll2.append(fld[1])

        # Texture fields:
        fldAll2.append("texture")
        fldAll2.append("lieutex")

        # Select component-horizon data for ALL components that have horizon data. Lack of horizon data
        # will cause some components to be missing from the PWSL.
        #
        # Later on in the actual calculations for RZAWS, only the major-earthy components will be used. But
        # all components are in this table!

        whereClause = "mapunit.mukey = component.mukey and \
        component.cokey = chorizon.cokey and mapunit.objectid = 1"

        outputTable = os.path.join(outputDB, "QueryTable_HZ")
        PrintMsg(" \nCreating table " + outputTable, 0)
        arcpy.MakeQueryTable_management(['mapunit', 'component', 'chorizon'], queryTemp, "USE_KEY_FIELDS", "#", fldAll, whereClause)
        arcpy.CreateTable_management(outputDB, "QueryTable_HZ", queryTemp)
        arcpy.AddField_management(outputTable, "texture", "TEXT", "", "", "30", "texture")
        arcpy.AddField_management(outputTable, "lieutex", "TEXT", "", "", "254", "lieutex")
        arcpy.Delete_management(queryTemp)

        # Process dictionaries and use them to write out the new QueryTable_HZ table
        #
        # Open output table
        outFld2 = arcpy.Describe(outputTable).fields
        outFlds = list()
        for fld in outFld2:
            outFlds.append(fld.name)

        outFlds.pop(0)

        # Create empty lists to replace missing data
        missingCo = ["", None, None, None, None, None, None, None, None, None, None]
        missingHz = ["", None, None, None, None, None, None, None, None, None, None, None, None]
        missingTx = [None, None]

        # Save information on mapunits or components with bad or missing data
        #badMu = list()   # list of mapunits with no components
        muNoCo = list()
        dNoCo = dict()

        coNoHz = list()  # list of components with no horizons
        dNoHz = dict() # component data for those components in coNoHz

        arcpy.SetProgressor ("step", "Writing data to " + outputTable + "...", 0, len(muList), 1)

        with arcpy.da.InsertCursor(outputTable, fldAll2) as ocur:

            for mukey in muList:
                mrec = dMu[mukey]
                arcpy.SetProgressorPosition()

                try:
                    coVals = dCo[mukey]  # got component records for this mapunit

                    # Sort lists by comppct_r
                    coList = sorted(coVals, key = lambda x: int(x[1]))

                    for corec in coList:
                        cokey = int(corec[0])
                        #if mukey == tmukey:
                        #    PrintMsg("\t\tCO1: " + str(cokey) + ": " + str(corec), 1)

                        try:
                            hzVals = dHz[cokey]  # horizon records for this component
                            # Sort record by hzdept_r
                            hzList = sorted(hzVals, key = lambda x: int(x[3]))

                            for hzrec in hzList:
                                chkey = int(hzrec[0])

                                #if mukey == tmukey:
                                #    PrintMsg("\t\t\tHZ1: " + str(chkey) + ": " + str(hzrec), 1)

                                try:
                                    # Get horizon texture
                                    txrec = dTexture[chkey]

                                except:
                                    txrec = missingTx

                                # Combine all records and write to table
                                newrec = [mukey]
                                newrec.extend(mrec)
                                newrec.extend(corec)
                                newrec.extend(hzrec)
                                newrec.extend(txrec)
                                #if mukey == tmukey:
                                #    PrintMsg("\t\t\tHZ2: " + str(cokey) + ": " + str(corec), 1)

                                ocur.insertRow(newrec)

                        except KeyError:
                            # No horizon records for this component
                            comppct = corec[1]
                            compname = corec[3]
                            compkind = corec[4]
                            mjrcomp = corec[2]
                            #PrintMsg("Major compflag = " + str(corec), 1)

                            hzrec = missingHz
                            txrec = missingTx
                            newrec = [mukey]
                            newrec.extend(mrec)
                            newrec.extend(corec)
                            newrec.extend(hzrec)
                            newrec.extend(txrec)
                            #if mukey == tmukey:
                                #PrintMsg("\t\tNo horizon data for mapunit:component: " + str(mukey) + ":" + str(cokey), 1)
                                #PrintMsg("\t\t\tHZ3: " + str(cokey) + ": " + str(newrec), 1)
                            ocur.insertRow(newrec)

                            if not (compname in ["NOTCOM", "NOTPUB"] or compkind == 'Miscellaneous area'):
                                badComp = [mukey, str(cokey), compname, compkind, mjrcomp, str(comppct)]
                                coNoHz.append(str(cokey))   # add cokey to list of components with no horizon data
                                dNoHz[cokey] = badComp      # add component information to dictionary
                                #PrintMsg(" \nMissing horizon data: " + str(corec), 1)

                        except:
                            PrintMsg(" \nhzVals error for " + str(mukey) + ":" + str(cokey) + ": " + str(txrec), 2)
                            PrintMsg(" \n" + str(fldAll2), 1)
                            errorMsg()

                except:
                    # No component records for this map unit
                    corec = missingCo
                    hzrec = missingHz
                    txrec = missingTx
                    newrec = [mukey]
                    newrec.extend(mrec)
                    newrec.extend(corec)
                    newrec.extend(hzrec)
                    newrec.extend(txrec)
                    #if mukey == tmukey:
                        #PrintMsg("\t\tNo component data for mapunit:component: " + str(mukey) + ":" + str(cokey), 1)
                    ocur.insertRow(newrec)

                    if not  mrec[0] in ['NOTCOM', 'NOTPUB']:
                        # skip map units that should never have component data
                        #
                        muNoCo.append(str(mukey))
                        dNoCo[str(mukey)] = [mrec[0], mrec[1]] # Save map unit name for the report
                        #PrintMsg(" \n\n** No component data for " + str(mrec[1]), 2)


        # Run through QueryTbl_HZ table, checking for inconsistencies in horizon depths
        # Create a dictionary containing a list of top and bottom of each horizon in each component
        # dictionary key = cokey
        # list contains tuples of hzdept_r, hzdepb_r, hzname, mukey, compname, localphase

        dHZ = dict()

        # Exclude horizon data with null hzdep with whereclause
        wc = "hzdept_r is not null and hzdepb_r is not null"
        arcpy.ResetProgressor()
        arcpy.SetProgressorLabel("Looking for inconsistencies in horizon depths...")

        with arcpy.da.SearchCursor(outputTable, ['mukey', 'cokey','hzdept_r','hzdepb_r', 'hzname', 'compname', 'localphase', 'majcompflag'], wc) as cur:
            for rec in cur:
                mukey, cokey, top, bot, hzname, compname, localphase, majcomp = rec
                cokey = int(cokey)

                if cokey in dHZ:
                    vals = dHZ[cokey]
                    vals.append([top, bot, hzname, mukey, compname, localphase, majcomp])
                    dHZ[cokey] = vals

                else:
                    dHZ[cokey] = [[top, bot, hzname, mukey, compname, localphase, majcomp]]

        # Number of items in each component value = len(rec)
        # top of first horizon = rec[0][0]
        # bottom of last horizon = rec[len(rec) - 1][1]
        # component thickness #1 = rec[len(rec) - 1][1] - rec[0][0]
        # Read each entry in the dictionary and check for gaps and overlaps in the horizons
        #
        badCoHz = list()
        badHorizons = list()

        for cokey, vals in dHZ.items():
            # Process each component
            #
            hzSum = 0                     # itialize sum of horizon thicknesses
            lb = vals[0][0]               # initialize last bottom to the top of the first horizon
            localphase = vals[0][5]

            if localphase is None:
                localphase = ""

            else:
                localphase = " " + localphase

            for v in vals:
                # Process each horizon in the component record
                #
                # sum of bottom - top for each horizon
                hzSum += (v[1] - v[0])

        		# Check for consistency between tops and bottoms for each consecutive horizon
                if v[0] != lb:
                    diff = v[0] - lb
                    badCoHz.append(str(cokey))
                    badHorizons.append(v[3] + ", " + str(cokey) + ", " + v[4] + localphase + ", " + majcomp + ", " + str(v[2]) + ", " + str(v[0]) + ", " + str(diff) )

                lb = v[1] # update last bottom depth

        PrintMsg(" \nWriting component restrictions to " + outputCR, 0)
        arcpy.SetProgressor ("step", "Writing component restriction data to " + outputCR + "...", 0, len(dCr), 1)

        with arcpy.da.InsertCursor(outputCR, fldCr2) as ocur:
            for cokey, crrec in dCr.items():
                ocur.insertRow(crrec)

        # Save data issues to permanent files for later review
        #
        # Don't report these errors from this script tool. This code has been moved to another tool
        if 1 == 2:  # Skip this code. May use some of it to null out some valu table results.
        #
        #if len(muNoCo) > 0 or len(coNoHz) or len(badCoHz) > 0:
            PrintMsg(" \nCreating log file: " + logFile, 1)
            now = datetime.now()
            fh = open(logFile, "w")
            fh.write("\n" + inputDB + "\n")
            fh.write("\nProcessed on " + now.strftime('%A %x  %X') + "\n\n")
            fh.write("This log file containins record of any data inconsistencies found while processing the Valu2 table \n\n")
            fh.close()

            # Save data issues (mapunits with no components) to log file for later review
            # Some of these will be NOTCOMs
            if len(muNoCo) > 0:
                fh = open(logFile, "a")
                #fh.write(inputDB + "\n")
                fh.write("\nQuery for map units missing component data\n")
                fh.write("====================================================================================\n")
                fh.write("MUKEY IN ('" + "', '".join(muNoCo) + "') \n\n")
                fh.write("\n\nTable of map units missing component data\n")
                fh.write("\nMUKEY, MUSYM, MUNAME\n")
                for mukey in muNoCo:
                    fh.write(str(mukey) + ", " + dNoCo[mukey][0] + ", " + dNoCo[mukey][1] + "\n")

                fh.close()
                PrintMsg(" \nMap units missing component data (" + Number_Format(len(muNoCo), 0, True) + ") saved to:\t" + logFile, 0)

            # Save data issues (components with no horizons) to log file for later review
            # Note; these COKEYs will work with gSSURGO but not Soil Data Access
            # Some of these may be NOTCOMs
            if len(coNoHz) > 0:
                PrintMsg(" \nQuery for components missing horizon data (" + Number_Format(len(coNoHz), 0, True) + ") saved to:\t" + logFile, 0)
                fh = open(logFile, "a")
                fh.write("\n\nQuery for components with no horizon data\n")
                fh.write("====================================================================================\n")
                fh.write("COKEY IN ('" + "', '".join(coNoHz) + "') \n\n")

                fh.write("Table of map-unit\components missing horizon data\n")
                fh.write("\nMUKEY, COKEY, COMPNAME, COMPKIND, MAJCOMPFLAG, COMPPCT, COMPKIND\n")

                for mukey, compInfo in dNoHz.items():
                    # mukey, str(cokey), compname, compkind, mjrcomp, str(comppct)
                    mukey, cokey, compname, compkind, majcomp, comppct = compInfo
                    #compInfo = mukey, str(cokey), compname, compkind, str(comppct)
                    fh.write(str(mukey) + ", " + str(cokey) + ", " + compname + ", " + str(compkind) + ", "  + majcomp + ", " + str(comppct) + ", " + str(compkind) + "\n")

                fh.close()

            # These components have inconsistent horizon depths that overlap or gap.
            # Note; these COKEYs will work with gSSURGO but not Soil Data Access
            if len(badCoHz) > 0:
                PrintMsg(" \nComponents with horizon gaps or overlaps (" + Number_Format(len(badCoHz), 0, True) + ") saved to:\t" + logFile, 0)
                fh = open(logFile, "a")
                fh.write("\n\nQuery for components with horizon gaps or overlaps\n")
                fh.write("====================================================================================\n")
                fh.write("COKEY IN ('" + "', '".join(badCoHz) + "') \n\n")
                fh.write("\nTable of components with horizon gaps or overlaps\n")
                fh.write("MUKEY, COKEY, COMPNAME, MJRCOMP, HZNAME, HZDEPT, DIFF\n")

                for h in badHorizons:
                    fh.write(h + "\n")

                fh.close()

        arcpy.ResetProgressor()

        env.workspace = outputDB
        return True

    except MyError, e:
        # Example: raise MyError("this is an error message")
        PrintMsg(str(e) + " \n", 2)
        return False

    except:
        errorMsg()
        return False

## ===================================================================================
def CreateOutputTableMu(theMuTable, depthList, dPct):
    # Create the mapunit level table
    #
    try:
        # Create the output tables and add required fields

        try:
            # Try to handle existing output table if user has added it to ArcMap from a previous run
            if arcpy.Exists(theMuTable):
                arcpy.Delete_management(theMuTable)

        except:
            raise MyError, "Previous output table (" + theMuTable + ") is in use and cannot be removed"
            return False

        PrintMsg(" \nCreating new output table (" + os.path.basename(theMuTable) + ") for mapunit level data", 0)
        outputDB = os.path.dirname(theMuTable)
        tmpTable = os.path.join("IN_MEMORY", os.path.basename(theMuTable))

        #arcpy.CreateTable_management(outputDB, os.path.basename(theMuTable))
        arcpy.CreateTable_management("IN_MEMORY", os.path.basename(theMuTable))

        # Add fields for AWS
        for rng in depthList:
            # Create the AWS fields in a loop
            #
            td = rng[0]
            bd = rng[1]
            awsField = "AWS" + str(td) + "_" + str(bd)
            arcpy.AddField_management(tmpTable, awsField, "SHORT", "", "", "", awsField)

        for rng in depthList:
            # Create the AWS fields in a loop
            #
            td = rng[0]
            bd = rng[1]
            awsField = "TK" + str(td) + "_" + str(bd) + "A"
            arcpy.AddField_management(tmpTable, awsField, "FLOAT", "", "", "", awsField)

        arcpy.AddField_management(tmpTable, "MUSUMCPCTA", "SHORT", "", "", "")


        # Add Fields for SOC
        for rng in depthList:
            # Create the SOC fields in a loop
            #
            td = rng[0]
            bd = rng[1]
            socField = "SOC" + str(td) + "_" + str(bd)
            arcpy.AddField_management(tmpTable, socField, "LONG", "", "", "", socField)

        for rng in depthList:
            # Create the SOC fields in a loop
            #
            td = rng[0]
            bd = rng[1]
            socField = "TK" + str(td) + "_" + str(bd) + "S"
            arcpy.AddField_management(tmpTable, socField, "FLOAT", "", "", "", socField)

        arcpy.AddField_management(tmpTable, "MUSUMCPCTS", "SHORT", "", "", "")

        # Add fields for NCCPI
        arcpy.AddField_management(tmpTable, "NCCPI2CS", "FLOAT", "", "", "")
        arcpy.AddField_management(tmpTable, "NCCPI2SG", "FLOAT", "", "", "")
        arcpy.AddField_management(tmpTable, "NCCPI2CO", "FLOAT", "", "", "")
        arcpy.AddField_management(tmpTable, "NCCPI2ALL", "FLOAT", "", "", "")

        # Add fields for root zone depth and root zone available water supply
        arcpy.AddField_management(tmpTable, "PCTEARTHMC", "SHORT", "", "", "")
        arcpy.AddField_management(tmpTable, "ROOTZNEMC", "SHORT", "", "", "")
        arcpy.AddField_management(tmpTable, "ROOTZNAWS", "SHORT", "", "", "")
        # Add field for droughty soils
        arcpy.AddField_management(tmpTable, "DROUGHTY", "SHORT", "", "", "")

        # Add field for potential wetland soils
        arcpy.AddField_management(tmpTable, "PWSL1POMU", "SHORT", "", "", "")

        # Add field for mapunit-sum of ALL component-comppct_r values
        arcpy.AddField_management(tmpTable, "MUSUMCPCT", "SHORT", "", "", "")

        # Add Mukey field (primary key)
        arcpy.AddField_management(tmpTable, "MUKEY", "TEXT", "", "", "30", "MUKEY")

        # Convert IN_MEMORY table to a permanent table
        arcpy.CreateTable_management(outputDB, os.path.basename(theMuTable), tmpTable)
        # Add attribute indexes for key fields
        arcpy.AddIndex_management(theMuTable, "MUKEY", "Indx_ResMukey", "NON_UNIQUE", "NON_ASCENDING")

        arcpy.Delete_management(os.path.join("IN_MEMORY", os.path.basename(theMuTable)))

        # Reading from the original mapunit table, populate the output Valu2 table with mukey and musumpct
        #PrintMsg(" \n\tPopulating " + theMuTable + " with mukey values", 1)
        with arcpy.da.SearchCursor(os.path.join(outputDB, "mapunit"), ["mukey"]) as incur:
            outcur = arcpy.da.InsertCursor(theMuTable, ["mukey", "musumcpct"])
            for inrec in incur:
                mukey = inrec[0]
                try:
                    sumPct = dPct[mukey][0]

                except:
                    sumPct = 0
                inrec = [mukey, sumPct]
                outcur.insertRow(inrec)

        return True

    except MyError, e:
        # Example: raise MyError("this is an error message")
        PrintMsg(str(e) + " \n", 2)
        return False

    except:
        errorMsg()
        return False

## ===================================================================================
def CreateOutputTableCo(theCompTable, depthList):
    # Create the component level table
    # The new input field is created using adaptive code from another script.
    #
    try:
        # Create two output tables and add required fields
        try:
            # Try to handle existing output table if user has added it to ArcMap from a previous run
            if arcpy.Exists(theCompTable):
                arcpy.Delete_management(theCompTable)

        except:
            raise MyError, "Previous output table (" + theCompTable + ") is in use and cannot be removed"
            return False

        PrintMsg(" \nCreating new output table (" + os.path.basename(theCompTable) + ") for component level data", 0)

        outputDB = os.path.dirname(theCompTable)
        tmpTable = os.path.join("IN_MEMORY", os.path.basename(theCompTable))

        arcpy.CreateTable_management("IN_MEMORY", os.path.basename(theCompTable))

        # Add fields appropriate for the component level restrictions
        # mukey,cokey, compName, localphase, compPct, comppct, resdept, restriction

        arcpy.AddField_management(tmpTable, "COKEY", "TEXT", "", "", "30", "COKEY")
        arcpy.AddField_management(tmpTable, "COMPNAME", "TEXT", "", "", "60", "COMPNAME")
        arcpy.AddField_management(tmpTable, "LOCALPHASE", "TEXT", "", "", "40", "LOCALPHASE")
        arcpy.AddField_management(tmpTable, "COMPPCT_R", "SHORT", "", "", "", "COMPPCT_R")

        for rng in depthList:
            # Create the AWS fields in a loop
            #
            td = rng[0]
            bd = rng[1]
            awsField = "AWS" + str(td) + "_" + str(bd)
            arcpy.AddField_management(tmpTable, awsField, "FLOAT", "", "", "", awsField)


        for rng in depthList:
            # Create the AWS fields in a loop
            #
            td = rng[0]
            bd = rng[1]
            awsField = "TK" + str(td) + "_" + str(bd) + "A"
            arcpy.AddField_management(tmpTable, awsField, "FLOAT", "", "", "", awsField)

        arcpy.AddField_management(tmpTable, "MUSUMCPCTA", "SHORT", "", "", "")

        for rng in depthList:
            # Create the SOC fields in a loop
            #
            td = rng[0]
            bd = rng[1]
            awsField = "SOC" + str(td) + "_" + str(bd)
            arcpy.AddField_management(tmpTable, awsField, "FLOAT", "", "", "")

        for rng in depthList:
            # Create the rest of the SOC thickness fields in a loop
            #
            td = rng[0]
            bd = rng[1]
            awsField = "TK" + str(td) + "_" + str(bd) + "S"
            arcpy.AddField_management(tmpTable, awsField, "FLOAT", "", "", "")
            arcpy.AddField_management(tmpTable, "MUSUMCPCTS", "SHORT", "", "", "")

        # Root Zone and root zone available water supply
        arcpy.AddField_management(tmpTable, "PCTEARTHMC", "SHORT", "", "", "")
        arcpy.AddField_management(tmpTable, "ROOTZNEMC", "SHORT", "", "", "")
        arcpy.AddField_management(tmpTable, "ROOTZNAWS", "SHORT", "", "", "")
        arcpy.AddField_management(tmpTable, "RESTRICTION", "TEXT", "", "", "254", "RESTRICTION")

        # Droughty soils
        arcpy.AddField_management(tmpTable, "DROUGHTY", "SHORT", "", "", "")

        # Add field for potential wetland soils
        arcpy.AddField_management(tmpTable, "PWSL1POMU", "SHORT", "", "", "")

        # Add primary key field
        arcpy.AddField_management(tmpTable, "MUKEY", "TEXT", "", "", "30", "MUKEY")

        # Convert IN_MEMORY table to a permanent table
        arcpy.CreateTable_management(outputDB, os.path.basename(theCompTable), tmpTable)

        # add attribute indexes for key fields
        arcpy.AddIndex_management(theCompTable, "MUKEY", "Indx_Res2Mukey", "NON_UNIQUE", "NON_ASCENDING")
        arcpy.AddIndex_management(theCompTable, "COKEY", "Indx_ResCokey", "UNIQUE", "NON_ASCENDING")

        # populate table with mukey values
        #PrintMsg(" \n\tPopulating " + theCompTable + " with basic component values", 1)
        with arcpy.da.SearchCursor(os.path.join(outputDB, "component"), ["mukey", "cokey", "compname", "localphase", "comppct_r"]) as incur:
            outcur = arcpy.da.InsertCursor(theCompTable, ["mukey", "cokey", "compname", "localphase", "comppct_r"])
            for inrec in incur:
                outcur.insertRow(inrec)

    except MyError, e:
        # Example: raise MyError("this is an error message")
        PrintMsg(str(e) + " \n", 2)
        return False

    except:
        errorMsg()
        return False


## ===================================================================================
def CheckTexture(mukey, cokey, desgnmaster, om, texture, lieutex, taxorder, taxsubgrp):
    # Is this an organic horizon? Look at desgnmaster and OM first. If those
    # don't help, look at chtexturegrp.texture next.
    #
    # if True: Organic, exclude from root zone calculations unless it is 'buried'
    # if False: Mineral, include in root zone calculations
    #
    # 01-26-2015
    #
    # According to Bob, if TAXORDER = 'Histosol' and DESGNMASTER = 'O' or 'L' then it should NOT be included in the RZAWS calculations
    #
    # If desgnmast = 'O' or 'L' and not (TAXORDER = 'Histosol' OR TAXSUBGRP like 'Histic%') then exclude this horizon from all RZAWS calcualtions.
    #
    # lieutext values: Slightly decomposed plant material, Moderately decomposed plant material,
    # Bedrock, Variable, Peat, Material, Unweathered bedrock, Sand and gravel, Mucky peat, Muck,
    # Highly decomposed plant material, Weathered bedrock, Cemented, Gravel, Water, Cobbles,
    # Stones, Channers, Parachanners, Indurated, Cinders, Duripan, Fragmental material, Paragravel,
    # Artifacts, Boulders, Marl, Flagstones, Coprogenous earth, Ashy, Gypsiferous material,
    # Petrocalcic, Paracobbles, Diatomaceous earth, Fine gypsum material, Undecomposed organic matter

    # According to Bob, any of the 'decomposed plant material', 'Muck, 'Mucky peat, 'Peat', 'Coprogenous earth' LIEUTEX
    # values qualify.
    #
    # This function does not determine whether the horizon might be a buried organic. That is done in CalcRZAWS1.
    #

    lieuList = ['Slightly decomposed plant material', 'Moderately decomposed plant material', \
    'Highly decomposed plant material', 'Undecomposed plant material', 'Muck', 'Mucky peat', \
    'Peat', 'Coprogenous earth']
    txList = ["CE", "COP-MAT", "HPM", "MPM", "MPT", "MUCK", "PDOM", "PEAT", "SPM", "UDOM"]

    try:

        if str(taxorder) == 'Histosols' or str(taxsubgrp).lower().find('histic') >= 0:
            # Always treat histisols and histic components as having all mineral horizons
            #if mukey == tmukey:
            #    PrintMsg("\tHistisol or histic: " + cokey + ", " + str(taxorder) + ", " + str(taxsubgrp), 1)
            return False

        elif desgnmaster in ["O", "L"]:
            # This is an organic horizon according to CHORIZON.DESGNMASTER OR OM_R
            #if mukey == tmukey:
            #    PrintMsg("\tO: " + cokey + ", " + str(taxorder) + ", " + str(taxsubgrp), 1)
            return True

        #elif om > 19:
            # This is an organic horizon according to CHORIZON.DESGNMASTER OR OM_R
        #    if mukey == tmukey:
        #        PrintMsg("\tHigh om_r: " + cokey + ", " + str(taxorder) + ", " + str(taxsubgrp), 1)
        #    return True

        elif str(texture) in txList:
            # This is an organic horizon according to CHTEXTUREGRP.TEXTURE
            #if mukey == tmukey:
            #    PrintMsg("\tTexture: " + cokey + ", " + str(taxorder) + ", " + str(taxsubgrp), 1)
            return True

        elif str(lieutex) in lieuList:
            # This is an organic horizon according to CHTEXTURE.LIEUTEX
            #if mukey == tmukey:
            #    PrintMsg("\tLieutex: " + cokey + ", " + str(taxorder) + ", " + str(taxsubgrp), 1)
            return True

        else:
            # Default to mineral horizon if it doesn't match any of the criteria
            #if mukey == tmukey:
            #    PrintMsg("\tDefault mineral: " + cokey + ", " + str(taxorder) + ", " + str(taxsubgrp), 1)
            return False

    except MyError, e:
        # Example: raise MyError("this is an error message")
        PrintMsg(str(e) + " \n", 2)
        return False

    except:
        errorMsg()
        return False

## ===================================================================================
def CheckBulkDensity(sand, silt, clay, bd, mukey, cokey):
    # Bob's check for a dense layer
    # If sand, silt or clay are missing then we default to Dense layer = False
    # If the sum of sand, silt, clay are less than 100 then we default to Dense layer = False
    # If a single sand, silt or clay value is NULL, calculate it

    try:

        #if mukey == tmukey:
        #    PrintMsg("\tCheck for Dense: " + str(mukey) + ", " + str(cokey) + ", " + \
        #    str(sand) + ", " + str(silt) + ", " + str(clay) + ", " + str(bd), 1)

        txlist = [sand, silt, clay]

        if bd is None:
            # This is not a Dense Layer
            #if mukey == tmukey:
            #    PrintMsg("\tMissing bulk density", 1)
            return False

        if txlist.count(None) == 1:
            # Missing a single total_r value, calculate it
            if txlist[0] is None:
                sand = 100.0 - silt - clay

            elif silt is None:
                silt = 100.0 - sand - clay

            else:
                clay = 100.0 - sand - silt

            txlist = [sand, silt, clay]

        if txlist.count(None) > 0:
            # Null values for more than one, return False
            #if mukey == tmukey:
            #    PrintMsg("\tDense layer with too many null texture values", 1)
            return False

        if round(sum(txlist), 1) <> 100.0:
            # Cannot run calculation, default value is False
            #if mukey == tmukey:
            #    PrintMsg("\tTexture values do not sum to 100", 1)
            return False

        # All values required to run the Dense Layer calculation are available

        a = bd - ((( sand * 1.65 ) / 100.0 ) + (( silt * 1.30 ) / 100.0 ) + (( clay * 1.25 ) / 100.0))

        b = ( 0.002081 * sand ) + ( 0.003912 * silt ) + ( 0.0024351 * clay )

        if a > b:
            # This is a Dense Layer
            #if mukey == tmukey:
            #    PrintMsg("\tDense layer: a = " + str(a) + " and   b = " + str(b) + " and BD = " + str(bd), 1)

            return True

        else:
            # This is not a Dense Layer
            #if mukey == tmukey:
            #    PrintMsg("\tNot a Dense layer: a = " + str(a) + " and   b = " + str(b) + " and BD = " + str(bd), 1)

            return False

    except MyError, e:
        # Example: raise MyError("this is an error message")
        PrintMsg(str(e) + " \n", 2)
        return False

    except:
        errorMsg()
        return False

## ===================================================================================
def CalcRZDepth(inputDB, outputDB, theCompTable, theMuTable, maxD, dPct, dCR):
    #
    # Look at soil horizon properties to adjust the root zone depth.
    # This is in addition to the standard component restrictions
    #
    # Read the component restrictions into a dictionary, then read through the
    # QueryTable_Hz table, calculating the final component rootzone depth
    #
    # Components with COMPKIND = 'Miscellaneous area' or NULL are filtered out.
    # Components with no horizon data are assigned a root zone depth of zero.
    #
    # Horizons with NULL hzdept_r or hzdepb_r are filtered out
    # Horizons with hzdept_r => hzdepb_r are filtered out
    # O horizons or organic horizons from the surface down to the first mineral horizon
    # are filtered out.
    #
    # Horizon data below 150cm or select component restrictions are filtered out.
    # A Dense layer calculation is also included as an additional horizon-specific restriction.

    try:
        dComp = dict()      # component level data for all component restrictions
        dComp2 = dict()     # store all component level data plus default values
        coList = list()

        #bInput = CreateQueryTables1(inputDB, outputDB)
        queryTbl = os.path.join(outputDB, "QueryTable_Hz")

        #if bInput == False:
        #    raise MyError, "Problem querying input database"

        PrintMsg(" \nProcessing input table (" + queryTbl + ")")

        # Create dictionaries and lists
        dMapunit = dict()   # store mapunit weighted restriction depths

        # FIELDS LIST FOR INPUT TABLE
        # areasymbol, mukey, musym, muname, mukname,
        # cokey, compct, compname, compkind, localphase,
        # taxorder, taxsubgrp, ec, pH, dbthirdbar, hzname,
        # hzdesgn, hzdept, hzdepb, hzthk, sand,
        # silt, clay, om, reskind, reshard,
        # resdept, resthk, texture, lieutex

        # All reskind values: Strongly contrasting textural stratification, Lithic bedrock, Densic material,
        # Ortstein, Permafrost, Paralithic bedrock, Cemented horizon, Undefined, Fragipan, Plinthite,
        # Abrupt textural change, Natric, Petrocalcic, Duripan, Densic bedrock, Salic,
        # Human-manufactured materials, Sulfuric, Placic, Petroferric, Petrogypsic
        #
        # Using these restrictions:
        # Lithic bedrock, Paralithic bedrock, Densic bedrock, Fragipan, Duripan, Sulfuric

        # Other restrictions include pH < 3.5 and EC > 16

        resTbl = os.path.join(outputDB, "QueryTable_CR")
        crFlds = ["cokey","reskind", "reshard", "resdept_r"]
        sqlClause = (None, "ORDER BY cokey, resdept_r ASC")

        # ********************************************************
        #
        # Read the QueryTable_HZ and adjust the component restrictions for additional
        # issues such as pH, EC, etc.
        #
        # Save these new restriction values to dComp dictionary
        #
        # Only process major-earthy components...
        whereClause = "component.compkind <> 'Miscellaneous area' and component.compkind is not Null and component.majcompflag = 'Yes'"

        sqlClause = (None, "ORDER BY mukey, comppct_r DESC, cokey, hzdept_r ASC")
        curFlds = ["mukey", "cokey", "compname", "compkind", "localphase", "comppct_r", "taxorder", "taxsubgrp", "hzname", "desgnmaster", "hzdept_r", "hzdepb_r", "sandtotal_r", "silttotal_r", "claytotal_r", "om_r", "dbthirdbar_r", "ph1to1h2o_r", "ec_r", "awc_r", "texture", "lieutex"]
        resList = ['Lithic bedrock','Paralithic bedrock','Densic bedrock', 'Fragipan', 'Duripan', 'Sulfuric']

        lastCokey = "xxxx"
        lastMukey = 'xxxx'

        # Display status of processing input table containing horizon data and component restrictions
        inCnt = int(arcpy.GetCount_management(queryTbl).getOutput(0))

        if inCnt > 0:
            arcpy.SetProgressor ("step", "Processing input table...", 0, inCnt, 1)

        else:
            raise MyError, "Input table contains no data"

        with arcpy.da.SearchCursor(queryTbl, curFlds, whereClause, "", "", sqlClause) as cur:
            # Reading horizon-level data
            for rec in cur:

                # ********************************************************
                #
                # Read QueryTable_HZ record
                mukey, cokey, compName, compKind, localPhase, compPct, taxorder, taxsubgrp, hzname, desgnmaster, hzDept, hzDepb, sand, silt, clay, om, bd, pH, ec, awc, texture, lieutex = rec

                # Initialize component restriction depth to maxD
                dComp2[cokey] = [mukey, compName, localPhase, compPct, maxD, ""]

                if lastCokey != cokey:
                    # Accumulate a list of components for future use
                    lastCokey = cokey
                    coList.append(cokey)

                if hzDept < maxD:
                    # ********************************************************
                    # For horizons above the floor level (maxD), look for other restrictive
                    # layers based on horizon properties such as pH, EC and bulk density.
                    # Start with the top horizons and work down.

                    # initialize list of restrictions
                    resKind = ""
                    restriction = list()

                    bOrganic = CheckTexture(mukey, cokey, desgnmaster, om, texture, lieutex, taxorder, taxsubgrp)

                    if not bOrganic:
                        # calculate alternate dense layer per Dobos
                        bDense = CheckBulkDensity(sand, silt, clay, bd, mukey, cokey)

                        if bDense:
                            # use horizon top depth for the dense layer
                            restriction.append("Dense")
                            resDept = hzDept

                        # Not sure whether these horizon property checks should be skipped for Organic
                        # Bob said to only skip Dense Layer check, but VALU table RZAWS looks like all
                        # horizon properties were skipped.
                        #
                        # If we decide to skip EC and pH horizon checks for histosols/histic, use this query
                        # Example Pongo muck in North Carolina that have low pH but no other restriction
                        #
                        if str(taxorder) != 'Histosols' and str(taxsubgrp).lower().find('histic') == -1:
                            # Only non histosols/histic soils will be checked for pH or EC restrictive horizons
                            if pH <= 3.5 and pH is not None:
                                restriction.append("pH")
                                resDept = hzDept
                                #if mukey == tmukey:
                                #    PrintMsg("\tpH restriction at " + str(resDept) + "cm", 1)

                        if ec >= 16.0 and ec is not None:
                            # Originally I understood that EC > 12 is a restriction, but Bob says he is
                            # now using 16.
                            restriction.append("EC")
                            resDept = hzDept
                            #if mukey == tmukey:
                            #    PrintMsg("\tEC restriction at " + str(resDept) + "cm", 1)

                        #if bd >= 1.8:
                        #    restriction.append("BD")
                        #    resDept = hzDept

                        #if awc is None:
                        #    restriction.append("AWC")
                        #    resDept = hzDept

                    # ********************************************************
                    #
                    # Finally, check for one of the standard component restrictions
                    #
                    if cokey in dCR:
                        resDepth2, resKind = dCR[cokey]

                        if hzDept <= resDepth2 < hzDepb:
                            # This restriction may not be at the top of the horizon, thus we
                            # need to override this if one of the other restrictions exists for this
                            # horizon

                            if len(restriction) == 0:
                                # If this is the only restriction, set the restriction depth
                                # to the value from the corestriction table.
                                resDept = resDepth2

                            # Adding this restriction name to the list even if there are others
                            # May want to take this out later
                            restriction.append(resKind)

                    # ********************************************************
                    #
                    if len(restriction) > 0:
                        # Found at least one restriction for this horizon

                        if not cokey in dComp:
                            # if there are no higher restrictions for this component, save this one
                            # to the dComp dictionary as the top-most restriction
                            #
                            dComp[cokey] = [mukey, compName, localPhase, compPct, resDept, restriction]

                arcpy.SetProgressorPosition()

        arcpy.ResetProgressor()

        # Load restrictions from dComp into dComp2 so that there is complete information for all components

        for cokey in dComp2:
            try:
                dComp2[cokey] = dComp[cokey]

            except:
                pass

        # Return the dictionary containing restriction depths and the dictionary containing defaults
        return dComp2

    except MyError, e:
        # Example: raise MyError("this is an error message")
        PrintMsg(str(e) + " \n", 2)
        return dComp2

    except:
        errorMsg()
        return dComp2


## ===================================================================================
def GetCoRestrictions(outputDB, maxD, resList):
    #
    # Returns a dictionary of top component restrictions for root growth
    #
    # resList is a comma-delimited string of reskind values, surrounded by parenthesis
    #
    # Get component root zone depth from QueryTable_CR and load into dictionary (dCR)
    # This is NOT the final root zone depth. This information will be compared with the
    # horizon soil properties to determine the final root zone depth.

    try:
        rSQL = "resdept_r < " + str(maxD) + " and reskind in " + resList
        sqlClause = (None, "ORDER BY cokey, resdept_r ASC")
        resTbl = os.path.join(outputDB, "QueryTable_CR")
        #PrintMsg("\tGetting corestrictions matching: " + resList, 1)

        if not arcpy.Exists(resTbl):
            raise MyError, "Missing required input table (" + resTbl + ")"

        dRestrictions = dict()

        # Get the top component restriction from the sorted table
        with arcpy.da.SearchCursor(resTbl, ["cokey", "resdept_r", "reskind"], rSQL, "", "", sqlClause) as cur:
            for rec in cur:
                cokey, resDept, reskind = rec
                #PrintMsg("Restriction: " + str(rec), 1)

                if not cokey in dRestrictions:
                    dRestrictions[cokey] = resDept, reskind

        return dRestrictions

    except MyError, e:
        # Example: raise MyError("this is an error message")
        PrintMsg(str(e) + " \n", 2)
        return dict()

    except:
        errorMsg()
        return dict()


## ===================================================================================
def CalcRZAWS(inputDB, outputDB, td, bd, theCompTable, theMuTable, dRestrictions, maxD, dPct):
    # Create a component-level summary table
    # Calculate mapunit-weighted average for each mapunit and write to a mapunit-level table
    # Need to filter out compkind = 'Miscellaneous area' for RZAWS
    # dRestrictions[cokey] = [mukey, compName, localPhase, compPct, resDept, restriction]

    try:
        import decimal

        env.workspace = outputDB

        # Using the same component horizon table that has been
        queryTbl = os.path.join(outputDB, "QueryTable_Hz")

        numRows = int(arcpy.GetCount_management(queryTbl).getOutput(0))

        PrintMsg(" \nCalculating Root Zone AWS for " + str(td) + " to " + str(bd) + "cm...", 0)

        # QueryTable_HZ fields
        qFieldNames = ["mukey", "cokey", "comppct_r",  "compname", "localphase", "majcompflag", "compkind", "taxorder", "taxsubgrp", "desgnmaster", "om_r", "awc_r", "hzdept_r", "hzdepb_r", "texture", "lieutex"]

        #arcpy.SetProgressorLabel("Creating output tables using dominant component...")
        #arcpy.SetProgressor("step", "Calculating root zone available water supply..." , 0, numRows, 1)

        # Open edit session on geodatabase to allow multiple update cursors
        with arcpy.da.Editor(inputDB) as edit:

            # initialize list of components with horizon overlaps
            #badCo = list()

            # Output fields for root zone and droughty
            muFieldNames = ["mukey", "pctearthmc", "rootznemc", "rootznaws", "droughty"]
            muCursor = arcpy.da.UpdateCursor(theMuTable, muFieldNames)

            # Open component-level output table for updates
            #coCursor = arcpy.da.InsertCursor(theCompTable, coFieldNames)
            coFieldNames = ["mukey", "cokey", "compname", "localphase", "comppct_r", "pctearthmc", "rootznemc", "rootznaws", "restriction"]
            coCursor = arcpy.da.UpdateCursor(theCompTable, coFieldNames)

            # Process query table using cursor, write out horizon data for each major component
            postFix = "order by mukey, comppct_r DESC, cokey, hzdept_r ASC"
            iCnt = int(arcpy.GetCount_management(queryTbl).getOutput(0))

            # For root zone calculations, we only want earthy, major components
            #PrintMsg(" \nFiltering components in Query_HZ for CalcRZAWS1 function", 1)
            #
            # Major-Earthy Components
            #hzSQL = "component.compkind <> 'Miscellaneous area' and component.compkind is not NULL and component.majcompflag = 'Yes'"
            # All Components
            hzSQL = ""

            inCur = arcpy.da.SearchCursor(queryTbl, qFieldNames, hzSQL, "", False, (None, postFix))

            arcpy.SetProgressor("step", "Reading query table...",  0, iCnt, 1)

            # Create dictionaries to handle the mapunit and component summaries
            dMu = dict()
            dComp = dict()

            # I may have to pull the sum of component percentages out of this function?
            # It seems to work OK for the earthy-major components, but will not work for
            # the standard AWS calculations. Those 'Miscellaneous area' components with no horizon data
            # are excluded from the Query table because it does not support Outer Joins.
            #
            mCnt = 0
            #PrintMsg("\tmukey, cokey, comppct, top, bottom, resdepth, thickness, aws", 0)

            # TEST: keep list of cokeys as a way to track the top organic horizons
            skipList = list()

            for rec in inCur:
                # read each horizon-level input record from QueryTable_HZ ...
                #
                mukey, cokey, compPct, compName, localPhase, mjrFlag, cKind, taxorder, taxsubgrp, desgnmaster, om, awc, top, bot, texture, lieutex = rec

                if mjrFlag == "Yes" and cKind != "Miscellaneous area" and cKind is not None:

                    # For major-earthy components
                    # Get restriction information from dictionary

                    # For non-Miscellaneous areas with no horizon data, set hzdepth values to zero so that
                    # PWSL and Droughty will get populated with zeros instead of NULL.
                    if top is None and bot is None:

                        if not cokey in dComp:
                            dComp[cokey] = mukey, compName, localPhase, compPct, 0, 0, ""

                    try:
                        # mukey, compName, localPhase, compPct, resDept, restriction
                        # rDepth is the component restriction depth or calculated horizon restriction from CalcRZDepth1 function

                        # mukey, compName, localPhase, compPct, resDept, restriction] = dRestrictions
                        d1, d2, d3, d4, rDepth, restriction = dRestrictions[cokey]
                        cBot = min(rDepth, bot, maxD)  # 01-05-2015 Added maxD because I found 46 CONUS mapunits with a ROOTZNEMC > 150

                        #if mukey == tmukey and rDepth != 150:
                        #    PrintMsg("\tRestriction, " + str(mukey) + ", " + str(cokey) + ", " + str(rDepth) + ", " + str(restriction), 1)

                    except:
                        #errorMsg()
                        cBot = min(maxD, bot)
                        restriction = []
                        rDepth = maxD

                        #if mukey == tmukey:
                        #    PrintMsg("RestrictionError, " + str(mukey) + ", " + str(cokey) + ", " + str(rDepth) + ", " + str(restriction), 1)

                    bOrganic = CheckTexture(mukey, cokey, desgnmaster, om, texture, lieutex, taxorder, taxsubgrp)

                    #if mukey == tmukey and bOrganic:
                    #    PrintMsg("Organic: " + str(mukey) + ", " + str(cokey) )


                    # fix awc_r to 2 decimal places
                    if awc is None:
                        awc = 0.0

                    else:
                        awc = round(awc, 2)

                    # Reasons for skipping RZ calculations on a horizon:
                    #   1. Desgnmaster = O, L and Taxorder != Histosol and is at the surface
                    #   2. Do I need to convert null awc values to zero?
                    #   3. Below component restriction or horizon restriction level

                    if bOrganic and not cokey in skipList:
                        # Organic surface horizon - Not using this horizon in the calculations
                        useHz = False

                        #if mukey == tmukey:
                        #    PrintMsg("Organic, " + str(mukey) + ", " + str(cokey) + ", " + str(compPct) + ", " + str(desgnmaster) + ", " + taxorder  + ", " + str(top) + ", " + str(bot) + ", " + str(cBot)  + ", " + str(awc) + ", " + str(useHz), 1)

                    else:
                        # Mineral, Histosol, buried Organic, Bedrock or there is a horizon restriction (EC, pH - Using this horizon in the calculations
                        useHz = True
                        skipList.append(cokey)

                        # Looking for problems
                        #if mukey == tmukey:
                        #    PrintMsg("Mineral, " + str(mukey) + ", " + str(cokey)  + ", " + str(compPct) + ", " + str(desgnmaster) + ", " + str(taxorder) + ", " + str(top) + ", " + str(bot) + ", " + str(cBot) + ", " + str(awc)  + ", " + str(useHz), 1)

                        # Attempt to fix component with a surface-level restriction that might be in an urban soil
                        if not cokey in dComp and cBot == 0:
                            dComp[cokey] = mukey, compName, localPhase, compPct, 0, 0, restriction

                            # Looking for problems
                            #if mukey == tmukey:
                            #    PrintMsg("MUKEY2: " + str(mukey) + ", " + str(top) + ", " + str(bot) + ", " + str(cBot) + ", " + str(useHz), 1)

                    if top < cBot and useHz == True:
                        # If the top depth is less than the bottom depth, proceed with the calculation
                        # Calculate sum of horizon thickness and sum of component ratings for all horizons above bottom
                        hzT = cBot - top
                        aws = float(hzT) * float(awc) * 10.0

                        # Looking for problems
                        #if mukey == tmukey:
                        #    PrintMsg("MUKEY3: " + str(mukey) + ", " + str(top) + ", " + str(bot) + ", " + str(cBot) + ", " + str(useHz), 1)


                        if cokey in dComp:
                            # accumulate total thickness and total rating value by adding to existing component values
                            mukey, compName, localPhase, compPct, dHzT, dAWS, restriction = dComp[cokey]
                            dAWS = dAWS + aws
                            dHzT += hzT

                            dComp[cokey] = mukey, compName, localPhase, compPct, dHzT, dAWS, restriction

                            # Looking for problems
                            #if mukey == tmukey:
                            #    PrintMsg("MUKEY4: " + str(mukey) + ", " + str(top) + ", " + str(bot) + ", " + str(cBot) + ", " + str(useHz), 1)

                            #if dHzT > 150.0:
                            #    overlap = dHzT - 150
                                # Found some components where there are overlapping horizons
                                #badCo.append(str(cokey))

                        else:
                            # Create initial entry for this component using the first horizon
                            dComp[cokey] = mukey, compName, localPhase, compPct, hzT, aws, restriction

                            # Looking for problems
                            #if mukey == tmukey:
                            #    PrintMsg("MUKEY5: " + str(mukey) + ", " + str(top) + ", " + str(bot) + ", " + str(cBot) + ", " + str(useHz), 1)

                    else:
                        # Do not include this horizon in the rootzone calculations
                        pass

                else:
                    # Not a major-earthy component, so write out everything BUT rzaws-related data (last values)
                    dComp[cokey] = mukey, compName, localPhase, compPct, None, None, None, None

                arcpy.SetProgressorPosition()

                # end of processing major-earthy components

            arcpy.ResetProgressor()

            # get the total number of major-earthy components from the dictionary count
            iComp = len(dComp)

            # Read through the component-level data and summarize to the mapunit level

            if iComp > 0:
                PrintMsg(" \nSaving component average RZAWS to table... (" + str(iComp) + ")", 0 )
                arcpy.SetProgressor("step", "Saving component data...",  0, iComp, 1)
                iCo = 0 # count component records written to theCompTbl

                for corec in coCursor:
                    mukey, cokey, compName, localPhase, compPct, pctearthmc, rDepth, aws, restrictions = corec

                    try:
                        # get sum of component percent for the mapunit
                        pctearthmc = float(dPct[mukey][1])   # sum of comppct_r for all major components Test 2014-10-07

                        # get rootzone data from dComp
                        mukey1, compName1, localPhase1, compPct1, hzT, awc, restriction = dComp[cokey]

                    except:
                        pctearthmc = 0
                        hzT = None
                        rDepth = None
                        awc = None
                        restriction = []

                    # calculate component percentage adjustment
                    if pctearthmc > 0 and not awc is None:
                        # If there is no data for any of the component horizons, could end up with 0 for
                        # sum of comppct_r

                        adjCompPct = float(compPct) / float(pctearthmc)

                        # adjust the rating value down by the component percentage and by the sum of the usable horizon thickness for this component
                        aws = adjCompPct * float(awc) # component rating

                        if restriction is None:
                            restrictions = ''

                        elif len(restriction) > 0:
                            restrictions = ",".join(restriction)

                        else:
                            restrictions = ''

                        corec = mukey, cokey, compName, localPhase, compPct, pctearthmc, hzT, aws, restrictions

                        coCursor.updateRow(corec)
                        iCo += 1

                        # Weight hzT for ROOTZNEMC by component percent
                        hzT = (float(hzT) * float(compPct) / pctearthmc)

                        if mukey in dMu:
                            val1, val2, val3 = dMu[mukey]
                            dMu[mukey] = pctearthmc, (hzT + val2), (aws + val3)

                        else:
                            # first entry for map unit ratings
                            dMu[mukey] = pctearthmc, hzT, aws

                    else:
                        # Populate component level record for a component with no AWC
                        corec = mukey, cokey, compName, localPhase, compPct, None, None, None, ""
                        coCursor.updateRow(corec)
                        iCo += 1

                    arcpy.SetProgressorPosition()

                arcpy.ResetProgressor()

            else:
                raise MyError, "No component data in dictionary dComp"

            if len(dMu) > 0:
                PrintMsg(" \nSaving map unit average RZAWS to table...(" + str(len(dMu)) + ")", 0 )

            else:
                raise MyError, "No map unit information in dictionary dMu"

            # Save root zone available water supply and droughty soils to output map unit table
            #
            for murec in muCursor:
                mukey, pctearthmc, rootznemc, rootznaws, droughty = murec

                try:
                    rec = dMu[mukey]
                    pct, rootznemc, rootznaws = rec
                    pctearthmc = dPct[mukey][1]

                    if rootznemc > 150.0:
                        # This is a bandaid for components that have horizon problems such
                        # overlapping that causes the calculated total to exceed 150cm.
                        rootznemc = 150.0

                    rootznaws = round(rootznaws, 0)
                    rootznemc = round(rootznemc, 0)

                    if rootznaws > 152:
                        droughty = 0

                    else:
                        droughty = 1

                except:
                    pctearthmc = 0
                    rootznemc = None
                    rootznaws = None

                murec = mukey, pctearthmc, rootznemc, rootznaws, droughty
                muCursor.updateRow(murec)

                #if mukey == tmukey:
                    # values at this point seem to be correct
                #    fldnames = muCursor.fields
                #    PrintMsg(str(fldnames), 1)
                #    PrintMsg(str(murec), 1)

            # Save data issues to permanent files for later review
            #if len(badCo) > 0:
            #    fileCo = os.path.basename(inputDB)[:-4] + "_OverlappingHz.txt"
            #    fileCo = os.path.join(os.path.dirname(inputDB), fileCo)
            #    fh = open(fileCo, "w")
            #    fh.write(inputDB + "\n")
            #    fh.write("Components with overlapping horizons\n\n")
            #    fh.write("COKEY IN ('" + "', '".join(badCo) + "') \n")
            #    fh.close()
            #    PrintMsg(" \nComponents with overlapping horizons (" + Number_Format(len(badCo), 0, True) + ") saved to:\t" + fileCo, 0)

            PrintMsg("", 0)

            return True

    except MyError, e:
        # Example: raise MyError("this is an error message")
        PrintMsg(str(e) + " \n", 2)
        return False

    except:
        errorMsg()
        return False

## ===================================================================================
def CalcAWS(inputDB, outputDB, theCompTable, theMuTable, dPct):
    # Create a component-level summary table
    # Calculate the standard mapunit-weighted available waters supply for each mapunit and
    # add it to the map unit-level table.
    #
    # 12-08 I see that for mukey='2479901' my rating is

    try:
        # Using the same component horizon table that has been
        queryTbl = os.path.join(outputDB, "QueryTable_HZ")
        numRows = int(arcpy.GetCount_management(queryTbl).getOutput(0))

        # mukey, cokey, compPct,val, top, bot
        qFieldNames = ["mukey", "cokey", "comppct_r", "awc_r", "hzdept_r", "hzdepb_r"]

        # Track map units that are missing data
        missingList = list()
        minusList = list()

        PrintMsg(" \nCalculating standard available water supply for:", 0)

        for rng in depthList:
            # Calculating and updating just one AWS column at a time
            #
            td = rng[0]
            bd = rng[1]
            #outputFields = "AWS" + str(td) + "_" + str(bd), "TK" + str(td) + "_" + str(bd) + "A"

            # Open output table Mu...All in write mode
            muFieldNames = ["MUKEY", "MUSUMCPCTA", "AWS" + str(td) + "_" + str(bd), "TK" + str(td) + "_" + str(bd) + "A"]
            coFieldNames = ["COKEY", "AWS" + str(td) + "_" + str(bd), "TK" + str(td) + "_" + str(bd) + "A"]

            # Create dictionaries to handle the mapunit and component summaries
            dMu = dict()
            dComp = dict()
            dSum = dict()     # store sum of comppct_r and total thickness for the component
            dHz = dict()      # Trying a new dictionary that will s


            arcpy.SetProgressorLabel("Creating output tables using dominant component...")
            arcpy.SetProgressor("step", "Aggregating data for the dominant component..." , 0, numRows, 1)

            # Open edit session on geodatabase to allow multiple insert cursors
            with arcpy.da.Editor(inputDB) as edit:

                # Open output mapunit-level table in update mode
                # MUKEY, AWS
                muCursor = arcpy.da.UpdateCursor(theMuTable, muFieldNames)

                # Open output component-level table in write mode
                # MUKEY, AWS
                coCursor = arcpy.da.UpdateCursor(theCompTable, coFieldNames)

                # Process query table using a searchcursor, write out horizon data for each component
                # At this time, almost all components are being used! There is no filter.
                postFix = "order by mukey, comppct_r DESC, cokey, hzdept_r ASC"
                #hzSQL = "compkind is not null and hzdept_r is not null"  # prevent divide-by-zero errors
                hzSQL = "hzdept_r is not null"  # prevent divide-by-zero errors by skipping components with no horizons

                iCnt = int(arcpy.GetCount_management(queryTbl).getOutput(0))
                inCur = arcpy.da.SearchCursor(queryTbl, qFieldNames, hzSQL, "", False, (None, postFix))

                arcpy.SetProgressor("step", "Reading QueryTable_HZ ...",  0, iCnt, 1)

                for rec in inCur:
                    # read each horizon-level input record from the query table ...

                    mukey, cokey, compPct, awc, top, bot = rec

                    if awc is not None:

                        # Calculate sum of horizon thickness and sum of component ratings for all horizons above bottom
                        hzT = min(bot, bd) - max(top, td)   # usable thickness from this horizon

                        if hzT > 0:
                            aws = float(hzT) * float(awc) * 10

                            if not cokey in dComp:
                                # Create initial entry for this component using the first horiozon CHK
                                dComp[cokey] = (mukey, compPct, hzT, aws)

                            else:
                                # accumulate total thickness and total rating value by adding to existing component values  CHK
                                mukey, compName, dHzT, dAWS = dComp[cokey]
                                dAWS = dAWS + aws
                                dHzT = dHzT + hzT
                                dComp[cokey] = (mukey, compPct, dHzT, dAWS)

                    arcpy.SetProgressorPosition()

                # get the total number of major components from the dictionary count
                iComp = len(dComp)

                # Read through the component-level data and summarize to the mapunit level

                if iComp > 0:
                    PrintMsg("\t" + str(td) + " - " + str(bd) + "cm (" + Number_Format(iComp, 0, True) + " components)"  , 0)
                    arcpy.SetProgressor("step", "Saving map unit and component AWS data...",  0, iComp, 1)

                    for corec in coCursor:
                        # get component level data  CHK
                        cokey = corec[0]

                        if cokey in dComp:
                            dRec = dComp[corec[0]]
                            mukey, compPct, hzT, awc = dRec

                            # get sum of component percent for the mapunit  CHK
                            try:
                                # Value[0] is for all components,
                                # Value[1] is just for major-earthy components,
                                # Value[2] is all major components
                                # Value[3] is earthy components
                                sumCompPct = float(dPct[mukey][0])
                                #sumCompPct = float(dPct[mukey][1])

                            except:
                                # set the component percent to zero if it is not found in the
                                # dictionary. This is probably a 'Miscellaneous area' not included in the  CHK
                                # data or it has no horizon information.
                                sumCompPct = 0
                                #missingList.append("'" + mukey + "'")

                            # calculate component percentage adjustment
                            if sumCompPct > 0:
                                # If there is no data for any of the component horizons, could end up with 0 for
                                # sum of comppct_r
                                #PrintMsg(" \nMUKEY " + mukey + " - " + compName + " has zero percent Sum Comppct", 1)


                                #adjCompPct = float(compPct) / sumCompPct   # WSS method
                                adjCompPct = compPct / 100.0                # VALU table method

                                # adjust the rating value down by the component percentage and by the sum of the usable horizon thickness for this component
                                aws = round((adjCompPct * awc), 2) # component rating

                                corec[1] = aws
                                hzT = hzT * adjCompPct    # Adjust component share of horizon thickness by comppct
                                corec[2] = hzT             # This is new for the TK0_5A column
                                coCursor.updateRow(corec)

                                # Update component values in component dictionary   CHK
                                # Not sure what dComp is being used for ???
                                dComp[cokey] = mukey, compPct, hzT, aws

                                # Try to fix high mapunit aggregate HZ by weighting with comppct

                                # Testing new mapunit aggregation 09-08-2014
                                # Trying to replace dMu dictionary
                                if mukey in dMu:
                                    val1, val2, val3 = dMu[mukey]
                                    #dMu[mukey] = (compPct + val1, hzT + val2, aws + val2)
                                    compPct = compPct + val1
                                    hzT = hzT + val2
                                    aws = aws + val3

                                #else:
                                dMu[mukey] = (compPct, hzT, aws)
                                #if mukey == '2479892':


                        arcpy.SetProgressorPosition()

                    arcpy.ResetProgressor()

                else:
                    PrintMsg("\t" + Number_Format(iComp, 0, True) + " components for "  + str(td) + " - " + str(bd) + "cm", 1)

                # Write out map unit aggregated AWS
                #
                for murec in muCursor:
                    mukey = murec[0]

                    if mukey in dMu:
                        compPct, hzT, aws = dMu[mukey]
                        murec[1] = compPct
                        murec[2] = aws
                        murec[3] = round(hzT, 2)  # sometimes this ends up being 2 or 3X what it should
                        muCursor.updateRow(murec)

        if len(missingList) > 0:
            missingList = list(set(missingList))
            PrintMsg(" \nFollowing mapunits have no comppct_r: " + ", ".join(missingList), 1)

        PrintMsg("", 0)

        return True

    except MyError, e:
        # Example: raise MyError("this is an error message")
        PrintMsg(str(e) + " \n", 2)
        return False

    except:
        errorMsg()
        return False

## ===================================================================================
def CalcSOC(inputDB, outputDB, theCompTable, theMuTable, dPct, dFrags, depthList, dRestrictions, maxD):
    # Create a component-level summary table
    # Calculate the standard mapunit-weighted available waters supply for each mapunit and
    # add it to the map unit-level table
    #
    # Current problems: have  MuSOC with multiple records for mukey=2479888

    try:
        # Using the same component horizon table that has been
        queryTbl = os.path.join(outputDB, "QueryTable_HZ")
        numRows = int(arcpy.GetCount_management(queryTbl).getOutput(0))

        # mukey, cokey, compPct,val, top, bot
        #qFieldNames = ["mukey", "cokey", "comppct_r", "hzdept_r", "hzdepb_r", "om_r", "dbthirdbar_r"]
        qFieldNames = ["mukey","cokey","comppct_r","compname","localphase","chkey","om_r","dbthirdbar_r", "hzdept_r","hzdepb_r"]

        # Track map units that are missing data
        missingList = list()
        minusList = list()

        PrintMsg(" \nCalculating soil organic carbon for:", 0)

        for rng in depthList:
            # Calculating and updating just one SOC column at a time
            #
            td = rng[0]
            bd = rng[1]

            # Open output table Mu...All in write mode
            # I lately added the "MUSUMCPCTS" to the output. Need to check output because
            # it will be writing this out for every range. Lots more overhead.
            #

            # Create dictionaries to handle the mapunit and component summaries
            dMu = dict()
            dComp = dict()
            #dSumPct = dict()  # store the sum of comppct_r for each mapunit to use in the calculations
            dSum = dict()     # store sum of comppct_r and total thickness for the component
            dHz = dict()      # Trying a new dictionary that will s
            mCnt = 0

            arcpy.SetProgressorLabel("Creating output tables using dominant component...")
            arcpy.SetProgressor("step", "Aggregating data for the dominant component..." , 0, numRows, 1)

            # Open edit session on geodatabase to allow multiple insert cursors
            with arcpy.da.Editor(inputDB) as edit:

                # Open output mapunit-level table in update mode
                muFieldNames = ["MUKEY", "MUSUMCPCTS", "SOC" + str(td) + "_" + str(bd), "TK" + str(td) + "_" + str(bd) + "S"]
                muCursor = arcpy.da.UpdateCursor(theMuTable, muFieldNames)

                # Open output component-level table in write mode

                coFieldNames = ["COKEY", "SOC" + str(td) + "_" + str(bd), "TK" + str(td) + "_" + str(bd) + "S"]
                coCursor = arcpy.da.UpdateCursor(theCompTable, coFieldNames)

                # Process query table using a searchcursor, write out horizon data for each component
                # At this time, almost all components are being used! There is no filter.
                postFix = "order by mukey, comppct_r DESC, cokey, hzdept_r ASC"
                #hzSQL = "compkind is not null and hzdept_r is not null"  # prevent divide-by-zero errors
                hzSQL = "hzdept_r is not null"  # prevent divide-by-zero errors by skipping components with no horizons

                iCnt = int(arcpy.GetCount_management(queryTbl).getOutput(0))
                inCur = arcpy.da.SearchCursor(queryTbl, qFieldNames, hzSQL, "", False, (None, postFix))

                arcpy.SetProgressor("step", "Reading QueryTable_HZ ...",  0, iCnt, 1)

                for rec in inCur:
                    # read each horizon-level input record from the query table ...

                    mukey, cokey, compPct, compName, localPhase, chkey, om, db3, top, bot = rec

                    if om is not None and db3 is not None:
                        # Calculate sum of horizon thickness and sum of component ratings for
                        # that portion of the horizon that is with in the td-bd range
                        top = max(top, td)
                        bot = min(bot, bd)
                        om = round(om, 3)

                        try:
                            rz, resKind = dRestrictions[cokey]

                        except:
                            rz = maxD
                            resKind = ""

                        # Now check for horizon restrictions within this range
                        if top < rz < bot:
                            # restriction found in this horizon, use it to set a new depth
                            #PrintMsg("\t\t" + resKind + " restriction for " + mukey + ":" + cokey + " at " + str(rz) + "cm", 1)
                            cBot = rz

                        else:
                            cBot = min(rz, bot)

                        # Calculate initial usable horizon thickness
                        hzT = cBot - top

                        if hzT > 0 and top < cBot:
                            # get horizon fragment volume
                            try:
                                fragvol = dFrags[chkey]

                            except:
                                fragvol = 0.0
                                pass

                            # Calculate SOC using horizon thickness, OM, BD, FragVol, CompPct
                            soc =  ( (hzT * ( ( om * 0.58 ) * db3 )) / 100.0 ) * ((100.0 - fragvol) / 100.0) * ( compPct * 100 )

                            if not cokey in dComp:
                                # Create initial entry for this component using the first horizon CHK
                                dComp[cokey] = (mukey, compPct, hzT, soc)

                            else:
                                # accumulate total thickness and total rating value by adding to existing component values  CHK
                                mukey, compName, dHzT, dSOC = dComp[cokey]
                                dSOC = dSOC + soc
                                dHzT = dHzT + hzT
                                dComp[cokey] = (mukey, compPct, dHzT, dSOC)


                    arcpy.SetProgressorPosition()

                # get the total number of major components from the dictionary count
                iComp = len(dComp)

                # Read through the component-level data and summarize to the mapunit level
                #
                if iComp > 0:
                    PrintMsg("\t" + str(td) + " - " + str(bd) + "cm (" + Number_Format(iComp, 0, True) + " components)", 0)
                    arcpy.SetProgressor("step", "Saving map unit and component SOC data...",  0, iComp, 1)

                    for corec in coCursor:
                        # Could this be where I am losing minor components????
                        #
                        # get component level data  CHK
                        cokey = corec[0]

                        if cokey in dComp:
                            # get SOC-related data from dComp by cokey
                            # dComp soc = ( (hzT * ( ( om * 0.58 ) * db3 )) / 100.0 ) * ~
                            # ((100.0 - fragvol) / 100.0) * ( compPct * 100 )
                            mukey, compPct, hzT, soc = dComp[corec[0]]

                            # get sum of component percent for the mapunit (all components???)
                            # Value[0] is for all components,
                            # Value[1] is just for major-earthy components,
                            # Value[2] is all major components
                            # Value[3] is earthy components
                            try:
                                sumCompPct = float(dPct[mukey][0]) # Sum comppct for ALl components
                                # Sum of comppct for Major-earthy components
                                #sumCompPct = float(dPct[mukey][1])

                            except:
                                # set the component percent to zero if it is not found in the
                                # dictionary. This is probably a 'Miscellaneous area' not included in the  CHK
                                # data or it has no horizon information.
                                sumCompPct = 0.0
                                #missingList.append("'" + mukey + "'")

                            # calculate component percentage adjustment
                            if sumCompPct > 0:
                                # adjust the rating value down by the component percentage and by the sum of the usable horizon thickness for this component
                                #soc = round((adjCompPct * om), 2) # component rating

                                #soc = round((100.0 * om / sumCompPct), 0)  # Try adjusting soc down by the sum of the component percents
                                adjCompPct = float(compPct) / sumCompPct  #

                                # write the new component-level SOC data to the Co_VALU table

                                corec[1] = soc                      # Test



                                hzT = hzT * compPct / 100.0      # Adjust component share of horizon thickness by comppct
                                corec[2] = hzT             # This is new for the TK0_5A column
                                coCursor.updateRow(corec)

                                # Update component values in component dictionary   CHK
                                dComp[cokey] = mukey, compPct, hzT, soc

                                # Try to fix high mapunit aggregate HZ by weighting with comppct

                                # Testing new mapunit aggregation 09-08-2014
                                # Trying to replace dMu dictionary
                                if mukey in dMu:
                                    val1, val2, val3 = dMu[mukey]
                                    #dMu[mukey] = (compPct + val1, hzT + val2, aws + val2)
                                    compPct = compPct + val1

                                    hzT = hzT + val2
                                    soc = soc + val3

                                dMu[mukey] = (compPct, hzT, soc)  # SOC is a little low
                                #dMu[mukey] = (compPct, hzT, ( soc / adjCompPct))


                        arcpy.SetProgressorPosition()

                    arcpy.ResetProgressor()

                else:
                    PrintMsg("\t" + Number_Format(iComp, 0, True) + " components for "  + str(td) + " - " + str(bd) + "cm", 1)

                # Write out map unit aggregated AWS
                #
                for murec in muCursor:
                    mukey = murec[0]

                    if mukey in dMu:
                        compPct, hzT, soc = dMu[mukey]
                        murec[1] = compPct
                        murec[2] = round(soc, 0)
                        murec[3] = round(hzT, 0)  # this value appears to be low sometimes
                        muCursor.updateRow(murec)

        #if len(missingList) > 0:
        #    missingList = list(set(missingList))
        #    PrintMsg(" \nFollowing mapunits have no comppct_r: " + ", ".join(missingList), 1)

        #PrintMsg("", 0)

        return True

    except MyError, e:
        # Example: raise MyError("this is an error message")
        PrintMsg(str(e) + " \n", 2)
        return False

    except:
        errorMsg()
        return False

## ===================================================================================
def GetFragVol(inputDB):
    # Get the horizon summary of rock fragment volume (percent)
    # load sum of comppct_r into a dictionary by chkey. This
    # value will be used to reduce amount of SOC for each horizon
    # If not all horizons are not present in the dictionary, failover to
    # zero for the fragvol value.

    try:

        fragSQL = ""
        fragFlds = ["chkey", "fragvol_r"]

        dFrags = dict()

        with arcpy.da.SearchCursor(os.path.join(inputDB, "chfrags"), fragFlds, fragSQL) as fragCur:
            for rec in fragCur:
                chkey, fragvol = rec

                if chkey in dFrags:
                    # This horizon has a volume already
                    # Get the existing values
                    # limit total fragvol to 100 or we will get negative SOC values where there
                    # are problems with fragvol data
                    val = dFrags[chkey]
                    dFrags[chkey] = min(val + max(fragvol, 0), 100)

                else:
                    # this is the first component for this map unit
                    dFrags[chkey] = min(max(fragvol, 0), 100)

        # in the rare case where fragvol sum is greater than 100%, return 100
        return dFrags

    except MyError, e:
        # Example: raise MyError("this is an error message")
        PrintMsg(str(e) + " \n", 2)
        return dict()

    except:
        errorMsg()
        return dict()

## ===================================================================================
def GetSumPct(inputDB):
    # Get map unit - sum of component percent for all components and also for major-earthy components
    # load sum of comppct_r into a dictionary.
    # Value[0] is for all components,
    # Value[1] is just for major-earthy components,
    # Value[2] is all major components
    # Value[3] is earthy components
    #
    # Do I need to add another option for earthy components?
    # WSS and SDV use all components with data for AWS.

    try:
        pctSQL = "comppct_r is not null"
        pctFlds = ["mukey", "compkind", "majcompflag", "comppct_r"]

        dPct = dict()

        with arcpy.da.SearchCursor(os.path.join(inputDB, "component"), pctFlds, pctSQL) as pctCur:
            for rec in pctCur:
                mukey, compkind, flag, comppct = rec
                m = 0     # major component percent
                me = 0    # major-earthy component percent
                e = 0     # earthy component percent

                if flag == 'Yes':
                    # major component percent
                    m = comppct

                    if not compkind in  ["Miscellaneous area", ""]:
                        # major-earthy component percent
                        me = comppct
                        e = comppct

                    else:
                        me = 0

                elif not compkind in  ["Miscellaneous area", ""]:
                    e = comppct

                if mukey in dPct:
                    # This mapunit has a pair of values already
                    # Get the existing values from the dictionary
                    #pctAll, pctMjr = dPct[mukey] # all components, major-earthy
                    pctAll, pctME, pctMjr, pctE = dPct[mukey]
                    dPct[mukey] = (pctAll + comppct, pctME + me, pctMjr + m, pctE + e)

                else:
                    # this is the first component for this map unit
                    dPct[mukey] = (comppct, me, m, e)

        return dPct

    except MyError, e:
        # Example: raise MyError("this is an error message")
        PrintMsg(str(e) + " \n", 2)
        return dict()

    except:
        errorMsg()
        return dict()

## ===================================================================================
def MakeNCCPIQueryTable(inputDB, qTable):
    # create query table containing information from component and chorizon tables
    # return name of querytable. Failure returns an empty string for the table name.
    #
    # COINTERP.RULENAME CHOICES:
    #
    # 'NCCPI - National Commodity Crop Productivity Index (Ver 2.0)'
    # 'NCCPI - NCCPI Corn and Soybeans Submodel (II)'
    # 'NCCPI - NCCPI Cotton Submodel (II)'
    # 'NCCPI - NCCPI Small Grains Submodel (II)'
    #

    try:

        # Join chorizon table with component table
        inTables = [os.path.join(inputDB, "component"), os.path.join(inputDB, "cointerp")]

        # interphr is the fuzzy value
        theFields = [["COMPONENT.MUKEY", "MUKEY"], \
        ["COMPONENT.COKEY", "COKEY"], \
        ["COMPONENT.COMPPCT_R", "COMPPCT_R"], \
        ["COINTERP.RULENAME", "RULENAME"], \
        ["COINTERP.RULEDEPTH", "RULEDEPTH"], \
        ["COINTERP.INTERPHR", "INTERPHR"]]

        rule = 'NCCPI - National Commodity Crop Productivity Index (Ver 2.0)'
        theSQL = "COMPONENT.COMPPCT_R > 0 AND COMPONENT.MAJCOMPFLAG = 'Yes' AND COMPONENT.COKEY = COINTERP.COKEY  AND COINTERP.MRULENAME = '" + rule + "'"
        PrintMsg(" \nCalculating NCCPI weighted averages for all major components...", 0)

        # Things to be aware of with MakeQueryTable:
        # USE_KEY_FIELDS does not create OBJECTID field. Lack of OBJECTID precludes sorting on Mukey.
        # ADD_VIRTUAL_KEY_FIELD creates OBJECTID, but qualifies field names using underscore (eg. COMPONENT_COKEY)
        #
        arcpy.MakeQueryTable_management(inTables, qTable, "ADD_VIRTUAL_KEY_FIELD","",theFields, theSQL)

        if arcpy.Exists(qTable):
            iCnt = int(arcpy.GetCount_management(qTable).getOutput(0))
            if iCnt == 0:
                raise MyError, "Failed to retrieve NCCPI data"

        return True

    except MyError, e:
        # Example: raise MyError("this is an error message")
        PrintMsg(str(e) + " \n", 2)
        return False

    except:
        errorMsg()
        return False

## ===================================================================================
def CalcNCCPI(inputDB, qTable, dPct):
    #
    #
    try:
        # and write to Mu_NCCPI2 table
        #
        PrintMsg(" \nAggregating data to mapunit level...", 0)

        # Alternate component fields for all NCCPI values
        cFlds = ["MUKEY","COKEY","COMPPCT_R","COMPNAME","LOCALPHASE", "DUMMY"]
        mFlds = ["MUKEY","COMPPCT_R","COMPNAME","LOCALPHASE", "DUMMY"]

        # Create dictionary key as MUKEY:INTERPHRC
        # Need to look through the component rating class for ruledepth = 0
        # and sum up COMPPCT_R for each key value
        #
        dVals = dict()  # dictionary containing sum of comppct for each MUKEY:RATING combination

        # Get sum of component percent for each map unit. There are different options:
        #     1. Use all major components
        #     2. Use all components that have an NCCPI rating
        #     3. Use all major-earthy components. This one is not currently available.
        #     4. Use all components (that have a component percent)
        #

        # Query table fields
        qFields = ["COMPONENT_MUKEY", "COMPONENT_COKEY", "COMPONENT_COMPPCT_R", "COINTERP_RULEDEPTH", "COINTERP_RULENAME", "COINTERP_INTERPHR"]

        sortFields = "ORDER BY COMPONENT_COKEY ASC, COMPONENT_COMPPCT_R DESC"
        querytblSQL = "COMPONENT_COMPPCT_R IS NOT NULL"

        iCnt = int(arcpy.GetCount_management(qTable).getOutput(0))
        noVal = list()  # Get a list of components with no overall index rating

        PrintMsg(" \nReading query table with " + Number_Format(iCnt, 0, True) + " records...", 0)

        arcpy.SetProgressor("step", "Reading query table...", 0,iCnt, 1)

        with arcpy.da.SearchCursor(qTable, qFields, querytblSQL, "","", (None, sortFields)) as qCursor:

            for qRec in qCursor:
                # qFields = MUKEY, COKEY, COMPPCT_R, COMPNAME, LOCALPHASE, RULEDEPTH, RULENAME, INTERPHR
                mukey, cokey, comppct, ruleDepth, ruleName, fuzzyValue = qRec

                # Dictionary order:  All, CS, CT, SG
                if not mukey in dVals:
                    # Initialize mukey NCCPI values
                    dVals[mukey] = [None, None, None, None]

                if not fuzzyValue is None:

                    if ruleDepth == 0:
                        # This is NCCPI Overall Index
                        oldVal = dVals[mukey][0]

                        if oldVal is None:
                            dVals[mukey][0] = fuzzyValue * comppct

                        else:
                            dVals[mukey][0] = (oldVal + (fuzzyValue * comppct))

                    elif ruleName == "NCCPI - NCCPI Corn and Soybeans Submodel (II)":
                        oldVal = dVals[mukey][1]

                        if oldVal is None:
                            dVals[mukey][1] = fuzzyValue * comppct

                        else:
                            dVals[mukey][1] = (oldVal + (fuzzyValue * comppct))

                    elif ruleName == "NCCPI - NCCPI Cotton Submodel (II)":
                        oldVal = dVals[mukey][2]

                        if oldVal is None:
                            dVals[mukey][2] =  fuzzyValue * comppct

                        else:
                            dVals[mukey][2] = (oldVal + (fuzzyValue * comppct))

                    elif ruleName == "NCCPI - NCCPI Small Grains Submodel (II)":
                        oldVal = dVals[mukey][3]

                        if oldVal is None:
                            dVals[mukey][3] = fuzzyValue * comppct

                        else:
                            dVals[mukey][3] = (oldVal + (fuzzyValue * comppct))

                elif ruleName == "NCCPI - National Commodity Crop Productivity Index (Ver 2.0)":
                    # This component does not have an NCCPI rating
                    #PrintMsg(" \n" + mukey + ":" + cokey + ", " + str(comppct) + "% has no NCCPI rating", 1)
                    noVal.append("'" + cokey + "'")

                arcpy.SetProgressorPosition()
                #
                # End of query table iteration
                #
        #if len(noVal) > 0:
        #    PrintMsg(" \nThe following components had no NCCPI overall index: " + ", ".join(noVal), 1)

        iCnt = len(dVals)

        if iCnt > 0:

            PrintMsg(" \nSaving map unit weighted NCCPI data (" + Number_Format(iCnt, 0, True) + " records) to " + os.path.basename(theMuTable) + "..." , 0)
            # Write map unit aggregate data to Mu_NCCPI2 table
            #
            # theMuTable is a global variable. Need to check this out in the gSSURGO_ValuTable script

            with arcpy.da.UpdateCursor(theMuTable, ["mukey", "NCCPI2CS", "NCCPI2CO","NCCPI2SG", "NCCPI2ALL"]) as muCur:

                arcpy.SetProgressor("step", "Saving map unit weighted NCCPI data to VALU table...", 0, iCnt, 0)
                for rec in muCur:
                    mukey = rec[0]

                    try:
                        # Get output values from dVals and dPct dictionaries
                        val = dVals[mukey]
                        ovrall, cs, co, sg = val
                        sumPct = dPct[mukey][2]  # sum of major-earthy components
                        rec = mukey, round(cs / sumPct, 2), round(co / sumPct, 2), round(sg / sumPct, 2), round(ovrall / sumPct, 2)
                        muCur.updateRow(rec)
                        #PrintMsg(mukey + ", " + str(round(cs / sumPct, 2)) + ", " + str(round(co / sumPct, 2)) + ", " + str(round(sg / sumPct, 2)) + ", " + str(round(ovrall / sumPct, 2)), 1)

                    except:
                        # Miscellaneous map unit encountered with no comppct_r?
                        pass

                    arcpy.SetProgressorPosition()

            arcpy.Delete_management(qTable)
            return True

        else:
            raise MyError, "No NCCPI data processed"

    except MyError, e:
        # Example: raise MyError("this is an error message")
        PrintMsg(str(e) + " \n", 2)
        return False

    except:
        errorMsg()
        return False

## ===================================================================================
def CalcPWSL(inputDB, outputDB, theMuTable, dPct):
    # Get potential wet soil landscape rating for each map unit
    # Assuming that all components (with comppct_r) will be processed
    #
    # Sharon: I treat all map unit components the same, so if I find 1% water I think
    # it should show up as 1% PWSL.  If the percentage of water is >= 80% then I class
    # it into the water body category or 999.

    # Steve: In looking at some of Sharon's PWSL=999, I find some mapunits where the
    # COMPNAME='Water' at 100%, and HYDRIC='Unranked'. My Valu2 table PWSL is null. That isn't right.

    try:
        env.workspace = outputDB

        # Using the same component horizon table as always
        queryTbl = os.path.join(outputDB, "QueryTable_Hz")
        numRows = int(arcpy.GetCount_management(queryTbl).getOutput(0))
        PrintMsg(" \nCalculating Potential Wet Soil Landscapes using " + queryTbl + "...", 0)
        qFieldNames = ["mukey", "muname", "cokey", "comppct_r",  "compname", "localphase", "otherph", "majcompflag", "compkind", "hydricrating", "drainagecl"]

        pwSQL = ""
        compList = list()
        dMu = dict()

        drainList = ["Poorly drained", "Very poorly drained"]
        phaseList = ["drained", "undrained", "channeled", "protected", "ponded", "flooded"]
        iCnt = int(arcpy.GetCount_management(queryTbl).getOutput(0))
        lastCokey = 'xxx'
        arcpy.SetProgressor("step", "Reading query table table for wetland information...",  0, iCnt, 1)

        with arcpy.da.SearchCursor(queryTbl, qFieldNames, pwSQL) as pwCur:
            for rec in pwCur:
                mukey, muname, cokey, comppct_r,  compname, localphase, otherph, majcompflag, compkind, hydricrating, drainagecl = rec

                #if mukey == tmukey:
                #    PrintMsg(" \nGetting PWSL data for test map unit: " + str(mukey), 1)

                if cokey != lastCokey:
                    # only process first horizon record for each component

                    compList.append(cokey)
                    # Only check the first horizon record, really only need component level
                    # Not very efficient, should problably create a new query table
                    #
                    # Need to split up these tests so that None types can be handled

                    # Sharon says that if the hydricrating for a component is 'No', don't
                    # look at it any further. If it is unranked, go ahead and look at
                    # other properties.
                    #

                    pw = False

                    if muname == 'Water' or compname == 'Water':
                        # Check for water before looking at Hydric rating
                        # Probably won't catch everything. Waiting for Sharon's criteria.

                        if comppct_r >= 80:
                            # Flag this mapunit with a '999'
                            # Not necessarily catching map unit with more than one Water component that
                            # might sum to >= 80.
                            pw = False
                            dMu[mukey] = 999

                        else:
                            pw = True

                            try:
                                sumPct = dMu[mukey]

                                if sumPct != 999:
                                    dMu[mukey] = sumPct + comppct_r

                            except:
                                dMu[mukey] = comppct_r

                    elif hydricrating == 'No':
                        # Added this bit so that other properties cannot override hydricrating = 'No'
                        pw = False

                    elif hydricrating == 'Yes':
                        # This is always a Hydric component
                        # Get component percent and add to map unit total PWSL
                        pw = True
                        #if mukey == tmukey:
                        #    PrintMsg("\tHydric percent = " + str(comppct_r), 1)

                        try:
                            sumPct = dMu[mukey]

                            if sumPct != 999:
                                dMu[mukey] = sumPct + comppct_r

                        except:
                            dMu[mukey] = comppct_r

                    elif hydricrating == 'Unranked':
                        # Not sure how Sharon is handling NULL hydric
                        #
                        # Unranked hydric from here on down, looking at other properties such as:
                        #   Local phase
                        #   Other phase
                        #   Drainage class
                        #   Map unit name strings

                        #if localphase is not None:
                        if [d for d in phaseList if str(localphase).lower().find(d) >= 0]:
                            pw = True
                            #if mukey == tmukey:
                            #    PrintMsg("\tLocalPhase:" + localphase, 1)

                            try:
                                sumPct = dMu[mukey]
                                dMu[mukey] = sumPct + comppct_r

                            except:
                                dMu[mukey] = comppct_r


                        # otherphase
                        #if pw == False and otherph is not None:
                        elif [d for d in phaseList if str(otherph).lower().find(d) >= 0]:
                            pw = True
                            #if mukey == tmukey:
                            #    PrintMsg("\tOtherPhase:" + localphase, 1)

                            try:
                                sumPct = dMu[mukey]
                                dMu[mukey] = sumPct + comppct_r

                            except:
                                dMu[mukey] = comppct_r

                        # look for specific strings in the map unit name
                        #if pw == False and muname is not None:
                        elif [d for d in phaseList if muname.find(d) >= 0]:
                            pw = True
                            #if mukey == tmukey:
                            #    PrintMsg("\tMuname = " + muname, 1)

                            try:
                                sumPct = dMu[mukey]
                                dMu[mukey] = sumPct + comppct_r

                            except:
                                dMu[mukey] = comppct_r

                        elif str(drainagecl) in drainList:
                            pw = True
                            #if mukey == tmukey:
                            #    PrintMsg("\tDrainageClass = " + drainagecl, 1)

                            try:
                                sumPct = dMu[mukey]
                                dMu[mukey] = sumPct + comppct_r

                            except:
                                dMu[mukey] = comppct_r

                lastCokey = cokey # use this to skip the rest of the horizons for this component
                arcpy.SetProgressorPosition()

        if len(dMu) > 0:
            arcpy.SetProgressor("step", "Populating " + os.path.basename(theMuTable) + "...",  0, len(dMu), 1)

            # Populate the PWSL1POMU column in the map unit level table
            muFlds = ["mukey", "pwsl1pomu"]
            with arcpy.da.UpdateCursor(theMuTable, muFlds) as muCur:
                for rec in muCur:
                    mukey = rec[0]
                    try:
                        rec[1] = dMu[mukey]
                        muCur.updateRow(rec)

                    except:
                        pass

                    arcpy.SetProgressorPosition()

        arcpy.ResetProgressor()
        return True

    except MyError, e:
        # Example: raise MyError("this is an error message")
        PrintMsg(str(e) + " \n", 2)
        return False

    except:
        errorMsg()
        return False

## ===================================================================================
def StateNames():
    # Create dictionary object containing list of state abbreviations and their names that
    # will be used to name the file geodatabase.
    # For some areas such as Puerto Rico, U.S. Virgin Islands, Pacific Islands Area the
    # abbrevation is

    # NEED TO UPDATE THIS FUNCTION TO USE THE LAOVERLAP TABLE AREANAME. AREASYMBOL IS STATE ABBREV

    try:
        stDict = dict()
        stDict["AL"] = "Alabama"
        stDict["AK"] = "Alaska"
        stDict["AS"] = "American Samoa"
        stDict["AZ"] = "Arizona"
        stDict["AR"] = "Arkansas"
        stDict["CA"] = "California"
        stDict["CO"] = "Colorado"
        stDict["CT"] = "Connecticut"
        stDict["DC"] = "District of Columbia"
        stDict["DE"] = "Delaware"
        stDict["FL"] = "Florida"
        stDict["GA"] = "Georgia"
        stDict["HI"] = "Hawaii"
        stDict["ID"] = "Idaho"
        stDict["IL"] = "Illinois"
        stDict["IN"] = "Indiana"
        stDict["IA"] = "Iowa"
        stDict["KS"] = "Kansas"
        stDict["KY"] = "Kentucky"
        stDict["LA"] = "Louisiana"
        stDict["ME"] = "Maine"
        stDict["MD"] = "Maryland"
        stDict["MA"] = "Massachusetts"
        stDict["MI"] = "Michigan"
        stDict["MN"] = "Minnesota"
        stDict["MS"] = "Mississippi"
        stDict["MO"] = "Missouri"
        stDict["MT"] = "Montana"
        stDict["NE"] = "Nebraska"
        stDict["NV"] = "Nevada"
        stDict["NH"] = "New Hampshire"
        stDict["NJ"] = "New Jersey"
        stDict["NM"] = "New Mexico"
        stDict["NY"] = "New York"
        stDict["NC"] = "North Carolina"
        stDict["ND"] = "North Dakota"
        stDict["OH"] = "Ohio"
        stDict["OK"] = "Oklahoma"
        stDict["OR"] = "Oregon"
        stDict["PA"] = "Pennsylvania"
        stDict["PRUSVI"] = "Puerto Rico and U.S. Virgin Islands"
        stDict["RI"] = "Rhode Island"
        stDict["Sc"] = "South Carolina"
        stDict["SD"] ="South Dakota"
        stDict["TN"] = "Tennessee"
        stDict["TX"] = "Texas"
        stDict["UT"] = "Utah"
        stDict["VT"] = "Vermont"
        stDict["VA"] = "Virginia"
        stDict["WA"] = "Washington"
        stDict["WV"] = "West Virginia"
        stDict["WI"] = "Wisconsin"
        stDict["WY"] = "Wyoming"
        return stDict

    except:
        PrintMsg("\tFailed to create list of state abbreviations (CreateStateList)", 2)
        return stDict

## ===================================================================================
def UpdateMetadata(outputWS, target, surveyInfo):
    # Update metadata for target object (VALU1 table)
    #
    try:

        # Clear process steps from the VALU1 table. Mostly AddField statements.
        #
        # Different path for ArcGIS 10.2.2??
        #
        #
        if not arcpy.Exists(target):
            target = os.path.join(outputWS, target)

        # Remove geoprocessing history
        remove_gp_history_xslt = os.path.join(os.path.dirname(sys.argv[0]), "remove geoprocessing history.xslt")
        out_xml = os.path.join(env.scratchFolder, "xxClean.xml")

        if arcpy.Exists(out_xml):
            arcpy.Delete_management(out_xml)

        # Using the stylesheet, write 'clean' metadata to out_xml file and then import back in
        arcpy.XSLTransform_conversion(target, remove_gp_history_xslt, out_xml, "")
        arcpy.MetadataImporter_conversion(out_xml, os.path.join(outputWS, target))

        # Set metadata translator file
        dInstall = arcpy.GetInstallInfo()
        installPath = dInstall["InstallDir"]
        prod = r"Metadata/Translator/ARCGIS2FGDC.xml"
        mdTranslator = os.path.join(installPath, prod)

        # Define input and output XML files
        #mdExport = os.path.join(env.scratchFolder, "xxExport.xml")  # initial metadata exported from current data data
        xmlPath = os.path.dirname(sys.argv[0])
        mdExport = os.path.join(xmlPath, "gSSURGO_ValuTable.xml")  # template metadata stored in ArcTool folder
        mdImport = os.path.join(env.scratchFolder, "xxImport.xml")  # the metadata xml that will provide the updated info

        # Cleanup XML files from previous runs
        if os.path.isfile(mdImport):
            os.remove(mdImport)

        # Start editing metadata using search and replace
        #
        stDict = StateNames()
        st = os.path.basename(outputWS)[8:-4]

        if st in stDict:
            # Get state name from the geodatabase name
            mdState = stDict[st]

        else:
            mdState = ""

        # Set date strings for metadata, based upon today's date
        #
        d = datetime.now()
        #today = d.strftime('%Y%m%d')

        # Alternative to using today's date. Use the last SAVEREST date
        today = GetLastDate(outputWS)

        # Set fiscal year according to the current month. If run during January thru September,
        # set it to the current calendar year. Otherwise set it to the next calendar year.
        #
        if d.month > 9:
            fy = "FY" + str(d.year + 1)

        else:
            fy = "FY" + str(d.year)

        # Convert XML from template metadata to tree format
        tree = ET.parse(mdExport)
        root = tree.getroot()

        # new citeInfo has title.text, edition.text, serinfo/issue.text
        citeInfo = root.findall('idinfo/citation/citeinfo/')

        if not citeInfo is None:
            # Process citation elements
            # title
            #
            # PrintMsg("citeInfo with " + str(len(citeInfo)) + " elements : " + str(citeInfo), 1)
            for child in citeInfo:
                PrintMsg("\t\t" + str(child.tag), 0)
                if child.tag == "title":
                    child.text = os.path.basename(target).title()

                    if mdState != "":
                        child.text = child.text + " - " + mdState

                elif child.tag == "edition":
                    if child.text.find('xxFYxx') >= 0:
                        child.text = child.text.replace('xxFYxx', fy)
                    else:
                        PrintMsg(" \n\tEdition: " + child.text, 1)

                    if child.text.find('xxTODAYxx') >= 0:
                        child.text = child.text.replace('xxTODAYxx', today)

                elif child.tag == "serinfo":
                    for subchild in child.iter('issue'):
                        if subchild.text == "xxFYxx":
                            subchild.text = fy

                        if child.text.find('xxTODAYxx') >= 0:
                            child.text = child.text.replace('xxTODAYxx', today)


        # Update place keywords
        #PrintMsg("\tplace keywords", 0)
        ePlace = root.find('idinfo/keywords/theme')

        if ePlace is not None:
            for child in ePlace.iter('themekey'):
                if child.text == "xxSTATExx":
                    #PrintMsg("\tReplaced xxSTATExx with " + mdState)
                    child.text = mdState

                elif child.text == "xxSURVEYSxx":
                    #child.text = "The Survey List"
                    child.text = surveyInfo

        else:
            PrintMsg("\tsearchKeys not found", 1)

        idDescript = root.find('idinfo/descript')

        if not idDescript is None:
            for child in idDescript.iter('supplinf'):
                #id = child.text
                #PrintMsg("\tip: " + ip, 1)
                if child.text.find("xxTODAYxx") >= 0:
                    #PrintMsg("\t\tip", 1)
                    child.text = child.text.replace("xxTODAYxx", today)

                if child.text.find("xxFYxx") >= 0:
                    #PrintMsg("\t\tip", 1)
                    child.text = child.text.replace("xxFYxx", fy)

        if not idDescript is None:
            for child in idDescript.iter('purpose'):
                #ip = child.text
                #PrintMsg("\tip: " + ip, 1)
                if child.text.find("xxFYxx") >= 0:
                    #PrintMsg("\t\tip", 1)
                    child.text = child.text.replace("xxFYxx", fy)

                if child.text.find("xxTODAYxx") >= 0:
                    #PrintMsg("\t\tip", 1)
                    child.text = child.text.replace("xxTODAYxx", today)

        idAbstract = root.find('idinfo/descript/abstract')
        if not idAbstract is None:
            iab = idAbstract.text

            if iab.find("xxFYxx") >= 0:
                #PrintMsg("\t\tip", 1)
                idAbstract.text = iab.replace("xxFYxx", fy)
                #PrintMsg("\tAbstract", 0)

        # Use contraints
        #idConstr = root.find('idinfo/useconst')
        #if not idConstr is None:
        #    iac = idConstr.text
            #PrintMsg("\tip: " + ip, 1)
        #    if iac.find("xxFYxx") >= 0:
        #        idConstr.text = iac.replace("xxFYxx", fy)
        #        PrintMsg("\t\tUse Constraint: " + idConstr.text, 0)

        # Update credits
        eIdInfo = root.find('idinfo')

        if not eIdInfo is None:

            for child in eIdInfo.iter('datacred'):
                sCreds = child.text

                if sCreds.find("xxTODAYxx") >= 0:
                    #PrintMsg("\tdata credits1", 1)
                    sCreds = sCreds.replace("xxTODAYxx", today)

                if sCreds.find("xxFYxx") >= 0:
                    #PrintMsg("\tdata credits2", 1)
                    sCreds = sCreds.replace("xxFYxx", fy)

                child.text = sCreds
                #PrintMsg("\tCredits: " + sCreds, 1)

        #  create new xml file which will be imported, thereby updating the table's metadata
        tree.write(mdImport, encoding="utf-8", xml_declaration=None, default_namespace=None, method="xml")

        # import updated metadata to the geodatabase table
        arcpy.MetadataImporter_conversion(mdExport, target)
        arcpy.ImportMetadata_conversion(mdImport, "FROM_FGDC", target, "DISABLED")

        # delete the temporary xml metadata files
        if os.path.isfile(mdImport):
            os.remove(mdImport)

        #if os.path.isfile(mdExport):
        #    os.remove(mdExport)

        return True

    except:
        errorMsg()
        False

## ===================================================================================
## ====================================== Main Body ==================================
# Import modules
import os, sys, string, re, locale, arcpy, traceback, collections
from operator import itemgetter, attrgetter
import xml.etree.cElementTree as ET
from datetime import datetime
from arcpy import env

# Original input table fields:
# areasymbol, mukey, musym, muname, muname, cokey, compct, compname, compkind, localphase,
# taxorder, taxsubgrp, ec, pH, dbthirdbar, hzname, hzdesgn, hzdept, hzdepb, hzthk, sand, silt,
# clay, om, reskind, reshard, resdept, resthk, texture, lieutex

# Compkind values for major components: Miscellaneous area, Series, Taxon above family,
# Family, Taxadjunct, Variant

# chtexturegrp.texture: 3,468 unique values, too many to list here. Bob's queries only
# look at excluding 'COP-MAT' texture from the Dense Layer calculation. He identifies the rest
# of the organic horizons using lieutext.

# lieutext values: Slightly decomposed plant material, Moderately decomposed plant material,
# Bedrock, Variable, Peat, Material, Unweathered bedrock, Sand and gravel, Mucky peat, Muck,
# Highly decomposed plant material, Weathered bedrock, Cemented, Gravel, Water, Cobbles,
# Stones, Channers, Parachanners, Indurated, Cinders, Duripan, Fragmental material, Paragravel,
# Artifacts, Boulders, Marl, Flagstones, Coprogenous earth, Ashy, Gypsiferous material,
# Petrocalcic, Paracobbles, Diatomaceous earth, Fine gypsum material, Undecomposed organic matter
#
# reskind values: Strongly contrasting textural stratification, Lithic bedrock, Densic material,
# Ortstein, Permafrost, Paralithic bedrock, Cemented horizon, Undefined, Fragipan, Plinthite,
# Abrupt textural change, Natric, Petrocalcic, Duripan, Densic bedrock, Salic,
# Human-manufactured materials, Sulfuric, Placic, Petroferric, Petrogypsic

# Dobos - NASIS lieutext matches: mpm, mpt, muck, peat, spm, udom, pdom, hpm
#
# Paul's organic horizon filters:  chtexturegrp.texture <> 'SPM', <>'UDOM'???, NOT LIKE: '%MPT%', '%MUCK', '%PEAT%' (from SVI query)
#
# Desgnmaster filter:  hzdesn IN ["O", "O'", "L"]
#
# Taxonomic Order filter:  upper(taxorder) LIKE 'HISTOSOLS'. Is this being used?
# Perhaps I should be using taxsubgrp.lower() like '%histic%' or use both!!!

#
# compkind filter for earthy components:  <> 'Miscellaneous area'


try:
    # Diagnostic map unit mukey:
    tmukey = 'xx'

    # Create geoprocessor object
    #gp = arcgisscripting.create(9.3)
    arcpy.OverwriteOutput = True

    inputDB = arcpy.GetParameterAsText(0)            # Input gSSURGO database

    # Set location for temporary tables
    #outputDB = "IN_MEMORY"
    outputDB = env.scratchGDB

    # Name of mapunit level output table (global variable)
    theMuTable = os.path.join(inputDB, "Valu2")

    # Name of component level output table (global variable)
    theCompTable = os.path.join(inputDB, "Co_VALU")

    # Set output workspace to same as the input table
    #env.workspace = os.path.dirname(arcpy.Describe(queryTbl).catalogPath)
    env.workspace = inputDB

    # Save record of any issues to a text file
    logFile = os.path.basename(inputDB)[:-4] + "_Problems.txt"
    logFile = os.path.join(os.path.dirname(inputDB), logFile)

    # Get the mapunit - sum of component percent for calculations
    dPct = GetSumPct(inputDB)
    if len(dPct) == 0:
        raise MyError, ""

    # Create initial set of query tables used for RZAWS, AWS and SOC
    if CreateQueryTables(inputDB, outputDB, 150.0) == False:
        raise MyError, ""

    # Create permanent output tables for the map unit and component levels
    depthList = [(0,5), (5, 20), (20, 50), (50, 100), (100, 150), (150, 999), (0, 20), (0, 30), (0, 100), (0, 150), (0, 999)]

    if CreateOutputTableMu(theMuTable, depthList, dPct) == False:
        raise MyError, ""

    if CreateOutputTableCo(theCompTable, depthList) == False:
        raise MyError, ""

    # Store component restrictions for root growth in a dictionary
    resListAWS = "('Lithic bedrock','Paralithic bedrock','Densic bedrock', 'Densic material', 'Fragipan', 'Duripan', 'Sulfuric')"
    dRZRestrictions = GetCoRestrictions(outputDB, 150.0, resListAWS)

    # Find the top restriction for each component, both from the corestrictions table and the horizon properties
    dComp2 = CalcRZDepth(inputDB, outputDB, theCompTable, theMuTable, 150.0, dPct, dRZRestrictions)

    # Calculate root zone available water capacity using a floor of 150cm or a root restriction depth
    #
    # dComp2[cokey] = [mukey, compName, localPhase, compPct, resDept, restriction]
    if CalcRZAWS(inputDB, outputDB, 0.0, 150.0, theCompTable, theMuTable, dComp2, 150.0, dPct) == False:
        raise MyError, ""

    # Calculate standard available water supply
    if CalcAWS(inputDB, outputDB, theCompTable, theMuTable, dPct) == False:
        raise MyError, ""

    # Run SOC calculations
    # Seems to be a problem with SOC calculations, numbers are high
    maxD = 999.0
    # Get bedrock restrictions for SOC  and write them to the output tables
    resListSOC = "('Lithic bedrock', 'Paralithic bedrock', 'Densic bedrock')"
    dSOCRestrictions = GetCoRestrictions(outputDB, maxD, resListSOC)

    # Store horizon fragment volumes in a dictionary and use in the root zone SOC calculations
    dFrags = GetFragVol(inputDB)

    if len(dFrags) == 0:
        raise MyError, "No fragment volume information"

    # Calculate soil organic carbon for all the different depth ranges
    depthList = [(0,5), (5, 20), (20, 50), (50, 100), (100, 150), (150, 999), (0, 20), (0, 30), (0, 100), (0, 150), (0, 999)]
    if CalcSOC(inputDB, outputDB, theCompTable, theMuTable, dPct, dFrags, depthList, dSOCRestrictions, maxD) == False:
        raise MyError, ""

    # Calculate NCCPI
    # Create query table using component and chorizon tables
    nccpiTbl = "NCCPI_Table"
    if  MakeNCCPIQueryTable(inputDB, nccpiTbl) == False:
        raise MyError, ""

    if CalcNCCPI(inputDB, nccpiTbl, dPct) == False:
        raise MyError, ""

    if CalcPWSL(inputDB, outputDB, theMuTable, dPct) == False:
        raise MyError, ""

    PrintMsg(" \nAll calculations complete", 0)

    # Create metadata for the VALU table
    # Query the output SACATALOG table to get list of surveys that were exported to the gSSURGO
    #
    saTbl = os.path.join(inputDB, "sacatalog")
    expList = list()
    queryList = list()

    with arcpy.da.SearchCursor(saTbl, ["AREASYMBOL", "SAVEREST"]) as srcCursor:
        for rec in srcCursor:
            expList.append(rec[0] + " (" + str(rec[1]).split()[0] + ")")
            queryList.append("'" + rec[0] + "'")

    surveyInfo = ", ".join(expList)
    queryInfo = ", ".join(queryList)

    # Update metadata for the geodatabase and all featureclasses
    PrintMsg(" \nUpdating " + os.path.basename(theMuTable) + " metadata...", 0)
    bMetadata = UpdateMetadata(inputDB, theMuTable, surveyInfo)

    PrintMsg("\tMetadata complete", 0)

    PrintMsg(" \nValu2 table complete for " + inputDB + " \n ", 0)

except MyError, e:
    # Example: raise MyError("this is an error message")
    PrintMsg(str(e) + " \n", 2)

except:
    errorMsg()
