''' ---------------------------------------------------------------------------
Name: Gridded PMP Tool Python Script

Script Version: 5

Python Version: 3.11.8

ArcGIS Version: ArcGIS Pro 3.3

Author: Applied Weather Associates

Usage:  The tool is designed to be executed within an ArcGIS Pro Project.

Description:
    This tool calculates PMP depths for a given drainage basin for the
specified durations.  PMP point values are calculated (in inches) for each
grid point (spaced at 90 arc-second intervals) over the project domain. The
points are converted to gridded PMP datasets for each duration.  Applicable
Storm-based and critically stakced temporal distributions are applied to the PMP.
Gridded ASCII or netCDF4 format are produced for the incremental temporally 
distributed rainfall.
---------------------------------------------------------------------------'''

###########################################################################
## import Python modules

import sys
import arcpy
import os
import traceback
from arcpy import env
import arcpy.analysis as an
import arcpy.management as dm
import arcpy.conversion as con
import numpy as np
import pandas as pd
import netCDF4 as nc
import xarray as xr
import datetime
from pandas import ExcelFile
import matplotlib.pyplot as plt
from heapq import nlargest
import math
from arcpy.sa import *

env.overwriteOutput = True                                                # Set overwrite option
env.addOutputsToMap = False

###########################################################################
## get input parameters

basin = arcpy.GetParameter(0)                                                   # get AOI Basin Shapefile
#home = arcpy.GetParameterAsText(1)                                              # get location of 'PMP' Project Folder
outLocation = arcpy.GetParameterAsText(1)
locDurations = arcpy.GetParameter(2)                                            # get local storm durations (string)
genDurations = arcpy.GetParameter(3)                                            # get general storm durations (string)
tropDurations = arcpy.GetParameter(4)                                           # get tropical storm durations (string)
useBasinArea = arcpy.GetParameter(5)

weightedAve = arcpy.GetParameter(7)						# get option to apply weighted average (boolean)
#outputTable = arcpy.GetParameter(9) 						# get file path for basin average summary table
includeSubbasin = arcpy.GetParameter(8) 					# get option add subbasin averages (boolean)
subbasinIDfield = arcpy.GetParameterAsText(9) 					# Subbasin ID field from AOI Basin Shapefile
ddChart = arcpy.GetParameter(10) 						# get option to create depth-duration chart(boolean)
includeTemporal = arcpy.GetParameter(11)                        # get option to apply temporal patterns
exportASCII = arcpy.GetParameter(12)                                            # get option to Export temporal patterns to ASCII (boolean) 
exportNetCDF = arcpy.GetParameter(13)                                            # get option to Export temporal patterns to ASCII (boolean) 

scriptPath = os.getcwd()
home = os.path.dirname(scriptPath)
arcpy.AddMessage("Path to PMP_Evaluation_Tool....  " + home)


dadGDB = home + "\\Input\\DAD_Tables.gdb"                                       # location of DAD tables
adjFactGDB = home + "\\Input\\Storm_Adj_Factors.gdb"                            # location of feature datasets containing total adjustment factors 
arcpy.AddMessage("\nDAD Tables geodatabase path:  " + dadGDB)
arcpy.AddMessage("Storm Adjustment Factor geodatabase path:  " + adjFactGDB)

basAveTables = []                                                               # global list of Basin Average Summary tables     
gridRes = "0.025"                                                               # global output grid resolution

def pmpAnalysis(aoiBasin, stormType, durList):

    ###########################################################################
    ## Create PMP Point Feature Class from points within AOI basin and add fields
    def createPMPfc():

        arcpy.AddMessage("\nCreating feature class: 'PMP_Points' in Scratch.gdb...")
        dm.MakeFeatureLayer(home + "\\Input\\Non_Storm_Data.gdb\\Vector_Grid", "vgLayer")               # make a feature layer of vector grid cells
        dm.SelectLayerByLocation("vgLayer", "INTERSECT", aoiBasin)                                      # select the vector grid cells that intersect the aoiBasin polygon
        dm.MakeFeatureLayer(home + "\\Input\\Non_Storm_Data.gdb\\Grid_Points", "gpLayer")               # make a feature layer of grid points
        dm.SelectLayerByLocation("gpLayer", "HAVE_THEIR_CENTER_IN", "vgLayer")                          # select the grid points within the vector grid selection
        con.FeatureClassToFeatureClass("gpLayer", env.scratchGDB, "PMP_Points")                         # save feature layer as "PMP_Points" feature class
        arcpy.AddMessage("(" + str(dm.GetCount("gpLayer")) + " grid points will be analyzed)\n")
        
        # Add PMP Fields
        for dur in durList:
            arcpy.AddMessage("\t...adding field: PMP_" + str(dur))
            dm.AddField(env.scratchGDB + "\\PMP_Points", "PMP_" + dur, "DOUBLE")

        # Add STORM Fields (this string values identifies the driving storm by SPAS ID number)
        for dur in durList:
            arcpy.AddMessage("\t...adding field: STORM_" + str(dur))
            dm.AddField(env.scratchGDB + "\\PMP_Points", "STORM_" + dur, "TEXT", "", "", 16, "Storm ID " + dur + "-hour")

        # Add STNAME Fields (this string values identifies the driving storm by SPAS ID number)
        for dur in durList:
            arcpy.AddMessage("\t...adding field: STNAME_" + str(dur))
            dm.AddField(env.scratchGDB + "\\PMP_Points", "STNAME_" + dur, "TEXT", "", "", 50, "Storm Name " + dur + "-hour")
        
        return

    def getAOIarea(aoiBasin, useBasinArea):
        if useBasinArea is True:
            dm.AddGeometryAttributes(aoiBasin, "AREA_GEODESIC", "", "SQUARE_MILES_US", 102039)
            totalArea = 0
            with arcpy.da.SearchCursor(aoiBasin, ['AREA_GEO']) as cursor:
                for row in cursor:
                    arcpy.AddMessage("\nAdd Polygon Geometry:  " + str(row[0]) + " sqmi")
                    totalArea += row[0]
            dm.DeleteField(aoiBasin, ['AREA_GEO'])
            if round(totalArea,1) < 1:
                aoiArea = round(totalArea, 1)
                arcpy.AddMessage("\n***Area used for PMP analysis: " + str(aoiArea) + " sqmi***")
                return (aoiArea)
            else:
                aoiArea = round(totalArea, 0)
                arcpy.AddMessage("\n***Area used for PMP analysis: " + str(aoiArea) + " sqmi***")
        if useBasinArea is False:
            aoiArea = arcpy.GetParameter(6)                                                     # Enable a constant area size
            aoiArea = round(aoiArea, 0)
            arcpy.AddMessage("\n***Area used for PMP analysis: " + str(aoiArea) + " sqmi***")

        return aoiArea
    ###########################################################################
    ##  Define dadLookup() function:
    ##  The dadLookup() function determines the DAD value for the current storm
    ##  and duration according to the basin area size.  The DAD depth is interpolated
    ##  linearly between the two nearest areal values within the DAD table.
    def dadLookup(stormLayer, duration, area):                  # dadLookup() accepts the current storm layer name (string), the current duration (string), and AOI area size (float)
        #arcpy.AddMessage("\t\tfunction dadLookup() called.")
        durField = "H_" + duration                              # defines the name of the duration field (eg., "H_06" for 6-hour)
        dadTable = dadGDB + "\\" + stormLayer
        rows = arcpy.SearchCursor(dadTable)
        
        try:       
            row = rows.next()                                       # Sets DAD area x1 to the value in the first row of the DAD table.
            x1 = row.AREASQMI
            y1 = row.getValue(durField)
            xFlag = "FALSE"                                         # xFlag will remain false for basins that are larger than the largest DAD area.
        except RuntimeError:                                        # return if duration does not exist in DAD table
            return
        
        row = rows.next()
        i = 0
        while row:                                                  # iterates through the DAD table - assiging the bounding values directly above and below the basin area size
            i += 1
            if row.AREASQMI < area:
                x1 = row.AREASQMI
                y1 = row.getValue(durField)
            else:
                xFlag = "TRUE"                                      # xFlag is switched to "TRUE" indicating area is within DAD range
                x2 = row.AREASQMI
                y2 = row.getValue(durField)
                break
            
            row = rows.next()
        del row, rows, i

        if xFlag == "FALSE":
            x2 = area                                           # If x2 is equal to the basin area, this means that the largest DAD area is smaller than the basin and the resulting DAD value must be extrapolated.            
            arcpy.AddMessage("\t\tThe basin area size: " + str(area) + " sqmi is greater than the largest DAD area: " + str(x1) + " sqmi.\n\t\tDAD value is estimated by extrapolation.")
            y = x1 / x2 * y1                                    # y (the DAD depth) is estimated by extrapolating the DAD area to the basin area size.
            return y                                            # The extrapolated DAD depth (in inches) is returned.

        # arcpy.AddMessage("\nArea = " + str(area) + "\nx1 = " + str(x1) + "\nx2 = " + str(x2) + "\ny1 = " + str(y1) + "\ny2 = " + str(y2))

        x = area                                                # If the basin area size is within the DAD table area range, the DAD depth is interpolated 
        deltax = x2 - x1                                        # to determine the DAD value (y) at area (x) based on next lower (x1) and next higher (x2) areas.
        deltay = y2 - y1
        diffx = x - x1

        y = y1 + diffx * deltay / deltax

        if x < x1:
            arcpy.AddMessage("\t\tThe basin area size: " + str(area) + " sqmi is less than the smallest DAD table area: " + str(x1) + " sqmi.\n\t\tDAD value is estimated by extrapolation.")
            
        return y                                                # The interpolated DAD depth (in mm) is returned.

    ###########################################################################
    ##  Define updatePMP() function:
    ##  This function updates the 'PMP_XX_' and 'STORM_XX' fields of the PMP_Points
    ##  feature class with the largest value from all analyzed storms stored in the
    ##  pmpValues list.
    def updatePMP(pmpValues, stormID, duration):                                                    # Accepts four arguments: pmpValues - largest adjusted rainfall for current duration (float list); stormID - driver storm ID for each PMP value (text list); and duration (string)  
        stormDict = {}                                                                              # Defines the storm list dictionary with 'SPAS_ID' as key and 'STORM_TXT' as value
        dictFields = ['SPAS_ID', 'STORM_TXT']
        stormList = home + "\\Input\\Non_Storm_Data.gdb\\Storm_List"
        with arcpy.da.SearchCursor(stormList, dictFields) as cursor:
            for row in cursor:
                stormDict[row[0]] = row[1]
        
        pmpfield = "PMP_" + duration
        stormfield = "STORM_" + duration
        stormTextField = "STNAME_" + duration    
        gridRows = env.scratchGDB + "\\PMP_Points"                                                  # iterates through PMP_Points rows
        i = 0
        with arcpy.da.UpdateCursor(gridRows, (pmpfield, stormfield, stormTextField)) as cursor:
            for row in cursor:
                row[0] = pmpValues[i]                                                               # Sets the PMP field value equal to the Max Adj. Rainfall value (if larger than existing value).
                row[1] = stormID[i]                                                                 # Sets the storm ID field to indicate the driving storm event
                if not stormID[i] == "STORM":                                                       # Conditional loop unless no controlling storm exists (ie, id of "STORM") to update table with storm text
                    if stormID[i].endswith("_LOC") or stormID[i].endswith("_GEN") or stormID[i].endswith("_TRO"):   # remove hybrid storm suffixes
                        stormNum = stormID[i][:-4]
                    else:
                        stormNum = stormID[i]                                                       
                    stormText = stormDict[stormNum]                                                 # Lookup 'STORM_TXT' in dictionary using stormNum (ie, 'SPAS_ID') as key
                    row[2] = stormText                                                              # Sets the storm text field ("STNAME_XX") to the storm text value
                cursor.updateRow(row)
                i += 1
        del row, gridRows, pmpfield, stormfield, stormTextField, i, dictFields, stormDict
        arcpy.AddMessage("\n\t" + duration + "-hour PMP values update complete. \n")
        return
    
    ###########################################################################        
    ##  The outputPMP() function produces raster GRID files for each of the PMP durations.
    ##  Aslo, a space-delimited PMP_Distribition.txt file is created in the 'Text_Output' folder.
    def outputPMP(type, area, outPath): 
        desc = arcpy.Describe(basin)
        basinName = desc.baseName
        pmpPoints = env.scratchGDB + "\\PMP_Points"                             # Location of 'PMP_Points' feature class which will provide data for output 

        outType = type[:1]
        outArea = str(int(round(area,0))) + "sqmi"
        outGDB = "PMP_"+ basinName + "_" + outArea +".gdb"                             
        # if not arcpy.Exists(outPath + "\\" + outGDB):                           # Check to see if PMP_XXXXX.gdb already exists
        #   arcpy.AddMessage("\nCreating output geodatabase '" + outGDB + "'")
        #   dm.CreateFileGDB(outPath, outGDB)
        arcpy.AddMessage("\nCreating output geodatabase '" + outGDB + "'")
        dm.CreateFileGDB(outPath, outGDB)
        arcpy.AddMessage("\nCopying PMP_Points feature class to " + outGDB + "...")
        con.FeatureClassToFeatureClass(pmpPoints, outPath + "\\" + outGDB, type + "_PMP_Points_" + basinName + "_" + outArea)
        pointFC = outPath + "\\" + outGDB + "\\" + type + "_PMP_Points_" + basinName + "_" + outArea    
        # addLayerMXD(pointFC)  # calls addLayerMDX function to add output to ArcMap session
        
        arcpy.AddMessage("\nBeginning PMP Raster Creation...")

        for dur in durList:                                                 # This code creates a raster GRID from the current PMP point layer
            durField = "PMP_" + dur
            outLoc = outPath + outGDB +"\\" + outType + "_" + dur + "_" + basinName + "_" + outArea
            arcpy.AddMessage("\n\tInput Path: " + pmpPoints)    
            arcpy.AddMessage("\tOutput raster path: " + outLoc)
            arcpy.AddMessage("\tField name: " + durField)
            con.FeatureToRaster(pmpPoints, durField, outLoc, gridRes)
            arcpy.AddMessage("\tOutput raster created...")               
        del durField, outLoc, dur

        arcpy.AddMessage("\nPMP Raster Creation complete.")

        if includeSubbasin:                                                 # Begin subbasin average calculations
            subbasinID = []
            with arcpy.da.SearchCursor(basin, subbasinIDfield) as cursor:   # Create list of subbasin ID names
                for row in cursor:
                    subbasinID.append(row[0])

            subIDtype = arcpy.ListFields(basin, subbasinIDfield)[0].type    # Define the datatype of the subbasin ID field

            if subIDtype != "String":                                       # Convert subbasin IDs to a string, if they are not already
                subbasinID = [str(i) for i in subbasinID]

            subNameLen = max(map(len, subbasinID))                          # Define the length of the longest subbasin ID
                
            # arcpy.AddMessage("\nList of subbasins...\n" + "\n".join(subbasinID))
            
            arcpy.AddMessage("\nCreating Subbasin Summary Table...")
            tableName = type + "_PMP_Subbasin_Average" + "_" + outArea
            tablePath = outPath + "\\" + outGDB + "\\" + tableName
            dm.CreateTable(outPath + "\\" + outGDB, tableName)              # Create blank table
      
            dm.AddField(tablePath, "STORM_TYPE", "TEXT", "", "", 10, "Storm Type")          # Create "Storm Type" field
            dm.AddField(tablePath, "SUBBASIN", "TEXT", "", "", subNameLen, "Subbasin")      # Create "Subbasin" field

            cursor = arcpy.da.InsertCursor(tablePath, "SUBBASIN")           # Create Insert cursor and add a blank row to the table for each subbasin
            for sub in subbasinID:
                cursor.insertRow([sub])
            del cursor, sub
            
            dm.CalculateField(tablePath, "STORM_TYPE", "'" + type + "'", "PYTHON_9.3")      # populate storm type field
                      
            i = 0
            for field in arcpy.ListFields(pmpPoints, "PMP_*"):              # Add fields for each PMP duration and calculate the subbasin averages
                fieldName = field.name
                arcpy.AddMessage("\n\tCalculating subbasin average for " + fieldName + " (weighted)...\n")                   
                dm.AddField(tablePath, fieldName, "DOUBLE", "", 2)          # Add duration field                
                subAveList = []
                for subbasin in subbasinID:                                 # Loop through each subbasin                  
                    if subIDtype != "String":                               # Define an SQL expression that specifies the current subbasin
                        sql_exp = """{0} = {1}""".format(arcpy.AddFieldDelimiters(basin, subbasinIDfield), subbasin)
                    else:
                        sql_exp = """{0} = '{1}'""".format(arcpy.AddFieldDelimiters(basin, subbasinIDfield), subbasin)   
                    dm.MakeFeatureLayer(basin, "subbasinLayer", sql_exp)
                    outLayer = outPath + "\\" + outGDB + "\\subbasin_" + str(subbasin)
                    subBasAve = basinAve("subbasinLayer", fieldName)        # Call the basAve() function passing the subbasin and duration field
                    arcpy.AddMessage("\tSubbasin average for " + str(subbasin) + ":  " + str(subBasAve) + '"') 
                    subAveList.append(subBasAve)                            # Add subbasin average to list
                p = 0
                with arcpy.da.UpdateCursor(tablePath, fieldName) as cursor: # Update the subbasin average summary table with the subbasin averages
                    for row in cursor:
                        row = subAveList[p]
                        cursor.updateRow([row])
                        p += 1
                i += 1
            arcpy.AddMessage("\nSubbasin summary table complete.")

        arcpy.AddMessage("\nCreating Basin Summary Table...")
        tableName = type + "_PMP_Basin_Average" + "_" + outArea
        tablePath = outPath + "\\" + outGDB + "\\" + tableName
        dm.CreateTable(outPath + "\\" + outGDB, tableName)          # Create blank table
        cursor = arcpy.da.InsertCursor(tablePath, "*")              # Create Insert cursor and add a blank row to the table
        cursor.insertRow([0])
        del cursor
        
        dm.AddField(tablePath, "STORM_TYPE", "TEXT", "", "", 30, "Storm Type")          # Create "Storm Type" field
        dm.CalculateField(tablePath, "STORM_TYPE", "'" + type + "'", "PYTHON_9.3")      # populate storm type field

        i = 0
        for field in arcpy.ListFields(pmpPoints, "PMP_*"):          # Add fields for each PMP duration and calculate the basin average
            fieldName = field.name
            fieldAve = basinAve(basin, fieldName)                   # Calls the basinAve() function - returns the average (weighted or not)
            dm.AddField(tablePath, fieldName, "DOUBLE", "", 2)      # Add duration field
            dm.CalculateField(tablePath, fieldName, fieldAve, "PYTHON_9.3")       # Assigns the basin average
            i += 1
        arcpy.AddMessage("\nSummary table complete.")
        basAveTables.append(tablePath)                              
        
        ##  The following lines export a .png image depth duration chart and PMP summary excel file to the output folder
        if ddChart:
            xValues = durList               #Get list of durations for chart
            xValues = [int(i) for i in xValues]     #Convert duration list to integers
            ax1 = plt.subplot2grid((1,1), (0,0))    #Create variable for subplot in chart
            yValues = []
            pmpFields = [field.name for field in arcpy.ListFields(tablePath, "PMP_*")] # Selects PMP fields for yValues
            with arcpy.da.SearchCursor(tablePath, pmpFields) as cursor:                # Adds PMP depths to yValues
                yValues = next(cursor)
            del cursor, pmpFields
            
            stormFields = [field.name for field in arcpy.ListFields(pmpPoints, "Storm_*")] # Selects Controlling Storm fields
            contStorms = []                             # List of controlling storms for a single duration
            listOfContStorms = []                       # List of controlling storms for all durations (list of lists)           
            i = 0                                       # iterator (for "Storm_*" fields)
            while i < len(stormFields):                      # iterates through controlling storm fields
                with arcpy.da.SearchCursor(pmpPoints, stormFields) as cursor:    # Search cursor returns list of unique controlling storms
                    contStorms = sorted({row[i] for row in cursor})
                listOfContStorms.append(contStorms)                         # Add unique storms for current duration to list of controlling stomrs             
                i += 1
            del cursor

            plt.plot(xValues,yValues)       #Creates chart
            plt.xlabel('Storm Duration in Hours')
            plt.ylabel('Rainfall Depth in Inches')
            plt.title(basinName + " (" + outArea + ") " + type + ' Storm Basin Average PMP\nDepth Duration Chart') 
            ax1.grid(True)              #Creates grid lines in chart
            yTop = max(yValues) + 1
            ax1.set_ylim(top = yTop)        #Sets y axis values to match depths +1 1
            ax1.set_xticks(xValues)     #Sets x axis values to match durations
            plt.rcParams["figure.figsize"] = (10,6)
            plt.tight_layout()
            plt.savefig(outPath + "\\" + basinName + "_" + outArea + "_"+ type + "_Depth_Duration_Chart.png")  #Save image
            plt.close()         #Close chart to remove from memory
            arcpy.AddMessage("\nDepth Duration Chart exported to output folder.")
            del xValues, yValues, #df, dfLimited
            return
        return

    ###########################################################################        
    ##  The basin() returns the basin average PMP value for a given duration field.
    ##  If the option for a weighted average is checked in the tool parameter the script
    ##  will weight the grid point values based on proportion of area inside the basin.
    def basinAve(aoiBasin, pmpField):
        pmpPoints = env.scratchGDB + "\\PMP_Points"                                                         # Path of 'PMP_Points' scratch feature class       
        if weightedAve:
            #arcpy.AddMessage("\tCalculating sub-basin average for " + pmpField + "(weighted)...")
            vectorGridClip = env.scratchGDB + "\\VectorGridClip"                                            # Path of 'VectorGridClip' scratch feature class
                
            dm.MakeFeatureLayer(home + "\\Input\\Non_Storm_Data.gdb\\Vector_Grid", "vgLayer")                # make a feature layer of vector grid cells
            dm.SelectLayerByLocation("vgLayer", "INTERSECT", aoiBasin)                                      # select the vector grid cells that intersect the aoiBasin polygon

            an.Clip("vgLayer", aoiBasin, vectorGridClip)                                    # clips aoi vector grid to basin
            dm.AddField(pmpPoints, "WEIGHT", "DOUBLE")                                                      # adds 'WEIGHT' field to PMP_Points scratch feature class
            dm.MakeFeatureLayer(vectorGridClip, "vgClipLayer")                                              # make a feature layer of basin clipped vector grid cells
            dm.MakeFeatureLayer(pmpPoints, "pmpPointsLayer")                                                # make a feature layer of PMP_Points feature class

            dm.AddJoin("pmpPointsLayer", "ID", "vgClipLayer", "ID")                                     # joins PMP_Points and vectorGridBasin tables
            dm.CalculateField("pmpPointsLayer", "WEIGHT", "!vectorGridClip.Shape_Area!", "PYTHON_9.3")      # Calculates basin area proportion to use as weight for each grid cell.
            dm.RemoveJoin("pmpPointsLayer", "vectorGridClip")
            dm.SelectLayerByLocation("pmpPointsLayer", "INTERSECT", "vgLayer")

            with arcpy.da.UpdateCursor("pmpPointsLayer",['WEIGHT']) as cursor:                          ## set 'weight' "nan" cells to 0
                for row in cursor:
                    if row[0] == None:
                        row[0] = 0
                        cursor.updateRow(row)
            del cursor            
                             
            na = arcpy.da.TableToNumPyArray("pmpPointsLayer",(pmpField, 'WEIGHT'))                          # Assign pmpPoints values and weights to Numpy array (na)
            wgtAve = np.ma.average(na[pmpField], weights=na['WEIGHT'])                                         # Calculate weighted average with Numpy average
            del na
            return round(wgtAve, 2)

        else:
            if includeSubbasin:
                #arcpy.AddMessage("\tCalculating sub-basin average for " + pmpField + "(non-weighted)...")
                vectorGridClip = env.scratchGDB + "\\VectorGridClip"                                            # Path of 'VectorGridClip' scratch feature class
                    
                dm.MakeFeatureLayer(home + "\\Input\\Non_Storm_Data.gdb\\Vector_Grid", "vgLayer")                # make a feature layer of vector grid cells
                dm.SelectLayerByLocation("vgLayer", "INTERSECT", aoiBasin)                                      # select the vector grid cells that intersect the aoiBasin polygon

                dm.MakeFeatureLayer(pmpPoints, "pmpPointsLayer")                                                # make a feature layer of PMP_Points feature class

                dm.SelectLayerByLocation("pmpPointsLayer", "INTERSECT", "vgLayer")
                                          
                na = arcpy.da.TableToNumPyArray("pmpPointsLayer", pmpField)                                     # Assign pmpPoints values and weights to Numpy array (na)
                fieldAve = np.average(na[pmpField])                                                             # Calculates aritmetic mean
                del na
                return round(fieldAve, 2)            
            
            else:
                arcpy.AddMessage("\tCalculating basin average for " + pmpField + "(not weighted)...")
                na = arcpy.da.TableToNumPyArray(pmpPoints, pmpField)                                            # Assign pmpPoints values to Numpy array (na)                     
                fieldAve = np.average(na[pmpField])                                                             # Calculates aritmetic mean
                del na
                return round(fieldAve, 2)


    ###########################################################################        
    ##  This basinZone() function returns a list containing transposition zone ID
    ##  (as an integer)


    def	basinZone(bas):	 ## This function returns the basin location transposition zone
        tempBasin = env.scratchGDB + "\\tempBasin"
        tempCentroid = env.scratchGDB + "\\tempCentroid"
        joinFeat = home + "\\Input\\Non_Storm_Data.gdb\\Vector_Grid"
        joinOutput = env.scratchGDB + "\\joinOut"
        dm.Dissolve(bas, tempBasin)
        desc = arcpy.Describe(tempBasin)
        sr = desc.spatialReference
        #dm.FeatureToPoint(tempBasin, tempCentroid, "INSIDE")

        dm.CreateFeatureclass(env.scratchGDB,"tempCentroid","POINT",spatial_reference = sr)
        with arcpy.da.InsertCursor(tempCentroid, "SHAPE@XY") as iCur:
            with arcpy.da.SearchCursor(tempBasin,"SHAPE@") as sCur:
                for sRow in sCur:
                    cent = sRow[0].centroid          # get the centroid
                    iCur.insertRow([(cent.X,cent.Y)])# write it to the new feature class
        
        an.SpatialJoin(tempCentroid, joinFeat, joinOutput)
        centZone = arcpy.da.SearchCursor(joinOutput, ("ZONE",)).next()[0]
        centDivide = arcpy.da.SearchCursor(joinOutput, ("DIVIDE",)).next()[0]
        del tempBasin, tempCentroid, joinFeat, joinOutput, desc, sr
        return (centZone, centDivide)


        ###########################################################################        
    ##  The temporalToAscii() function applies 
    ##
    def temporalToAscii(temporalTable, distributionFields, outPath, outArea, pmpDur, stormType):

        basinPMP = outPath + "\\" + stormType + "_PMP_Basin_Average_" + outArea                                # Location of basin average PMP table        pmpDur = 
        # pmpMax = arcpy.da.SearchCursor(basinPMP, ("PMP_120",)).next()[0]
        dm.CreateFolder(outLocation + "\\" + stormType, "ASCII_" + basinName + "_" + outArea)

        for distField in distributionFields:
            arcpy.AddMessage("\nConverting to ASCII GRIDs: " + distField)
            dm.CreateFolder(outLocation + "\\" + stormType + "\\ASCII_" + basinName + "_" + outArea, "PMP_" + distField)
            ascPath = outLocation + "\\" + stormType + "\\" + "ASCII_" + basinName + "_" + outArea + "\\" + "PMP_" + distField
            whereClause = distField + " IS NOT NULL"
            pmpMax = max([cur[0] for cur in arcpy.da.SearchCursor(temporalTable, distField, whereClause)])
            arcpy.AddMessage("\tMax PMP: " + str(pmpMax))
            accumFacts = []
            incFacts = []
            incFacts.append(arcpy.da.SearchCursor(temporalTable, distField).next()[0] / pmpMax)
            # arcpy.AddMessage("\nInitial Incremental: " + str(incFacts[0]))

            i = 0            
            with arcpy.da.SearchCursor(temporalTable, distField) as cursor:
                for row in cursor:
                    try:
                        accumFacts.append(row[0] / pmpMax)
                        # arcpy.AddMessage("\naccum: " + str(row[0] / pmpMax))
                        if i == 0:
                            # arcpy.AddMessage("Skiprow")
                            i += 1
                            continue
                        incFacts.append(accumFacts[i] - accumFacts[i - 1])
                        # arcpy.AddMessage("i: " + str(i) + "\n\tAccum: " + str(accumFacts[i]) + "\n\tInc: " + str(incFacts[i]))
                        i += 1
                    except:
                        continue
            del cursor, row, i
            # arcpy.AddMessage("\nAccum. Factors: " + str(accumFacts))
            # arcpy.AddMessage("\nInc. Factors: " + str(incFacts))

            i = 1

            with arcpy.da.SearchCursor(temporalTable, distField) as cursor:      ## Calculate number of timesteps based on distribution length
                n_timesteps = 0
                for row in cursor:
                    if row[0] is not None:
                        n_timesteps += 1
                        
            intDur = int(n_timesteps / (60 / timestep))    ## Calculate duration based on number timesteps divided by timestep length (in hours)
            pmpDur = str(intDur).zfill(2)                  ## Convert duration to text with leading zero
            pmpRaster = Raster(outPath + "\\" + stormType[:1] + "_" + pmpDur + "_" + basinName + "_" + outArea)
            for inc in incFacts:
                outRaster = outGDB + "\\PMP_" + distField + "_" + str(i).zfill(3)
                incRaster = pmpRaster * inc
                ascName = "PMP_" + basinName + "_" + distField + "_" + str(i).zfill(3)
                con.RasterToASCII(incRaster, ascPath + "\\" + ascName + ".asc")
                i += 1


    #     ###########################################################################        
    def createNetCDF(raster, start_date, n_timesteps, timeFreqMin, out_netcdf):
        # Get raster properties
        extent = raster.extent
        cell_size = round(raster.meanCellHeight,3)
        rows = raster.height
        cols = raster.width
        max_x = round(extent.XMax, 4) - cell_size / 2
        min_x = round(extent.XMin, 4) + cell_size / 2
        max_y = round(extent.YMax, 4) - cell_size / 2
        min_y = round(extent.YMin, 4) + cell_size / 2
        projection = raster.spatialReference

        # Create x and y coordinate arrays
        y_coords = np.flip(np.linspace(min_y, max_y, rows))    # flipped so array goes from top to bottom to match Raster format
        x_coords = np.linspace(min_x, max_x, cols)
        # arcpy.AddMessage("\ny_coords: " + str(y_coords) + "\nx_coords: " + str(x_coords)) 

        # creat xarray dataset with coordinates and time
        ds = xr.Dataset(coords={'latitude': y_coords,
                                'longitude': x_coords,
                                'time': pd.date_range(start=start_date, periods=n_timesteps, freq=timeFreqMin)})

        # Set time variable to standard calendar
        ds['time'].encoding['calendar'] = 'standard'

        # Add precipitation variable (initialize with zeros)
        ds['precipitation'] = xr.DataArray(np.zeros((n_timesteps, ds.latitude.size, ds.longitude.size)), dims=('time', 'latitude', 'longitude'), coords=ds.coords)
        # ds['precipitation'] = (('time', 'latitude', 'longitude'), np.full((n_timesteps, rows, cols), -9999, dtype=np.float32))
        ds['precipitation'].attrs['units'] = 'inches'
        ds['precipitation'].encoding['_FillValue'] = -9999

        # Add elevation variable (initialize with zeros)
        ds['elevation'] = xr.DataArray(np.zeros((ds.latitude.size, ds.longitude.size)), dims=('latitude', 'longitude'))
        ds['elevation'].attrs['units'] = 'ft'
        ds['elevation'].encoding['_FillValue'] = -9999

        # Set projection information
        ds.attrs['crs'] = projection.exportToString()

        # Save netCDF file
        ds.to_netcdf(out_netcdf)

        ###########################################################################        

    def temporalToNetCDF(temporalTable, distributionFields, outPath, outArea, pmpDur, stormType, timestep):
        dm.CreateFolder(outLocation + "\\" + stormType, "NetCDF_" + basinName + "_" + outArea)
        ncPath = outLocation + "\\" + stormType + "\\" + "NetCDF_" + basinName + "_" + outArea
        basinPMPPoints = outPath + "\\" + stormType + "_PMP_Points_" + basinName + "_" + outArea                                # Location of basin average PMP table

        ## Set start time based on seasonality
        if stormType == 'General' or stormType == 'Local' or stormType == 'Tropical':
            start_date = "2050-06-01 00:00:00"
#        elif stormType == 'CoolSeason':
#            start_date = "2050-03-01 00:00:00"

        ## Iterate through temporal distribiutions and create netCDF files for each
        for distField in distributionFields:
            ### Create output netCDF for temporal distributions  
            arcpy.AddMessage("\n\tCreating netCDF file for " + distField + " temporal distribution...")
            with arcpy.da.SearchCursor(temporalTable, distField) as cursor:      ## Calculate number of timesteps based on distribution length
                n_timesteps = 0
                for row in cursor:
                    if row[0] is not None:
                        n_timesteps += 1
            intDur = int(n_timesteps / (60 / timestep))    ## Calculate duration based on number timesteps divided by timestep length (in hours)
            pmpDur = str(intDur).zfill(2)                  ## Convert duration to text with leading zero
            pmpRaster = Raster(outPath + "\\" + stormType[:1] + "_" + pmpDur + "_" + basinName + "_" + outArea)         
            #arcpy.AddMessage("\tPMP Duration: " + str(pmpDur))
            out_netcdf = ncPath + "\\" + "PMP_" + distField + ".nc"   
            delta_time = datetime.timedelta(minutes=timestep)
            #n_timesteps = dm.GetCount(temporalTable)

            arcpy.AddMessage("\t\tCreating netCDF file: " + "PMP_" + distField + ".nc")
            createNetCDF(pmpRaster, start_date, n_timesteps, delta_time, out_netcdf)    # creates netCDF based on spatial template 

            ## Produce a list of incremental adjustment factors for current distribution            
            whereClause = distField + " IS NOT NULL"
            pmpMax = max([cur[0] for cur in arcpy.da.SearchCursor(temporalTable, distField, whereClause)])
            # arcpy.AddMessage("\tMax PMP: " + str(pmpMax))
            accumFacts = []
            incFacts = []
            incFacts.append(arcpy.da.SearchCursor(temporalTable, distField).next()[0] / pmpMax)
            i = 0            
            with arcpy.da.SearchCursor(temporalTable, distField) as cursor:
                for row in cursor:
                    try:
                        accumFacts.append(row[0] / pmpMax)
                        # arcpy.AddMessage("\naccum: " + str(row[0] / pmpMax))
                        if i == 0:
                            # arcpy.AddMessage("Skiprow")
                            i += 1
                            continue
                        incFacts.append(accumFacts[i] - accumFacts[i - 1])
                        # arcpy.AddMessage("i: " + str(i) + "\n\tAccum: " + str(accumFacts[i]) + "\n\tInc: " + str(incFacts[i]))
                        i += 1
                    except:
                        continue
            del cursor, row, i

            ## open netCDF xarray dataset for appending
            ds = xr.open_dataset(out_netcdf)       # open netCDF      
            
            arcpy.AddMessage("\t\tAdding precipitation values to netCDF...")
            i = 0
            for inc in incFacts:
                incRaster = pmpRaster * inc

                # Convert raster to NumPy array
                rasArray = arcpy.RasterToNumPyArray(incRaster, nodata_to_value=-9999)

                # Check if the raster dimensions match the netCDF dimensions
                if rasArray.shape[0] != ds.dims['latitude'] or rasArray.shape[1] != ds.dims['longitude']:
                    raise ValueError("Raster dimensions do not match netCDF dimensions.")
            
                # # Write data to the xarray precip variable at i timestep
                pptvar = ds.variables['precipitation']
                pptvar[i,:,:] = rasArray
                i += 1
                del rasArray

            ## Add elevation grid values
            arcpy.AddMessage("\t\tAdding elevation values to netCDF...")
            elevField = 'ELEV_FT'
            con.FeatureToRaster(basinPMPPoints, elevField, r'memory\elevgrid', gridRes)      ## Create temporary elevation raster as netCDF grid template
            elevArray = arcpy.RasterToNumPyArray(r'memory\elevgrid', nodata_to_value=-9999)
            dm.Delete(r'memory\elevgrid')
            elevVar = ds.variables['elevation']
            elevVar[:,:] = elevArray

            ## Write to netCDF file
            arcpy.AddMessage("\t\tSaving netcdf file: " + out_netcdf)
            ds.to_netcdf(out_netcdf)    # write xarray to netCDF
            dm.Delete("memory")

    ###########################################################################        
    ##  The temporalCritStacked() function applies the critically stacked
    ##  temporal distributions scenarios.  The function accepts the storm type,
    ##  output .gdb path, AOI area size, PMP duration string (hours), and
    ##  integer timestep duration (minutes). The function outputs a gdb table.       
    def temporalCritStacked(stormType, outPath, area, duration, timestep, tablePath, write):                                          # Function applied Critically Stacked temporal distribution
        basinPMP = outPath + "\\" + stormType + "_PMP_Basin_Average_" + area                                                        # Location of basin average PMP table
        # if stormType == "Local" and duration == "06":                                                                               # These conditional statements define the field name based on storm type, PMP duration, and timestep duration 
        #     csField = "LS_" + duration + "_HOUR_" + str(timestep) + "MIN_CRIT_STACKED"
        # elif stormType == "General":
        #     csField = "GS_" + duration + "_HOUR_" + str(timestep) + "MIN_CRIT_STACKED"
        # elif stormType == "Tropical":
        #     csField = "TS_" + duration + "_HOUR_" + str(timestep) + "MIN_CRIT_STACKED"
        # else:
        #     arcpy.AddMessage("\n***Stacked distribution not available for storm type: " + stormType)
        #     return
        csField = "Critically_Stacked_" + duration + "h"
        arcpy.AddMessage("\n\tApplying " + duration + "-hour " + str(timestep) + "-min " + csField + " Temporal Distribution...")
        # tableName = "Temporal_Distribution_" + duration + "hr_" + str(timestep) + "min_Crit_Stacked"                                # Output table name
        # tablePath = outPath + "\\" + tableName                                                                                      # Output table full path
        pmpFields = [field.name for field in arcpy.ListFields(basinPMP, "PMP*")]                                                    # Gets the "PMP_XX" field names from the basin avereage PMP table       
        if duration == "06":                                                                                                        # These conditional statements define the key durations needed to build the critically stacked patterns for the following durations...
            keyDurations = [1, 2, 3, 4, 5, 6]                                                                                       
        elif duration == "12":
            keyDurations = [1, 2, 3, 4, 5, 6, 12]            
        elif duration == "24":
            keyDurations = [1, 2, 3, 4, 5, 6, 12, 24]
        elif duration == "48":
            keyDurations = [1, 2, 3, 4, 5, 6, 12, 24, 48]
        elif duration == "72":
            keyDurations = [1, 2, 3, 4, 5, 6, 12, 24, 48, 72]
        elif duration == "96":
            keyDurations = [1, 2, 3, 4, 5, 6, 12, 24, 48, 72, 96]
        elif duration == "120":
            keyDurations = [1, 2, 3, 4, 5, 6, 12, 24, 48, 72, 96, 120]
        elif duration == "144":
            keyDurations = [1, 2, 3, 4, 5, 6, 12, 24, 48, 72, 96, 120, 144]
        elif duration == "168":
            keyDurations = [1, 2, 3, 4, 5, 6, 12, 24, 48, 72, 96, 120, 144, 168]     
        else:
            arcpy.AddMessage("\n\t\t***Critically stacked temporal distribution not available for " + duration + "-hour duration.***")
            return        
        timestepLen = int(duration) * 60 // timestep                                                                                 # number of rows in output table     
        xValues = [0]
        for i in keyDurations:                                                                                                      # defines the known x-values (xp) to be used in the interpolation
            xVal = i * timestepLen / int(duration)
            xValues.append(xVal)
        del i, xVal        
        yValues = [0]
        d = 0
        for i in keyDurations:                                                                                                      # defines the known y-values (fp) to be used in the interpolation
            pmpDepth = arcpy.da.SearchCursor(basinPMP, pmpFields).next()[d]
            yValues.append(pmpDepth)
            d += 1
        del d, i, pmpDepth

        x = np.arange(0, timestepLen + 1, 1)                                                                                        # defines the x points at which to interpolate values
        xp = np.asarray(xValues)                                                                                                    # np.asarray converts lists into numpy arrays
        fp = np.asarray(yValues)
        y = np.interp(x, xp, fp)

        # arcpy.AddMessage("\nyValues: " + str(yValues))
        # arcpy.AddMessage("\nxValues: " + str(xValues))
        # arcpy.AddMessage("\nxpValues: " + str(xp))
        # arcpy.AddMessage("\nfpValues: " + str(fp))
        # arcpy.AddMessage("\ni: " + str(y))
        
        inc = []
        prevDepth = 0
        i = 0
        for depth in np.nditer(y):                                                                                                  # populates incremental depths list 'inc' with y array
            inc.append(depth - prevDepth)
            prevDepth = depth
            i += 1
        del i, prevDepth
        periods = int(duration)                                                                                                     # defines number of periods (known hours) as the duration     
        periodLen = 60 // timestep                                                                                                  # defines number of timesteps (minutes) in each period
        ranks = []
        stackRank = 1
        i = 0
        while i < periods:                                                                                                          # populates list 'ranks' with a rank integer, one entry per period
            ranks.append(stackRank)
            stackRank += periodLen
            i += 1
        del i

        orderRanks = []
        orderRanks.insert(0, ranks.pop(0))
        for i in range(timestepLen // periodLen):                                                                                   ## orders the ranks according to critically stacked pattern.  Pulls
            if ranks:                                                                                                               ## (pop()) the first rank from the ranks list and places it in the orderRanks
                orderRanks.insert(0, ranks.pop(0))                                                                                  ## list. Places next two ranks at the beginning of the list
            if ranks:                                                                                                               ## and the following at the end of the list.  Repeats until ranks is empty.
                orderRanks.insert(0, ranks.pop(0))
            if ranks:
                orderRanks.append(ranks.pop(0))
        del i
        orderRanks += [orderRanks.pop(0)]
        if orderRanks[0] == max(orderRanks):                                                                                        ## Moves last rank to the end of of orderRanks list.
            arcpy.AddMessage("\n*** moving first rank to last")
            orderRanks.append(orderRanks.pop(max))      
        orderInc = []
        n = 0
        for i in range(periods):                                                                                                    # gets the nth largest increment where n is the ordered Rank.
            for q in range(periodLen):
                nthLargest = nlargest(orderRanks[n], inc)[-1]
                orderInc.append(nthLargest)
            n += 1
        del n, i, q
        # arcpy.AddMessage("\nPeriods: " + str(periods))
        # arcpy.AddMessage("\norderRanks: " + str(orderRanks))
        # arcpy.AddMessage("\nxorderInc: " + str(orderInc))
        if write == False:
            # arcpy.AddMessage("\n***Returning orderInc***")
            return csField, orderInc

        cumulative = []
        prevInc = 0
        step = 1
        # env.extent = arcpy.Extent(-127.8125, 52.9125, -123.8625, 54.5125)
        # arcpy.AddMessage("\n***OrderInc Sum:" + str(sum(orderInc)))
        # arcpy.AddMessage("\n***OrderInc:" + str(orderInc))
        for i in orderInc:                                                                                                       # Converts the incremental depths to cumulatove depths and places in cumulative list
            value = round(i + prevInc, 4)  
            cumulative.append(value)
            prevInc = i + prevInc
            incFactor = i / sum(orderInc)
            step += 1
            i += 1
        del i, prevInc        
        timesteps = x.tolist()                                                                                                      # Converts the timesteps array (x) to a list then removes the first zero entry
        timesteps.pop(0)
        minutes = []
        minutesInc = timestep
        for i in range(timestepLen):                                                                                                # Constructs the minutes list to be used in output column based on timestep interval
            minutes.append(minutesInc)
            minutesInc += timestep
        del i  
        # dm.CreateTable(outPath, tableName)                                                                                        # Create the output geodatabase table
        # dm.AddField(tablePath, "TIMESTEP", "DOUBLE")                                                                              # Create "TIMESTEP" field
        # dm.AddField(tablePath, "MINUTES", "DOUBLE")                                                                               # Create "MINUTES" field
        dm.AddField(tablePath, csField, "DOUBLE")                                                                                   # Create cumulated rainfall field     
        # zipped = zip(timesteps, minutes, cumulative)                                                                                # Zip up lists of output items.
        # arcpy.AddMessage("\n\tZipped: " + str(zipped)) 
        # fields = ('TIMESTEP', 'MINUTES', csField)                                                                                   # Output table field names        
        # arcpy.AddMessage("\n\tApplying temporal distribution for: " + csField)
        i = 0
        with arcpy.da.UpdateCursor(tablePath, csField) as cursor:                                                    # Cursor to populate output Critically Stacked table
            for row in cursor:
                try:
                    # arcpy.AddMessage("\n\tRow Value: " + str(cumulative[i]))
                    row[0] = cumulative[i]
                    cursor.updateRow(row)
                    i += 1
                except IndexError:                                                # Exits updatecursor once distribution factor is NULL
                    break        

        return csField, orderInc

    ###########################################################################        
    ##  The temporalDistHuff2() function applies the Huff Curve temporal distributions
    ##  scenarios.  The function accepts the storm type as the first argument, the temporal
    ##  distribution table (containing the temporal pattern factors) as the second
    ##  argument, and the output GDB, areaSize, and output table path as additional arguments.
    def temporalDistHuff02(stormType, temporalDistTable, outPath, areaSize, outTable):   
        basinPMP = outPath + "\\" + stormType + "_PMP_Basin_Average_" + areaSize                                # Location of basin average PMP table

        arcpy.AddMessage("\n\t***" + stormType + " Storm - PMP Temporal Distributions***")
        oneHour = arcpy.da.SearchCursor(basinPMP, ("PMP_01",)).next()[0]                    # Gets 1-hour PMP depth
        twoHour = arcpy.da.SearchCursor(basinPMP, ("PMP_02",)).next()[0]                    # Gets 2-hour PMP depth

        StormDistributionList = [field.name for field in arcpy.ListFields(temporalDistTable, "*2*")]
        arcpy.AddMessage("\n\t2-hour Huff Distribution Field Names: " + str(StormDistributionList))          

        for distribution in StormDistributionList:
            arcpy.AddMessage("\n\tApplying temporal distribution for: " + distribution)
            dm.AddField(outTable, distribution, "DOUBLE")
            distVals = [row[0] for row in arcpy.da.SearchCursor(temporalDistTable, distribution)]
            accumPMP = 0
            # i = 6
            # with arcpy.da.UpdateCursor(outTable, distribution) as cursor:             # Cursor to distribute 24h rainfall
            #     for row in cursor:
            #         try:
            #             accumPMP =  pmp * distVals[i]
            #             row[0] = round(accumPMP, 4)
            #             cursor.updateRow(row)
            #             i += 1
            #         except IndexError:                                                # Exits updatecursor once distribution factor is NULL
            #             break
            #     del row, cursor
            i = 6
            with arcpy.da.UpdateCursor(outTable, [distribution, "TIMESTEP"]) as cursor:                   # Cursor to apply temporal factor to 2-hour PMP
                for row in cursor:
                    if row[1] <= 6:                                             # Constrain update to rows 1-6                               
                        accumPMP += (twoHour - oneHour) / 12
                        row[0] = accumPMP
                        cursor.updateRow(row)
                    if row[1] > 6 and row[1] <= 18:                               # Constrain update to rows 7-18 
                        accumPMP += oneHour * distVals[i]
                        row[0] = accumPMP
                        cursor.updateRow(row)
                        i += 1
                    if row[1] > 18 and row[1] <= 24:                               # Constrain update to rows 19-24 
                        accumPMP +=  (twoHour - oneHour) / 12
                        row[0] = accumPMP
                        cursor.updateRow(row)
                del row, cursor
        del stormType, basinPMP
        return

    ###########################################################################        
    ##  The temporalDistHuff6() function applies the Huff Curve temporal distributions
    ##  scenarios.  The function accepts the storm type as the first argument, the temporal
    ##  distribution table (containing the temporal pattern factors) as the second
    ##  argument, and the output GDB, areaSize, and output table path as additional arguments.
    def temporalDistHuff06(stormType, temporalDistTable, outPath, areaSize, outTable):   
        basinPMP = outPath + "\\" + stormType + "_PMP_Basin_Average_" + areaSize                                # Location of basin average PMP table

        arcpy.AddMessage("\n\t***" + stormType + " Storm - PMP Temporal Distributions***")
        pmp = arcpy.da.SearchCursor(basinPMP, ("PMP_06",)).next()[0]

        StormDistributionList = [field.name for field in arcpy.ListFields(temporalDistTable, "*6*")]
        arcpy.AddMessage("\n\t6-hour Huff Distribution Field Names: " + str(StormDistributionList))          

        for distribution in StormDistributionList:
            arcpy.AddMessage("\n\tApplying temporal distribution for: " + distribution)
            dm.AddField(outTable, distribution, "DOUBLE")
            distVals = [row[0] for row in arcpy.da.SearchCursor(temporalDistTable, distribution)]
            accumPMP = 0
            i = 0
            with arcpy.da.UpdateCursor(outTable, distribution) as cursor:             # Cursor to distribute 24h rainfall
                for row in cursor:
                    if i < len(distVals):
                        accumPMP =  pmp * distVals[i]
                        row[0] = round(accumPMP, 4)
                        cursor.updateRow(row)
                        i += 1
                    else:
                        i += 1
                del row, cursor

        del stormType, basinPMP
        return

    ###########################################################################        
    ##  The temporalDistHuff12() function applies the Huff Curve temporal distributions
    ##  scenarios.  The function accepts the storm type as the first argument, the temporal
    ##  distribution table (containing the temporal pattern factors) as the second
    ##  argument, and the output GDB, areaSize, and output table path as additional arguments.
    def temporalDistHuff12(stormType, temporalDistTable, outPath, areaSize, outTable):   
        basinPMP = outPath + "\\" + stormType + "_PMP_Basin_Average_" + areaSize                                # Location of basin average PMP table

        arcpy.AddMessage("\n\t***" + stormType + " Storm - PMP Temporal Distributions***")

        StormDistributionList = [field.name for field in arcpy.ListFields(temporalDistTable, "*12*")]
        arcpy.AddMessage("\n\t12-hour Huff Distribution Field Names: " + str(StormDistributionList))          

        largest6 = arcpy.da.SearchCursor(basinPMP, ("PMP_06",)).next()[0]                          # Calculate largest 6-hour period PMP
        second6 = (arcpy.da.SearchCursor(basinPMP, ("PMP_12",)).next()[0] - largest6)/2            # Calculate the next largest 6-hr period PMP and divide by 2

        arcpy.AddMessage("\n\tLargest 6-hour Period: " + str(largest6))
        arcpy.AddMessage("\tFirst 3-hour: " + str(second6))
        arcpy.AddMessage("\tLast 3-hour: " + str(second6))

        for distribution in StormDistributionList:                                   # Loops thourgh each 12-hour temporal distribution
            arcpy.AddMessage("\n\tApplying temporal distribution for: " + distribution)
            dm.AddField(outTable, distribution, "DOUBLE")
            distVals = [row[0] for row in arcpy.da.SearchCursor(temporalDistTable, distribution)]
            arcpy.AddMessage("\t\tFirst 3-hour Period...")
            accumPMP = 0
            with arcpy.da.UpdateCursor(outTable, [distribution, "TIMESTEP"]) as cursor:             # Cursor to evenly distribute half of 2nd largest 6-hr into first 3 hours
                for row in cursor:
                    if row[1] <= 36:                                        # Leave loop once a row containing a temporal dist. factor (ie, first 3h period) is reached
                        accumPMP += second6 / 36
                        row[0] = accumPMP
                        cursor.updateRow(row)
                del row, cursor

            arcpy.AddMessage("\t\tLargest 6-hour Period...")
            i = 36
            with arcpy.da.UpdateCursor(outTable, [distribution, "TIMESTEP"]) as cursor:             # Cursor to apply temporal factors to largest 6-hour PMP
                for row in cursor:
                    if row[1] > 36 and row[1] <= 108:                               # Constrain update to timesteps 37-108 (middle 6hr period)
                        accumPMP = (largest6 * distVals[i]) + second6
                        row[0] = accumPMP
                        cursor.updateRow(row)
                        i += 1
                del row, cursor
                
            arcpy.AddMessage("\t\tLast 3-hour Period...")
            # whereClause = distribution + " IS NULL"
            with arcpy.da.UpdateCursor(outTable, [distribution, "TIMESTEP"]) as cursor:              # Cursor to evenly distribute half of 2nd largest 6-hr into last 3 hours
                for row in cursor:
                    if row[1] > 108 and row[1] <= 144:                               # Constrain update to timesteps 109-144 (middle 6hr period)
                        accumPMP += second6 / 36
                        row[0] = accumPMP
                        cursor.updateRow(row)
                del row, cursor, accumPMP
        return

    ###########################################################################        
    ##  The temporalDistHuff24() function applies the Huff Curve temporal distributions
    ##  scenarios.  The function accepts the storm type as the first argument, the temporal
    ##  distribution table (containing the temporal pattern factors) as the second
    ##  argument, and the output PMP GDB, areaSize, and output temporal distribution 
    ##  table path as additional arguments.
    def temporalDistHuff24(stormType, temporalDistTable, outPath, areaSize, outTable):                                   # General Storm Temporal Distributions Application Function
        basinPMP = outPath + "\\" + stormType + "_PMP_Basin_Average_" + areaSize                                # Location of basin average PMP table

        arcpy.AddMessage("\n\t***" + stormType + " Storm - PMP Temporal Distributions***")
        pmp = arcpy.da.SearchCursor(basinPMP, ("PMP_24",)).next()[0]

        StormDistributionList = [field.name for field in arcpy.ListFields(temporalDistTable, "*24*")]
        arcpy.AddMessage("\n\t24-hour Huff Distribution Field Names: " + str(StormDistributionList))          

        for distribution in StormDistributionList:
            arcpy.AddMessage("\n\tApplying temporal distribution for: " + distribution)
            dm.AddField(outTable, distribution, "DOUBLE")
            distVals = [row[0] for row in arcpy.da.SearchCursor(temporalDistTable, distribution)]
            accumPMP = 0
            i = 0
            with arcpy.da.UpdateCursor(outTable, distribution) as cursor:             # Cursor to distribute 24h rainfall
                for row in cursor:
                    # arcpy.AddMessage("\n\tDistVal: " + str(distVals[i])) 
                    try:                                                              
                        accumPMP =  pmp * distVals[i]
                        row[0] = round(accumPMP, 4)
                        cursor.updateRow(row)
                        i += 1
                    except IndexError:                                                # Exits updatecursor once distribution factor is NULL
                        break
                del row, cursor

        del stormType, basinPMP
        return

    ###########################################################################        
    ##  The temporalDistHuff24L() function applies the Local Storm Huff Curve temporal distributions
    ##  scenarios.  The function accepts the storm type as the first argument, the temporal
    ##  distribution table (containing the temporal pattern factors) as the second
    ##  argument, and the output PMP GDB, areaSize, and output temporal distribution 
    ##  table path as additional arguments.
    def temporalDistHuff24L(stormType, temporalDistTable, outPath, areaSize, outTable):                                   # General Storm Temporal Distributions Application Function
        basinPMP = outPath + "\\" + stormType + "_PMP_Basin_Average_" + areaSize                                # Location of basin average PMP table

        arcpy.AddMessage("\n\t***" + stormType + " Storm - PMP Temporal Distributions***")

        StormDistributionList = [field.name for field in arcpy.ListFields(temporalDistTable, "*24*")]
        arcpy.AddMessage("\n\t24-hour Huff Distribution Field Names: " + str(StormDistributionList))          

        largest6 = arcpy.da.SearchCursor(basinPMP, ("PMP_06",)).next()[0]                          # Calculate largest 6-hour period PMP
        second6 = (arcpy.da.SearchCursor(basinPMP, ("PMP_12",)).next()[0] - largest6)/2            # Calculate the next largest 6-hr period PMP and divide by 2
        last12 = (arcpy.da.SearchCursor(basinPMP, ("PMP_24",)).next()[0] - arcpy.da.SearchCursor(basinPMP, ("PMP_12",)).next()[0])/2           # Calculate the reamaining 12hr period PMP and divide by 2

        for distribution in StormDistributionList:                                   # Loops thourgh each 24-hour temporal distribution
            arcpy.AddMessage("\n\tApplying temporal distribution for: " + distribution)
            dm.AddField(outTable, distribution, "DOUBLE")
            distVals = [row[0] for row in arcpy.da.SearchCursor(temporalDistTable, distribution)]

            arcpy.AddMessage("\t\tFirst section...")
            accumPMP = 0
            
            with arcpy.da.UpdateCursor(outTable, [distribution, "TIMESTEP"]) as cursor:             # Cursor to evenly distribute half of last 12-hr into first 6 hours
                for row in cursor:
                    if row[1] <= 72:                                        # Leave loop once a row containing a temporal dist. factor (ie, first 6h period) is reached
                        accumPMP +=  last12 / 72
                        row[0] = accumPMP
                        cursor.updateRow(row)
                del row, cursor                
            
            arcpy.AddMessage("\t\tSecond section...")
            with arcpy.da.UpdateCursor(outTable, [distribution, "TIMESTEP"]) as cursor:             # Cursor to evenly distribute half of 2nd largest 6-hr into 3 hours
                for row in cursor:
                    if row[1] > 72 and row[1] <= 108:                                        # Leave loop once a row containing a temporal dist. factor is reached
                        accumPMP +=  second6 / 36
                        row[0] = accumPMP
                        cursor.updateRow(row)
                del row, cursor

            i = 108
            arcpy.AddMessage("\t\tLargest 6-hour Period...")
            with arcpy.da.UpdateCursor(outTable, [distribution, "TIMESTEP"]) as cursor:             # Cursor to apply temporal factors to largest 6-hour PMP
                for row in cursor:
                    if row[1] > 108 and row[1] <= 180:                               # Constrain update to rows 108-180 (middle 6hr period)
                        # arcpy.AddMessage("\t\tDistval: " + type(distVals[i]) + "\tlargest6: " + type(largest6))
                        accumPMP = (largest6 * distVals[i]) + second6 + last12
                        row[0] = accumPMP
                        cursor.updateRow(row)
                        i += 1
                del row, cursor

            arcpy.AddMessage("\t\tLast section...")
            with arcpy.da.UpdateCursor(outTable, [distribution, "TIMESTEP"]) as cursor:             # Cursor to evenly distribute half of 2nd largest 6-hr into 3 hours
                for row in cursor:
                    if row[1] > 180 and row[1] <= 216:                                        # Leave loop once a row containing a temporal dist. factor is reached
                        accumPMP +=  second6 / 36
                        row[0] = accumPMP
                        cursor.updateRow(row)
                del row, cursor

            arcpy.AddMessage("\t\tLast 6-hour Period...")
            whereClause = distribution + " IS NULL"
            with arcpy.da.UpdateCursor(outTable, distribution, whereClause) as cursor:              # Cursor to evenly distribute half of 2nd largest 6-hr into last 6 hours
                for row in cursor:
                    accumPMP +=  last12 / 72
                    row[0] = accumPMP
                    cursor.updateRow(row)
                del row, cursor, accumPMP, whereClause
        return


    ###########################################################################        
    ##  The temporalDistHuff48() function applies the Huff Curve temporal distributions
    ##  scenarios.  The function accepts the storm type as the first argument, the temporal
    ##  distribution table (containing the temporal pattern factors) as the second
    ##  argument, and the output GDB, areaSize, and output table path as additional arguments.
    def temporalDistHuff48(stormType, temporalDistTable, outPath, areaSize, outTable):                                   # General Storm Temporal Distributions Application Function
        basinPMP = outPath + "\\" + stormType + "_PMP_Basin_Average_" + areaSize                                # Location of basin average PMP table

        arcpy.AddMessage("\n\t***" + stormType + " Storm - PMP Temporal Distributions***")
        # pmp = arcpy.da.SearchCursor(basinPMP, ("PMP_72",)).next()[0]

        largest24 = arcpy.da.SearchCursor(basinPMP, ("PMP_24",)).next()[0]                                                          # Calculate largest 72-hour period PMP
        second24 = arcpy.da.SearchCursor(basinPMP, ("PMP_48",)).next()[0] - largest24                   # Calculate next largest 24-hour period PMP
        halfSecond24 = round(second24 / 2.0,2)        

        arcpy.AddMessage("\n\t\tLargest 24-hour: " + str(largest24))
        arcpy.AddMessage("\t\tSecond largest 24-hour: " + str(second24))
        arcpy.AddMessage("\t\tHalf of second largest 24-hour: " + str(halfSecond24))

        StormDistributionList = [field.name for field in arcpy.ListFields(temporalDistTable, "*48*")]
        arcpy.AddMessage("\n\t48-hr Huff Distribution Field Names: " + str(StormDistributionList))  
        for distribution in StormDistributionList:
            arcpy.AddMessage("\n\tApplying temporal distribution for: " + distribution)
            dm.AddField(outTable, distribution, "DOUBLE")
            distVals = [row[0] for row in arcpy.da.SearchCursor(temporalDistTable, distribution)]
            # arcpy.AddMessage("\nDistVals: " + str(distVals))
            arcpy.AddMessage("\t\t...first 12-hour Period...")
            accumPMP = 0
            with arcpy.da.UpdateCursor(outTable, [distribution, "TIMESTEP"]) as cursor:             # Cursor to evenly distribute first 12hr block
                for row in cursor:
                    if row[1] <= 48:                                        # First 12-hour block
                        accumPMP +=  halfSecond24 / 48
                        # arcpy.AddMessage("\tAccumulated Rain: " + str(round(accumPMP, 3)))
                        row[0] = round(accumPMP, 4)
                        cursor.updateRow(row)
                del row, cursor
            arcpy.AddMessage("\t\t...middle 24-hour Period...")
            i = 48
            with arcpy.da.UpdateCursor(outTable, [distribution, "TIMESTEP"]) as cursor:             # Cursor to apply temporal factors to largest 24-hour PMP
                for row in cursor:
                    if row[1] > 48 and row[1] <= 144:                               # Middle 24-hour block
                        accumPMP = round((largest24 * distVals[i]) + halfSecond24, 3)
                        # arcpy.AddMessage("\tAccumulated Rain: " + str(accumPMP))
                        row[0] = round(accumPMP, 4)
                        cursor.updateRow(row)
                        i += 1
                del row, cursor                    
            arcpy.AddMessage("\t\t...final 12-hour Period...")
            # whereClause = distribution + " IS NULL"
            # with arcpy.da.UpdateCursor(outTable, distribution, whereClause) as cursor:              # Cursor to evenly distribute final 12hr block
            with arcpy.da.UpdateCursor(outTable, [distribution, "TIMESTEP"]) as cursor:              # Cursor to evenly distribute final 12hr block
                for row in cursor:
                    if row[1] > 144 and row[1] <= 192:                               # Final 12-hour block                       
                        accumPMP +=  halfSecond24 / 48
                        # arcpy.AddMessage("\tAccumulated Rain: " + str(round(accumPMP, 3)))
                        row[0] = round(accumPMP, 4)
                        cursor.updateRow(row)
                del row, cursor, accumPMP

        del stormType, basinPMP
        return

    ###########################################################################        
    ##  The temporalDistHuff72() function applies the Huff Curve temporal distributions
    ##  scenarios.  The function accepts the storm type as the first argument, the temporal
    ##  distribution table (containing the temporal pattern factors) as the second
    ##  argument, and the output GDB, areaSize, and output table path as additional arguments.
    def temporalDistHuff72(stormType, temporalDistTable, outPath, areaSize, outTable):                                   # General Storm Temporal Distributions Application Function
        basinPMP = outPath + "\\" + stormType + "_PMP_Basin_Average_" + areaSize                                # Location of basin average PMP table

        arcpy.AddMessage("\n\t***" + stormType + " Storm - PMP Temporal Distributions***")
        # pmp = arcpy.da.SearchCursor(basinPMP, ("PMP_72",)).next()[0]

        largest24 = arcpy.da.SearchCursor(basinPMP, ("PMP_24",)).next()[0]                                                          # Calculate largest 72-hour period PMP
        second24 = arcpy.da.SearchCursor(basinPMP, ("PMP_48",)).next()[0] - largest24                   # Calculate next largest 24-hour period PMP
        third24 = arcpy.da.SearchCursor(basinPMP, ("PMP_72",)).next()[0] - arcpy.da.SearchCursor(basinPMP, ("PMP_48",)).next()[0]   # Calculate 3rd-largest 24-hour period PMP

        arcpy.AddMessage("\n\t\tLargest 24-hour: " + str(largest24))
        arcpy.AddMessage("\t\tSecond largest 24-hour: " + str(second24))
        arcpy.AddMessage("\t\tThird largest 24-hour: " + str(third24))
        # arcpy.AddMessage("\t\t24-hour PMP: " + str(pmp24))
        # arcpy.AddMessage("\t\t48-hour PMP: " + str(pmp48))      
        # arcpy.AddMessage("\t\t120-hour PMP: " + str(pmp120))
        StormDistributionList = [field.name for field in arcpy.ListFields(temporalDistTable, "*72*")]
        arcpy.AddMessage("\n\t72-hr Huff Distribution Field Names: " + str(StormDistributionList))  
        for distribution in StormDistributionList:
            arcpy.AddMessage("\n\tApplying temporal distribution for: " + distribution)
            dm.AddField(outTable, distribution, "DOUBLE")
            distVals = [row[0] for row in arcpy.da.SearchCursor(temporalDistTable, distribution)]
            # arcpy.AddMessage("\nDistVals: " + str(distVals))
            arcpy.AddMessage("\t\t...first 24-hour Period...")
            accumPMP = 0
            with arcpy.da.UpdateCursor(outTable, [distribution, "TIMESTEP"]) as cursor:             # Cursor to evenly distribute 2nd largest hour
                for row in cursor:
                    if row[1] <= 96:                                        # Leave loop once a row containing a temporal dist. factor (ie, second 24h period) is reached
                        accumPMP +=  second24 / 96
                        # arcpy.AddMessage("\tAccumulated Rain: " + str(round(accumPMP, 3)))
                        row[0] = round(accumPMP, 4)
                        cursor.updateRow(row)
                del row, cursor
            arcpy.AddMessage("\t\t...middle 24-hour Period...")
            i = 96
            with arcpy.da.UpdateCursor(outTable, [distribution, "TIMESTEP"]) as cursor:             # Cursor to apply temporal factors to largest 24-hour PMP
                for row in cursor:
                    if row[1] > 96 and row[1] <= 192:                               # Constrain update to rows 96-192 (second 24hr period)
                        accumPMP = round((largest24 * distVals[i]) + second24, 3)
                        # arcpy.AddMessage("\tAccumulated Rain: " + str(accumPMP))
                        row[0] = round(accumPMP, 4)
                        cursor.updateRow(row)
                        i += 1
                del row, cursor                    
            arcpy.AddMessage("\t\t...final 24-hour Period...")
            whereClause = distribution + " IS NULL"
            with arcpy.da.UpdateCursor(outTable, distribution, whereClause) as cursor:              # Cursor to evenly distribute 3nd largest hour over remaining empty rows
                for row in cursor:
                    accumPMP +=  third24 / 96
                    # arcpy.AddMessage("\tAccumulated Rain: " + str(round(accumPMP, 3)))
                    row[0] = round(accumPMP, 4)
                    cursor.updateRow(row)
                del row, cursor, accumPMP, whereClause

        del stormType, basinPMP
        return
    
    def temporalDistControlStorm_06hr(stormType, temporalDistTable, outPath, areaSize, outTable):                         # Local Storm 6-hr Temporal Distributions Function
        basinPMP = outPath + "\\" + stormType + "_PMP_Basin_Average_" + areaSize
        basinPMPPoints = outPath + "\\" + stormType + "_PMP_Points_" + basinName + "_" + areaSize                                # Location of basin average PMP table
        pmp06 = arcpy.da.SearchCursor(basinPMP, ("PMP_06",)).next()[0]

        StormDistributionList = []
        with arcpy.da.SearchCursor(basinPMPPoints, "STORM_06") as cursor:       # Get list of unique controlling storms for specified PMP field
            StormDistributionList = sorted({row[0] for row in cursor})

        arcpy.AddMessage("\n\t6-hour Controlling Storm Distribution Field Names: " + str(StormDistributionList))          

        for distribution in StormDistributionList:                                   # Loops thourgh each 24-hour temporal distribution
            arcpy.AddMessage("\n\tApplying temporal distribution for: " + distribution)
            if distribution in [f.name for f in arcpy.ListFields(temporalDistTable)]:
                dm.AddField(outTable, distribution + "_6hr", "DOUBLE")
                distVals = [row[0] for row in arcpy.da.SearchCursor(temporalDistTable, distribution)]
                # arcpy.AddMessage("\nDistVals: " + str(distVals))
                accumPMP = 0
                i = 0
                with arcpy.da.UpdateCursor(outTable, distribution + "_6hr") as cursor:             # Cursor to distribute 120h rainfall
                    for row in cursor:
                        try:                                                # Exits updatecursor once distribution factor is NULL
                            accumPMP =  pmp06 * distVals[i]
                            row[0] = round(accumPMP, 4)
                            cursor.updateRow(row)
                            i += 1
                        except IndexError:                                                # Exits updatecursor once distribution factor is NULL
                            i += 1
                    del row, cursor
            else:
                arcpy.AddMessage("\n\t\t" + distribution + " temporal distribution not available.  Skipping...")     
        return

    def temporalDistControlStorm_24hr(stormType, temporalDistTable, outPath, areaSize, outTable):                         # Local Storm 6-hr Temporal Distributions Function
        basinPMP = outPath + "\\" + stormType + "_PMP_Basin_Average_" + areaSize
        basinPMPPoints = outPath + "\\" + stormType + "_PMP_Points_" + basinName + "_" + areaSize                                # Location of basin average PMP table
        pmp24 = arcpy.da.SearchCursor(basinPMP, ("PMP_24",)).next()[0]

        StormDistributionList = []
        with arcpy.da.SearchCursor(basinPMPPoints, "STORM_24") as cursor:       # Get list of unique controlling storms for specified PMP field
            StormDistributionList = sorted({row[0] for row in cursor})

        arcpy.AddMessage("\n\t24-hour Controlling Storm Distribution Field Names: " + str(StormDistributionList))          

        for distribution in StormDistributionList:                                   # Loops thourgh each 24-hour temporal distribution
            arcpy.AddMessage("\n\tApplying temporal distribution for: " + distribution)
            if distribution in [f.name for f in arcpy.ListFields(temporalDistTable)]:
                dm.AddField(outTable, distribution + "_24hr", "DOUBLE")
                distVals = [row[0] for row in arcpy.da.SearchCursor(temporalDistTable, distribution)]
                # arcpy.AddMessage("\nDistVals: " + str(distVals))
                accumPMP = 0
                i = 0
                with arcpy.da.UpdateCursor(outTable, distribution + "_24hr") as cursor:             # Cursor to distribute 120h rainfall
                    for row in cursor:
                        try:                                                # Exits updatecursor once distribution factor is NULL
                            accumPMP =  pmp24 * distVals[i]
                            row[0] = round(accumPMP, 4)
                            cursor.updateRow(row)
                            i += 1
                        except IndexError:                                                # Exits updatecursor once distribution factor is NULL
                            i += 1
                    del row, cursor
            else:
                arcpy.AddMessage("\n\t\t" + distribution + " temporal distribution not available.  Skipping...")    
        return

    def temporalDistControlStorm_72hr(stormType, temporalDistTable, outPath, areaSize, outTable):                         # Local Storm 6-hr Temporal Distributions Function
        basinPMP = outPath + "\\" + stormType + "_PMP_Basin_Average_" + areaSize
        basinPMPPoints = outPath + "\\" + stormType + "_PMP_Points_" + basinName + "_" + areaSize                                # Location of basin average PMP table
        pmp72 = arcpy.da.SearchCursor(basinPMP, ("PMP_72",)).next()[0]

        StormDistributionList = []
        with arcpy.da.SearchCursor(basinPMPPoints, "STORM_72") as cursor:       # Get list of unique controlling storms for specified PMP field
            StormDistributionList = sorted({row[0] for row in cursor})

        arcpy.AddMessage("\n\t72-hour Controlling Storm Distribution Field Names: " + str(StormDistributionList))          

        for distribution in StormDistributionList:                                   # Loops thourgh each 72-hour temporal distribution
            arcpy.AddMessage("\n\tApplying temporal distribution for: " + distribution)
            if distribution in [f.name for f in arcpy.ListFields(temporalDistTable)]:
                dm.AddField(outTable, distribution + "_72hr", "DOUBLE")
                distVals = [row[0] for row in arcpy.da.SearchCursor(temporalDistTable, distribution)]
                # arcpy.AddMessage("\nDistVals: " + str(distVals))
                accumPMP = 0
                i = 0
                with arcpy.da.UpdateCursor(outTable, distribution + "_72hr") as cursor:             # Cursor to distribute 120h rainfall
                    for row in cursor:
                        try:
                            accumPMP =  pmp72 * distVals[i]
                            row[0] = round(accumPMP, 4)
                            cursor.updateRow(row)
                            i += 1
                        except IndexError:                                                # Exits updatecursor once distribution factor is NULL
                            i += 1
                    del row, cursor
            else:
                arcpy.AddMessage("\n\t\t" + distribution + " temporal distribution not available.  Skipping...")    
        return

    def checkTemporal(stormType, outPath, temporalTableName, distributionFields, dur, areaSize):
        TemporalTable = outPath + "\\" + temporalTableName
        basinPMP = outPath + "\\" + stormType + "_PMP_Basin_Average_" + areaSize                                # Location of basin average PMP table
        pmpFields = [field.name for field in arcpy.ListFields(basinPMP, "PMP_*")]                               # PMP duration run       
        # temporalFields = [field.name for field in arcpy.ListFields(TemporalTable)]
        # table = arcpy.Describe(TemporalTable)
        # tableName = table.name

        pmp = []                                                                                                #Creates empty list and updates with PMP values for each duration run
        i = 0
        while i < len(pmpFields):
            with arcpy.da.SearchCursor(basinPMP,pmpFields) as cursor:
                for row in cursor:
                    pmp.append(row[i])
                    i += 1
        del i, cursor

        # checkTable = outPath + "\\Temporal_Distribution_Check_" + stormType
        checkTable = TemporalTable + "_Check"
        # arcpy.AddMessage("\nCheckTable: "  + checkTable)
        maxFields = []                                                                                          #Create Max fields for each duration
        checkFields = []                                                                                        #Create Check fields for each duration
        if arcpy.Exists(checkTable):
            with arcpy.da.InsertCursor(checkTable, "PATTERN") as cursor:
                for val in distributionFields:
                    cursor.insertRow([val])
            i = 0                                                                                                   #Populate fields
            for pmpField in pmpFields:
                with arcpy.da.UpdateCursor(checkTable, pmpField) as cursor:
                    for row in cursor:
                        row = pmp[i]
                        cursor.updateRow([row])
                    i += 1
            del i, cursor           
        else:
            checkTable = dm.CreateTable(outPath, temporalTableName + "_Check")                       #Creates table in output GDB, adds field, and populates field with distributions
            dm.AddField(checkTable, "PATTERN", "TEXT", "", "", 50)
            with arcpy.da.InsertCursor(checkTable, "PATTERN") as cursor:
                for val in distributionFields:
                    cursor.insertRow([val])
            for maxField in pmpFields:
                newField = maxField.replace("PMP","MAX")
                maxFields.append(newField)
            del newField
            for checkField in pmpFields:
                newField = checkField.replace("PMP","CHECK")
                checkFields.append(newField)
            del newField
            i = 0                                                                                                   #Populate fields
            for pmpField in pmpFields:
                dm.AddField(checkTable, pmpField, "DOUBLE", "", "", 50)
                dm.AddField(checkTable, maxFields[i], "DOUBLE", "", "", 50)
                dm.AddField(checkTable, checkFields[i], "TEXT", "", "", 50)
                with arcpy.da.UpdateCursor(checkTable, pmpField) as cursor:
                    for row in cursor:
                        row = pmp[i]
                        cursor.updateRow([row])
                    i += 1
            del i, cursor

        step = arcpy.da.SearchCursor(TemporalTable,("MINUTES",)).next()[0]   
        if step == 15:
            dic = {"01": 4, "02": 8, "03": 12, "04": 16, "05": 20, "06": 24, "12": 48, "24": 96, "48": 192, "72": 288, "96": 384, "120": 480}  # Dictionary to convert durations into 15-minute timesteps
        elif step == 5:
            dic = {"01": 12, "02": 24, "03": 36, "04": 48, "05": 60, "06": 72, "12": 144, "24": 288, "48": 576, "72": 864, "96": 1152, "120": 1440}
        elif step == 60:
            dic = {"01": 1, "02": 2, "03": 3, "04": 4, "05": 5, "06": 6, "12": 12, "24": 24, "48": 48, "72": 72, "96": 96, "120": 120, "144": 144, "168": 168}
        # arcpy.AddMessage(str(step) + " Minute distribution Pattern.....")

        maxFields = [field.name for field in arcpy.ListFields(checkTable, "MAX*")]
        i = 0                                        # Calculates incremental PMP depths from temporal distribution and gets maximum rainfall for each duration run
        d = durList.index(dur) + 1
        for dur in durList[:d]:
            k = dic[dur]
            p = 0                                    # Skip first 3 fields in temporaltable (objectID, Timesteps, minutes)
            for distribution in distributionFields:
                # arcpy.AddMessage("\nCheck Dist Field: " + distribution)
                incPMP = []
                previousRow = 0
                # with arcpy.da.SearchCursor(TemporalTable, distributionFields) as cursor:
                with arcpy.da.SearchCursor(TemporalTable, distributionFields) as cursor:
                    for row in cursor:
                        try:
                            increment = row[p] - previousRow
                            previousRow = row[p]
                            incPMP.append(increment)
                        except:
                            continue
                # arcpy.AddMessage("\nincPMP: " + str(incPMP))
                na = np.array(incPMP)
                sumList = np.convolve(na,np.ones(k))
                maxPMP = max(sumList)
                maximumPMP = math.trunc(maxPMP * 10 ** 2.0) / 10 ** 2.0
                p += 1
                with arcpy.da.UpdateCursor(checkTable, ["PATTERN", maxFields[i]]) as cursor:                # Updates table with max values
                    for row in cursor:
                        if row[0] == distribution:
                            row[1] = maximumPMP
                            cursor.updateRow(row)
            i += 1
        del i, k, cursor
        with arcpy.da.UpdateCursor(checkTable, '*') as cursor:                              # Compares PMP values to max values for each duration.  If PMP values are larger update check field with PASS if not FAIL
            for row in cursor:
                rec = dict(zip(cursor.fields, row))
                arcpy.AddMessage("\n\n\tChecking temporally distributed depth-durations against PMP: " + rec['PATTERN'] + "\n")
                for k, v in rec.items():
                    if not k.startswith('PMP_'):
                        continue
                    _, n = k.split('_')
                    try:                                                            # This try/except skips comparisons for additional durations not present in current temporal pattern
                        mx = rec['MAX_{}'.format(n)]
                        rec['CHECK_{}'.format(n)] = 'EXCEED' if v/mx-1 < -0.05 else 'OK'
                    except:
                        arcpy.AddMessage("\n\tDuration not present...")
                        continue
                    if rec['CHECK_{}'.format(n)] == 'EXCEED':
                        arcpy.AddMessage("\t" + str(n) + "-hour \n\t\tPMP value is... " + str(v) + "  \n\t\tmax rainfall value is..." + str(mx) + "\n\t\tThis distribution.... " + rec['CHECK_{}'.format(n)]+ "\n\t\t***Max " + str(n) + "-hour values exceed PMP values for this temporal distribution***")
                    else:
                        arcpy.AddMessage("\t" + str(n) + "-hour \n\t\tPMP value is... " + str(v) + "  \n\t\tmax rainfall value is..." + str(mx) + "\n\t\tThis distribution.... " + rec['CHECK_{}'.format(n)])
                cursor.updateRow([rec[k] for k in cursor.fields])
        del cursor, k, v, rec
        return    

    ###########################################################################
    ##  This portion of the code iterates through each storm feature class in the
    ##  'Storm_Adj_Factors' geodatabase (evaluating the feature class only within
    ##  the Local, General, or Tropical feature dataset).  For each duration,
    ##  at each grid point within the aoi basin, the transpositionality is
    ##  confirmed.  Then the DAD precip depth is retrieved and applied to the
    ##  total adjustement factor to yield the total adjusted rainfall.  This
    ##  value is then sent to the updatePMP() function to update the 'PMP_Points'
    ##  feature class.
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~##   


    desc = arcpy.Describe(basin)                                                        # Check to ensure AOI input shape is a Polygon. If not - exit. 
    basinShape = desc.shapeType
    if desc.shapeType == "Polygon":
        arcpy.AddMessage("\nBasin shape type: " + desc.shapeType)
    else:
        arcpy.AddMessage("\nBasin shape type: " + desc.shapeType)
        arcpy.AddMessage("\nError: Input shapefile must be a polygon!\n")
        sys.exit()
    
    createPMPfc()                                                                       # Call the createPMPfc() function to create the PMP_Points feature class.

    env.workspace = adjFactGDB                                                          # the workspace environment is set to the 'Storm_Adj_Factors' file geodatabase
    
    # dm.MakeFeatureLayer(home + "\\Input\\Non_Storm_Data.gdb\\Vector_Grid", "vgLayer")               # make a feature layer of vector grid cells
    # dm.SelectLayerByLocation("vgLayer", "INTERSECT", basin)                                      # select the vector grid cells that intersect the aoiBasin polygon

    aoiSQMI = round(getAOIarea(basin, useBasinArea),2)                                                     # Calls the getAOIarea() function to assign area of AOI shapefile to 'aoiSQKM'
    # if aoiSQMI > 100 and stormType is "Local":
    #     arcpy.AddMessage("\n***Warning - Local storm PMP depths only valid for basins 100 square miles or smaller***")
        
    stormList = arcpy.ListFeatureClasses("", "Point", stormType)                        # List all the total adjustment factor feature classes within the storm type feature dataset.
    for dur in durList:
        arcpy.AddMessage("\n*************************************************************\nEvaluating " + dur + "-hour duration...")

        pmpList = []
        driverList = []
        gridRows = arcpy.SearchCursor(env.scratchGDB + "\\PMP_Points")
        try:
            for row in gridRows:
                pmpList.append(0.0)                                                         # creates pmpList of empty float values for each grid point to store final PMP values
                driverList.append("STORM")                                                  # creates driverList of empty text values for each grid point to store final Driver Storm IDs
            del row, gridRows
        except UnboundLocalError:
            arcpy.AddMessage("\n***Error: No data present within basin/AOI area.***\n")
            sys.exit()

        env.workspace = adjFactGDB
        for storm in stormList[:]:
            arcpy.AddMessage("\n\tEvaluating storm: " + storm + "...") 
            dm.MakeFeatureLayer(storm, "stormLayer")                                    # creates a feature layer for the current storm
            dm.SelectLayerByLocation("stormLayer", "HAVE_THEIR_CENTER_IN", "vgLayer")   # examines only the grid points that lie within the AOI
            pmpField = "PMP_" + dur
            i = 0
            try:
                dadPrecip = round(dadLookup(storm, dur, aoiSQMI),3)
                arcpy.AddMessage("\t\t" + dur + "-hour DAD value:  " + str(dadPrecip) + "in")
            except TypeError:                                                           # In no duration exists in the DAD table - move to the next storm
                arcpy.AddMessage("\t***Duration '" + str(dur) + "-hour' is not present for " + str(storm) + ".***\n")
                continue    
            arcpy.AddMessage("\t\tComparing " + storm + " adjusted rainfall values against current driver values...")
            transCounter = 0                                                    # Counter for number of grid points transposed to

            for row in arcpy.da.SearchCursor("stormLayer", ['ID','TAF','TRANS']):
                if row[2] == 1:                                              # Only continue if grid point is transpositionable ('1' is transpostionable, '0' is not).
                    try:                                                        # get total adj. factor if duration exists
                        # arcpy.AddMessage("\nRow ID: " + str(row[0]))
                        transCounter += 1
                        adjRain = round(dadPrecip * row[1],2)
                        if adjRain > pmpList[i]:
                            pmpList[i] = adjRain
                            driverList[i] = storm
                    except RuntimeError:
                        arcpy.AddMessage("\t\t   *Warning*  Total Adjusted Raifnall value falied to set for row " + str(row[0]))
                        break
                    del adjRain 
                i += 1
            if transCounter == 0:
                arcpy.AddMessage("\t\tStorm not transposable to basin. Removing " + storm + " from list...\n")
                stormList.remove(storm)
            else:
                arcpy.AddMessage("\t\tTransposed to " + str(transCounter) + "/" + str(i) + " grid points...\n")
            del row, transCounter
        del storm, dadPrecip
        updatePMP(pmpList, driverList, dur)              # calls function to update "PMP Points" feature class  
    del dur, pmpList, stormList
    
    arcpy.AddMessage("\n'PMP_Points' Feature Class 'PMP_XX' fields update complete for all '" + stormType + "' storms.")
    
    outputPMP(stormType, aoiSQMI, outputPath)               # calls outputPMP() function

    basinName = desc.baseName
    outArea = str(int(round(aoiSQMI,0))) + "sqmi"
    outGDB = outLocation + "\\" + stormType + "\\PMP_" + basinName + "_" + outArea + ".gdb"
    
    if includeTemporal:
        centroidLocation = basinZone(basin)
        arcpy.AddMessage("\nBasin Centroid Transposition Zone: " + str(centroidLocation[0]) + "\nBasin Centroid side of Appalachian Divide: " + str(centroidLocation[1]))
#        if arcpy.CheckExtension("Spatial") == "Available":
        arcpy.CheckOutExtension("Spatial")

        ### Apply the 5/15-min Huff and Crit Stacked Distributions
        arcpy.AddMessage("\n\n***Applying Huff Temporal Distributions***")
        ## Define durations and timesteps to be used for temporal distributions
        if stormType == 'General' or stormType == 'Tropical':
            tempDurs = ["24", "48", "72"]
            timestep = 15 
            tableName = stormType + "_Temporal_Distributions_" + str(timestep) + "min"                   # Output table name
            tablePath = outGDB + "\\" + tableName                               # Output table full path
            dm.CreateTable(outGDB, tableName)                                   # Create the output geodatabase table
            dm.AddField(tablePath, "TIMESTEP", "DOUBLE")                        # Create "TIMESTEP" field
            dm.AddField(tablePath, "MINUTES", "DOUBLE")                         # Create "MINUTES" field
        elif stormType == 'Local':
            tempDurs = ["02", "06", "12", "24"]
            timestep = 5
            tableName = stormType + "_Temporal_Distributions_" + str(timestep) + "min"                   # Output table name
            tablePath = outGDB + "\\" + tableName                               # Output table full path
            dm.CreateTable(outGDB, tableName)                                   # Create the output geodatabase table
            dm.AddField(tablePath, "TIMESTEP", "DOUBLE")                        # Create "TIMESTEP" field
            dm.AddField(tablePath, "MINUTES", "DOUBLE")                         # Create "MINUTES" field

        ## Add TIMESTEP and MINUTES fields to output temporal distribution table and add 'timestepLen' number of rows
        timestepLen = int(tempDurs[-1]) * 60 // timestep       # get the length of rows in temporal dist tables                                                                          # number of rows in output table     
        xValues = [0]
        x = np.arange(0, timestepLen + 1, 1)                                                                                        # Create cumulated rainfall field     
        timesteps = x.tolist()                                                                                                      # Converts the timesteps array (x) to a list then removes the first zero entry
        timesteps.pop(0)
        minutes = []
        minutesInc = timestep
        for i in range(timestepLen):                                                                                                # Constructs the minutes list to be used in output column based on timestep interval
            minutes.append(minutesInc)
            minutesInc += timestep
        del i  
        zipped = zip(timesteps, minutes)
        fields = ('TIMESTEP', 'MINUTES')
        with arcpy.da.InsertCursor(tablePath, fields) as cursor:                                                                    # Cursor to populate output Critically Stacked table
            for i in zipped:
                cursor.insertRow(i)
        del cursor, i

        ## iterate through applicable temporal durations and apply apply appropriate temporal patterns:
        for tempDur in tempDurs:    
            if stormType == 'General':
                if tempDur == "24":
                    temporalDistTable = home + "\\Input\\Non_Storm_Data.gdb\\GS_TEMPORAL_DISTRIBUTIONS_24HR_" + str(centroidLocation[1])
                    temporalDistHuff24(stormType, temporalDistTable, outGDB, outArea, tablePath)
                    temporalCritStacked(stormType, outGDB, outArea, tempDur, timestep, tablePath, True)                #Calls Crit Stacked temporal function with a x-min timestep
                if tempDur == "48":
                    temporalDistTable = home + "\\Input\\Non_Storm_Data.gdb\\GS_TEMPORAL_DISTRIBUTIONS_48HR_" + str(centroidLocation[1])
                    temporalDistHuff48(stormType, temporalDistTable, outGDB, outArea, tablePath)
                    temporalCritStacked(stormType, outGDB, outArea, tempDur, timestep, tablePath, True)                #Calls Crit Stacked temporal function with a x-min timestep
                if tempDur == "72":
                    temporalDistTable = home + "\\Input\\Non_Storm_Data.gdb\\GS_TEMPORAL_DISTRIBUTIONS_72HR_" + str(centroidLocation[1])
                    temporalDistHuff72(stormType, temporalDistTable, outGDB, outArea, tablePath)
                    temporalCritStacked(stormType, outGDB, outArea, tempDur, timestep, tablePath, True)                #Calls Crit Stacked temporal function with a x-min timestep
            if stormType == 'Tropical':
                if tempDur == "24":
                    temporalDistTable = home + "\\Input\\Non_Storm_Data.gdb\\TS_TEMPORAL_DISTRIBUTIONS_24HR_" + str(centroidLocation[1])
                    temporalDistHuff24(stormType, temporalDistTable, outGDB, outArea, tablePath)
                    temporalCritStacked(stormType, outGDB, outArea, tempDur, timestep, tablePath, True)                #Calls Crit Stacked temporal function with a x-min timestep
                if tempDur == "48":
                    temporalDistTable = home + "\\Input\\Non_Storm_Data.gdb\\TS_TEMPORAL_DISTRIBUTIONS_48HR_" + str(centroidLocation[1])
                    temporalDistHuff48(stormType, temporalDistTable, outGDB, outArea, tablePath)
                    temporalCritStacked(stormType, outGDB, outArea, tempDur, timestep, tablePath, True)                #Calls Crit Stacked temporal function with a x-min timestep
                if tempDur == "72":
                    temporalDistTable = home + "\\Input\\Non_Storm_Data.gdb\\TS_TEMPORAL_DISTRIBUTIONS_72HR_" + str(centroidLocation[1])
                    temporalDistHuff72(stormType, temporalDistTable, outGDB, outArea, tablePath)
                    temporalCritStacked(stormType, outGDB, outArea, tempDur, timestep, tablePath, True)                #Calls Crit Stacked temporal function with a x-min timestep               
            elif stormType == 'Local':
                if tempDur == "02":
                    temporalDistTable = home + "\\Input\\Non_Storm_Data.gdb\\LS_TEMPORAL_DISTRIBUTIONS_02HR"
                    temporalDistHuff02(stormType, temporalDistTable, outGDB, outArea, tablePath)                    
                if tempDur == "06":
                    temporalDistTable = home + "\\Input\\Non_Storm_Data.gdb\\LS_TEMPORAL_DISTRIBUTIONS_06HR_" + str(centroidLocation[1])
                    temporalDistHuff06(stormType, temporalDistTable, outGDB, outArea, tablePath)
                    temporalCritStacked(stormType, outGDB, outArea, tempDur, timestep, tablePath, True)  
                if tempDur == "12":
                    temporalDistTable = home + "\\Input\\Non_Storm_Data.gdb\\LS_TEMPORAL_DISTRIBUTIONS_12HR_" + str(centroidLocation[1])
                    temporalDistHuff12(stormType, temporalDistTable, outGDB, outArea, tablePath)
                    temporalCritStacked(stormType, outGDB, outArea, tempDur, timestep, tablePath, True)
                if tempDur == "24":
                    temporalDistTable = home + "\\Input\\Non_Storm_Data.gdb\\LS_TEMPORAL_DISTRIBUTIONS_24HR_" + str(centroidLocation[1])
                    temporalDistHuff24L(stormType, temporalDistTable, outGDB, outArea, tablePath)
                    temporalCritStacked(stormType, outGDB, outArea, tempDur, timestep, tablePath, True)

        distributionFields = [field.name for field in arcpy.ListFields(tablePath)][3:]
        checkTemporal(stormType, outGDB, tableName, distributionFields, tempDur, outArea)       # Run intermediate duration temporal check function
        if exportASCII:
            arcpy.AddMessage("\nProducing ASCII-format temporal distribution output for " + tableName)                
            temporalToAscii(tablePath, distributionFields, outGDB, outArea, tempDur, stormType)    # Export ASCII gridded output
        if exportNetCDF:
            arcpy.AddMessage("\nProducing netCDF temporal distribution output for " + tableName)
            temporalToNetCDF(tablePath, distributionFields, outGDB, outArea, tempDur, stormType, timestep)    # Export netCDF gridded output

        ### Apply the 60-min Controlling Storm Distributions
        arcpy.AddMessage("\n\n***Applying Controlling Storm Temporal Distributions***")
        ## Define durations and timesteps to be used for temporal distributions
        if stormType == 'General' or stormType == 'Tropical':
            tempDurs = ["24", "72"]
            timestep = 60 
            tableName = stormType + "_Temporal_Distributions_" + str(timestep) + "min_Controlling_Storms"                   # Output table name
            tablePath = outGDB + "\\" + tableName                               # Output table full path
            dm.CreateTable(outGDB, tableName)                                   # Create the output geodatabase table
            dm.AddField(tablePath, "TIMESTEP", "DOUBLE")                        # Create "TIMESTEP" field
            dm.AddField(tablePath, "MINUTES", "DOUBLE")                         # Create "MINUTES" field

        elif stormType == 'Local':
            tempDurs = ["06", "24"]
            timestep = 60
            tableName = stormType + "_Temporal_Distributions_" + str(timestep) + "min_Controlling_Storms"                   # Output table name
            tablePath = outGDB + "\\" + tableName                               # Output table full path
            dm.CreateTable(outGDB, tableName)                                   # Create the output geodatabase table
            dm.AddField(tablePath, "TIMESTEP", "DOUBLE")                        # Create "TIMESTEP" field
            dm.AddField(tablePath, "MINUTES", "DOUBLE")                         # Create "MINUTES" field

        ## Add TIMESTEP and MINUTES fields to output temporal distribution table and add 'timestepLen' number of rows
        timestepLen = int(tempDurs[-1]) * 60 // timestep       # get the length of rows in temporal dist tables                                                                          # number of rows in output table     
        xValues = [0]
        x = np.arange(0, timestepLen + 1, 1)                                                                                        # Create cumulated rainfall field     
        timesteps = x.tolist()                                                                                                      # Converts the timesteps array (x) to a list then removes the first zero entry
        timesteps.pop(0)
        minutes = []
        minutesInc = timestep
        for i in range(timestepLen):                                                                                                # Constructs the minutes list to be used in output column based on timestep interval
            minutes.append(minutesInc)
            minutesInc += timestep
        del i  
        zipped = zip(timesteps, minutes)
        fields = ('TIMESTEP', 'MINUTES')
        with arcpy.da.InsertCursor(tablePath, fields) as cursor:                                                                    # Cursor to populate output Critically Stacked table
            for i in zipped:
                cursor.insertRow(i)
        del cursor, i

        ## iterate through applicable temporal durations and apply apply appropriate temporal patterns:
        for tempDur in tempDurs:    
            if stormType == 'General' or stormType == 'Tropical':
                if tempDur == "72":
                    temporalDistTable = home + "\\Input\\Non_Storm_Data.gdb\\CONTROLLING_STORM_TEMPORAL_DISTRIBUTIONS_72"
                    temporalDistControlStorm_72hr(stormType, temporalDistTable, outGDB, outArea, tablePath)
                if tempDur == "24":
                    temporalDistTable = home + "\\Input\\Non_Storm_Data.gdb\\CONTROLLING_STORM_TEMPORAL_DISTRIBUTIONS_24"
                    temporalDistControlStorm_24hr(stormType, temporalDistTable, outGDB, outArea, tablePath)
            elif stormType == 'Local':                 
                if tempDur == "06":
                    temporalDistTable = home + "\\Input\\Non_Storm_Data.gdb\\CONTROLLING_STORM_TEMPORAL_DISTRIBUTIONS_06"
                    temporalDistControlStorm_06hr(stormType, temporalDistTable, outGDB, outArea, tablePath)
                if tempDur == "24":
                    temporalDistTable = home + "\\Input\\Non_Storm_Data.gdb\\CONTROLLING_STORM_TEMPORAL_DISTRIBUTIONS_24"
                    temporalDistControlStorm_24hr(stormType, temporalDistTable, outGDB, outArea, tablePath)

        distributionFields = [field.name for field in arcpy.ListFields(tablePath)][3:]
        if distributionFields:
            checkTemporal(stormType, outGDB, tableName, distributionFields, tempDur, outArea)       # Run intermediate duration temporal check function
        if exportASCII:
            arcpy.AddMessage("\nProducing ASCII-format temporal distribution output for " + tableName)
            temporalToAscii(tablePath, distributionFields, outGDB, outArea, tempDur, stormType)    # Export ASCII gridded output
        if exportNetCDF:
            arcpy.AddMessage("\nProducing netCDF temporal distribution output for " + tableName)
            temporalToNetCDF(tablePath, distributionFields, outGDB, outArea, tempDur, stormType, timestep)    # Export netCDF gridded output

#        else:
#            arcpy.AddError("***Spatial Analyst extension required for temporal distribtion.  Exiting...***")
#            arcpy.AddMessage(arcpy.GetMessages(0))
#            raise LicenseError
#            # sys.exit(0)
#        arcpy.CheckInExtension("Spatial")

    i = 0                                                                                   #Creates CSV files of all output tables
    csvPath = outLocation + "\\" + stormType + "\\CSV_" + desc.baseName + "_" + outArea + "\\" 
    if not arcpy.Exists(outLocation + "\\" + stormType + "\\CSV_" + desc.baseName + "_" + outArea):
        arcpy.CreateFolder_management(outLocation + "\\" + stormType + "\\", "CSV_" + desc.baseName + "_" + outArea)
    arcpy.AddMessage("\n\t...Creating output tables as CSV files.. ")
    env.workspace = outGDB
    outTables = arcpy.ListTables()
    arcpy.AddMessage("...Tables: " + str(outTables))
    for t in outTables:
        arcpy.TableToTable_conversion(t, csvPath, outTables[i] + ".csv")
        i += 1
    csvFiles = os.listdir(csvPath)
    for file in csvFiles:
        if file.endswith(".xml"):
            os.remove(os.path.join(csvPath,file))
    del csvFiles, outTables, aoiSQMI #,stormType

    return

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~##     


if locDurations:
    type = "Local"
    durations = locDurations
    dm.CreateFolder(outLocation, type)
    outputPath = outLocation + "\\Local\\"  
    arcpy.AddMessage("\nRunning PMP analysis for storm type: " + type)
    pmpAnalysis(basin, type, durations)          # Calls the pmpAnalysis() function to calculate the local storm PMP
    arcpy.AddMessage("\nLocal storm analysis complete...\n*********************************************************************************************************")

if genDurations:
    type = "General"
    durations = genDurations
    dm.CreateFolder(outLocation, type)
    outputPath = outLocation + "\\General\\"  
    arcpy.AddMessage("\nRunning PMP analysis for storm type: " + type)
    pmpAnalysis(basin, type, durations)          # Calls the pmpAnalysis() function to calculate the general storm PMP
    arcpy.AddMessage("\nGeneral storm analysis complete...\n*********************************************************************************************************")

if tropDurations:
    type = "Tropical"
    durations = tropDurations
    dm.CreateFolder(outLocation, type)
    outputPath = outLocation + "\\Tropical\\"  
    arcpy.AddMessage("\nRunning PMP analysis for storm type: " + type)
    pmpAnalysis(basin, type, durations)          # Calls the pmpAnalysis() function to calculate the Tropical storm PMP
    arcpy.AddMessage("\nTropical storm analysis complete...\n*********************************************************************************************************")

