#! /usr/bin/env python
#
# ===========================================================
#
# Alexander R Davies
#
# Ocean Exploration, Remote Sensing, and Biogeography Lab
# School of Marine Science and Policy
# College of Earth, Ocean, and Environment
# University of Delaware
# ardavies@udel.edu
#
# Latest Update: 01/13/2014
#
# NOTES:    1) Does all the OSCAR Data Processing
#
#               
#
# ===========================================================
# Import Modules
# ===========================================================
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from datetime import datetime
from mpl_toolkits.basemap import Basemap, addcyclic
from mpl_toolkits.basemap import shiftgrid, cm
from scipy.ndimage.filters import minimum_filter, maximum_filter
from netCDF4 import Dataset
import glob
import os
from scipy.io import netcdf
import matplotlib.transforms as mtransforms
from matplotlib.patches import FancyBboxPatch
# from matplotlib.transforms import Bboxcx
from matplotlib.path import Path
import matplotlib.patches as patches
import math as ma
import csv
from matplotlib.ticker import MaxNLocator
from matplotlib import rc, rcParams
from numpy.random import uniform, seed
from matplotlib.mlab import griddata
from numpy import arange,array,ones#,random,linalg
from pylab import plot,show
from scipy import interpolate
from matplotlib.colors import LogNorm
from matplotlib.backends.backend_pdf import PdfPages
#
# ===========================================================
# Function to input delX, delY, return bearing and magnitude
# ===========================================================
#
def find_bearing_magnitude(X, Y):
    import math as m
    if ((Y >= 0) and (X >=0 )):
        bearing_rads = m.atan(abs(X)/abs(Y))
        bearing = m.degrees(bearing_rads)
        magnitude = abs(X)/m.sin(bearing_rads)

    elif ((Y < 0) and (X >=0)):
        bearing_rads = m.atan(abs(Y)/abs(X))
        bearing = m.degrees(bearing_rads) + 90
        magnitude = abs(Y)/m.sin(bearing_rads)

    elif ((Y < 0) and (X < 0)):
        bearing_rads = m.atan(abs(X)/abs(Y))
        bearing = m.degrees(bearing_rads) + 180
        magnitude = abs(X)/m.sin(bearing_rads)

    elif ((Y >= 0) and (X < 0)):
        bearing_rads = m.atan(abs(Y)/abs(X))
        bearing = m.degrees(bearing_rads) + 270
        magnitude = abs(Y)/m.sin(bearing_rads)
    else:
        print "#### ERROR IN: FIND BEARING AND MAG FUNCTION ####"

    return bearing_rads, bearing, magnitude
#
# ===========================================================
# Function Determine new Lat/Lon from Old + Bearing + dist
# ===========================================================
#
def destpoint(lat1, lon1, brng, d):
    import math

    R = 6374 #Radius of the Earth
    # brng = Bearing converted to radians.
    # d = Distance in km
    # lat1 = Current lat point converted to radians
    # lon1 = Current lon point converted to radians

    lat2 = math.asin(math.sin(lat1)*math.cos(d/R) + math.cos(lat1)*math.sin(d/R)*math.cos(brng))
    lon2 = lon1 + math.atan2(math.sin(brng)*math.sin(d/R)*math.cos(lat1), math.cos(d/R)-math.sin(lat1)*math.sin(lat2))
    
    lat22 = math.degrees(lat2)
    lon22 = math.degrees(lon2)

    return lat22, lon22
#
# ===========================================================
# Function to read csv files
# ===========================================================
#
def csvread(filename):
    with open(filename, 'rb') as f:
        reader = csv.reader(f, delimiter=',')
        nrow = 0;
        for row in reader:
            nrow = nrow +1
            ncol = len(row)
    array = np.zeros([ncol,nrow]) 
    with open(filename, 'rb') as f:
        reader = csv.reader(f, delimiter=',')
        counter = 0
        for row in reader:
            for jj in range(0,ncol):
                array[jj,counter] = row[jj]
            counter = counter + 1
    return array, nrow, ncol   
#
# ===========================================================
# Function to Find the nearest value
# ===========================================================
#
def find_nearest(array,value):
    idx = (np.abs(array-value)).argmin()
    return array[idx], idx
#
# ===========================================================
# Function to Calculate Distances on a Sphere in kilometers
# ===========================================================
#
def pos2dist(pos):
    # Pos = [lon1,lat1,lon2,lat2]
    import math as m
    R_aver = 6374
    pi = 3.14
    deg2rad = pi/180

    lag3 = pos[1] * deg2rad
    lon3 = pos[0] * deg2rad
    lag4 = pos[3] * deg2rad
    lon4 = pos[2] * deg2rad
    dist = R_aver * m.acos(m.cos(lag3)*m.cos(lag4)*m.cos(lon3-lon4) + m.sin(lag3)*m.sin(lag4))
    # Returns a distance in KM
    return dist
#
#
# ===================================================================================
# Function for Bilinear Interpolation
# ===================================================================================
#
def bilinear_interpolation(x, y, points):
    #
    # Order Points vy x, then by y
    points = sorted(points)               # order points by x, then by y
    (x1, y1, q11), (_x1, y2, q12), (x2, _y1, q21), (_x2, _y2, q22) = points
    #
    # Are the points a rectangle?
    if x1 != _x1 or x2 != _x2 or y1 != _y1 or y2 != _y2:
        raise ValueError('points do not form a rectangle')
    if not x1 <= x <= x2 or not y1 <= y <= y2:
        raise ValueError('(x, y) not within the rectangle')
    #
    # Bilinear Interp and return value.
    return (q11 * (x2 - x) * (y2 - y) +
            q21 * (x - x1) * (y2 - y) +
            q12 * (x2 - x) * (y - y1) +
            q22 * (x - x1) * (y - y1)
           ) / ((x2 - x1) * (y2 - y1) + 0.0)
#
def findmaxs(data, depths, arraylength, rowlengths):
    #
    # Array of Maxes
    maxes = np.zeros((3,int(arraylength)))
    #
    #
    for r1 in range(0,int(arraylength)):        
        if r1 == 0:     
            start = 0
            dumby = np.zeros(int(rowlengths[r1]))
            r22 = start
            for r2 in range(0, int(rowlengths[r1])):
                r22 = start + r2
                dumby[r2] = abs(data[r22])
            mx = np.amax(dumby)

            r33 = start
            for r3 in range(0, int(rowlengths[r1])):
                r33 = start + r3
                if (abs(data[r33]) == mx):
                    break
            maxes[1,1] = data[r33]
            maxes[0,1] = depths[r33]
            maxes[2,1] = r33
            start = r22 + 1
            
        else:
            dumby = np.zeros(int(rowlengths[r1]))
            r22 = start
            for r2 in range(0, int(rowlengths[r1])):
                r22 = start + r2
                dumby[r2] = abs(data[r22])
            mx = np.amax(dumby)
            
            r33 = start
            for r3 in range(0, int(rowlengths[r1])):
                r33 = start + r3
                if (abs(data[r33]) == mx):
                    break
            maxes[0,r1] = depths[r33]
            maxes[1,r1] = data[r33]
            maxes[2,r1] = r33
            start = r22 + 1
    return maxes
#
# ===================================================================================
# Initial Set-up to read float data
# ===================================================================================
#
# Get to the float data...
os.chdir('/data/orbprocess_mail/alex/Oct02_Dec25_2013/data/type/')
#
# Compile a list of the filenames and sort them based on the time they were generated
fullnames = glob.glob('Full*')
gradnames = glob.glob('Grad*')
infonames = glob.glob('Text*')
fullnames.sort(key = lambda x: os.path.getmtime(x))
gradnames.sort(key = lambda x: os.path.getmtime(x))
infonames.sort(key = lambda x: os.path.getmtime(x))
#
# Array Initialization
arraylen =  len(infonames)
floatdate = np.zeros(arraylen)
floatlat = [None]*arraylen
floatlon = [None]*arraylen
day = np.zeros(arraylen)
month=[None]*arraylen
year=[None]*arraylen
daystr=[None]*arraylen
listdate =[None]*arraylen
#
# ===========================================================
# Processing Float Time and Location Data 
# ===========================================================
#
# Retrive the Information for Each Data Profile
for d in range(0,int(arraylen)):
    #
    # Get the Date Info
    f = open('/data/orbprocess_mail/alex/Oct02_Dec25_2013/data/type/' + infonames[d], 'r')
    while 1:
        line = f.readline()
        if line[:4] == "Date":
            month[d] = str(line[5:7])
            day[d] = int(line[8:10])
            daystr[d] = str(line[8:10])
            year[d] = str(line[11:15])
            break
    #
    # Get the LaT and Lon
    f = open('/data/orbprocess_mail/alex/Oct02_Dec25_2013/data/type/' + infonames[d], 'r')
    while 1:
        line = f.readline()
        if line[:3] == "Lat":
            floatlat[d] = str(line[5:12])
            break
    f = open('/data/orbprocess_mail/alex/Oct02_Dec25_2013/data/type/' + infonames[d], 'r')
    while 1:
        line = f.readline()
        if line[:3] == "Lon":
            floatlon[d] = str(line[5:12])
            break
    #
    # Get the DOY for each profile
    if (year[d] == '2013'):
        yr = 0
    elif (year[d] == '2014'):
        yr = 365
    elif (year[d] == '2015'):
        yr = 365*2
    if(month[d] == '01'):
        mon = 0
    elif(month[d] == '02'):
        mon = 31
    elif(month[d] == '03'):
        mon = 31 + 28
    elif(month[d] == '04'):
        mon = (31 + 28 + 31)
    elif(month[d] == '05'):
        mon = (31 + 28 + 31 + 30)
    elif(month[d] == '06'):
        mon = (31 + 28 + 31 + 30 + 31)
    elif(month[d] == '07'):
        mon = (31 + 28 + 31 + 30 + 31 + 30)
    elif(month[d] == '08'):
        mon = (31 + 28 + 31 + 30 + 31 + 30 + 31)
    elif(month[d] == '09'):
        mon = (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31)
    elif(month[d] == '10'):
        mon = (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30)
    elif(month[d] == '11'):
        mon = (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31)
    elif(month[d] == '12'):
        mon = (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31 + 30)
    else:
        print 'invalid month'
    floatdate[d] = yr + mon + day[d]
# Get Back to start directory
os.chdir('/home/ardavies/satdata/OSCAR')
#
# ===================================================================================
# Reading the U and V Currents from NETCDF file
# ===================================================================================
#
# Change the longitude format and make lat/lon floats
for i in range(0,arraylen):
    floatlon[i] = 360 + float(floatlon[i])
    floatlat[i] = float(floatlat[i])
#
# ncdump file header to .txt file
os.system('ncdump -h oscar-third-52ec05e6dea57.nc > oscar-third-52ec05e6dea57.txt')
#
# open the netcdf file to read
ncf = netcdf.netcdf_file('oscar-third-52ec05e6dea57.nc','r')
#
# Getting the netcdf data into the corresponding arrays
datadates = ncf.variables['time'][:].squeeze()
v = ncf.variables['v'][:].squeeze()
u = ncf.variables['u'][:].squeeze()
rawlon = ncf.variables['lon'][:].squeeze() # degrees E
rawlat = ncf.variables['lat'][:].squeeze() # degrees N
#
# Array Initailization
numdates = len(u[:,0,0])
lonlen =  len(rawlon)
latlen =  len(rawlat)
lat = np.zeros([latlen,lonlen])
lon = np.zeros([latlen,lonlen])
uavg = np.zeros([numdates-1,latlen,lonlen])
vavg = np.zeros([numdates-1,latlen,lonlen])
#
# Getting lat/lon's in i,j arrays
for i in range(0,lonlen):
    lon[:,i] = rawlon[i]
for j in range(0,latlen):
    lat[j,:] = rawlat[j]
#
# Average Velocities
for k in range(0,numdates-1):
    for i in range(0,lonlen):
        for j in range(0,latlen):
            uavg[k,j,i] = (u[k,j,i] + u[k+1,j,i])/2
            vavg[k,j,i] = (v[k,j,i] + v[k+1,j,i])/2
#
# DOY of 5 day avg Oscar data
centerdates = np.zeros(numdates)
for k in range(0,numdates):
    centerdates[k] = 275 + datadates[k]
#
# ===================================================================================
# Daily Averaged OSCAR Data
# ===================================================================================
#
# Daily Averaged OSCAR Data Array Initialization
oscardates = np.linspace(centerdates[0],centerdates[numdates-1],(centerdates[numdates-1]-centerdates[0])+1)
#
udaily = np.zeros([len(oscardates)-5,latlen,lonlen])
vdaily = np.zeros([len(oscardates)-5,latlen,lonlen])
#
for ii in range(0,len(oscardates)-5):
    #
    # Get the closest OSCAR data date to the current day in the loop. k is index of the closest centerdate
    closestdate, k = find_nearest(centerdates, oscardates[ii])
    #
    # If the closest data day and the current day in the loop are the same, the current day take 100% of the u and v velocities from the oscar data day
    if (closestdate == oscardates[ii]):
        udaily[ii,:,:] = u[k,:,:]
        vdaily[ii,:,:] = v[k,:,:]
    #
    # What if k = 0?  We must "weigh" the u and v currents on the current day of the loop between centerdates[k=0] and centerdates[(k=0) +1] as k-1 doesn't exist 
    elif (k == 0):
        #
        # next closest date must be at k+1
        nextclosestdate = centerdates[k+1]
        kk = k + 1
        #
        # The u and v components at the current day in the loop are weighed by how close they are to the oscar data dates.
        udaily[ii,:,:] = u[k,:,:]*(abs(oscardates[ii] - nextclosestdate)/abs(nextclosestdate-closestdate)) + u[kk,:,:]*(abs(oscardates[ii] - closestdate)/abs(nextclosestdate-closestdate))
        vdaily[ii,:,:] = v[k,:,:]*(abs(oscardates[ii] - nextclosestdate)/abs(nextclosestdate-closestdate)) + v[kk,:,:]*(abs(oscardates[ii] - closestdate)/abs(nextclosestdate-closestdate))
    #
    # If the dates do not match exactly, and k doesn't equal 0 then we need to figure out if the next closest oscar data dates are before or after centerdates[k]
    else:
        #
        # Next closest day to the current day of the loop is AFTER centerdates[k]
        if ((abs(oscardates[ii]-centerdates[k+1]))<=(abs(oscardates[ii]-centerdates[k-1]))):
            nextclosestdate = centerdates[k+1]
            kk = k+1
        #
        # Next closest day to the current day of the loop is BEFORE centerdates[k]
        else: 
            nextclosestdate = centerdates[k-1]
            kk = k-1
        #
        # The u and v components at the current day in the loop are weighed by how close they are to the oscar data dates.
        udaily[ii,:,:] = u[k,:,:]*(abs(oscardates[ii] - nextclosestdate)/abs(nextclosestdate-closestdate)) + u[kk,:,:]*(abs(oscardates[ii] - closestdate)/abs(nextclosestdate-closestdate))
        vdaily[ii,:,:] = v[k,:,:]*(abs(oscardates[ii] - nextclosestdate)/abs(nextclosestdate-closestdate)) + v[kk,:,:]*(abs(oscardates[ii] - closestdate)/abs(nextclosestdate-closestdate))
#
# ===================================================================================
# Mean Kinetic Energy
# ===================================================================================
#
import math as ma
mke5day = np.zeros([numdates-1,latlen,lonlen])
mke = np.zeros([len(oscardates)-5,latlen,lonlen])
#
# MKE from Raw OSCAR Data
for k in range(0,numdates-1):
    for i in range(0,lonlen):
        for j in range(0,latlen):
            mke5day[k,j,i] = 0.5*((u[k,j,i]*100)**2 + (v[k,j,i]*100)**2)
#
# MKE for Daily Weighed OSCAR Currents
for k in range(0,len(oscardates)-5):
    for i in range(0,lonlen):
        for j in range(0,latlen):
            mke[k,j,i] = 0.5*((udaily[k,j,i]*100)**2 + (vdaily[k,j,i]*100)**2)
#
# ===================================================================================
# Find Nearest Lat's and Lon's to float location for U and V components
# ===================================================================================
#
latvalue1 = np.zeros([arraylen])
latvalue2 = np.zeros([arraylen])
lonvalue1 = np.zeros([arraylen])
lonvalue2 = np.zeros([arraylen])
#
latindex1 = np.zeros([arraylen])
latindex2 = np.zeros([arraylen])
lonindex1 = np.zeros([arraylen])
lonindex2 = np.zeros([arraylen])
#
for jj in range(0,arraylen):
    #
    # Get the closest OSCAR data latitude with respect to the float latitude
    latvalue1[jj], latindex1[jj] = find_nearest(lat[:,0], floatlat[jj])
    #
    # Find the next closest OSCAR data latitude with respect to the float latitude
    upperlat = lat[latindex1[jj]+1,0]
    lowerlat = lat[latindex1[jj]-1,0]
    #
    # Closest lat + 0.33 degrees?
    if (abs(floatlat[jj]-upperlat)<=abs(floatlat[jj]-lowerlat)):
        latindex2[jj] = latindex1[jj] + 1
    #
    # Closest lat - 0.33 degrees?
    else:
        latindex2[jj] = latindex1[jj] - 1
    latvalue2[jj] = lat[latindex2[jj],0]
    #
    # Get the closest OSCAR data longitude with respect to the float longitude
    lonvalue1[jj], lonindex1[jj] = find_nearest(lon[0,:], floatlon[jj])
    #
    # Find the next closest OSCAR data longitude with respect to the float longitude
    upperlon = lon[0,lonindex1[jj]+1]
    lowerlon = lon[0,lonindex1[jj]-1]
    #
    # Closest lon + 0.33 degrees?
    if (abs(floatlon[jj]-upperlon)<=abs(floatlon[jj]-lowerlon)):
        lonindex2[jj] = lonindex1[jj] + 1
    #
    # Closest lon - 0.33 degrees?
    else:
        lonindex2[jj] = lonindex1[jj] - 1
    lonvalue2[jj] = lon[0,lonindex2[jj]]
#
# ===================================================================================
# Original 5 Day Bi-Linear Interpolation
# ===================================================================================
#
# Initializing Arrays for Interpolation
floatU5day = np.zeros([arraylen])
floatV5day = np.zeros([arraylen])
floatMKE5day = np.zeros([arraylen])
#
# Loop through the float data
for jj in range(0,arraylen):
    #
    # Loop though the daily oscar data
    for k in range(0,numdates-3):
        if (centerdates[k]-2 <= floatdate[jj] <= centerdates[k] +2):
            #
            # Interpolate Linearly Along Closest Latitude (latvalue1)
            x1_closest5day  = abs(pos2dist([lonvalue1[jj],latvalue1[jj],floatlon[jj],latvalue1[jj]]))
            x2_closest5day  = abs(pos2dist([lonvalue2[jj],latvalue1[jj],floatlon[jj],latvalue1[jj]]))

            xtot_closest5day  = abs(pos2dist([lonvalue1[jj],latvalue1[jj],lonvalue2[jj],latvalue1[jj]]))
            Uinterp_closest5day  = u[k,latindex1[jj], lonindex1[jj]]*(x2_closest5day /xtot_closest5day ) + u[k, latindex1[jj], lonindex2[jj]]*(x1_closest5day /xtot_closest5day)
            Vinterp_closest5day  = v[k,latindex1[jj], lonindex1[jj]]*(x2_closest5day /xtot_closest5day ) + v[k, latindex1[jj], lonindex2[jj]]*(x1_closest5day /xtot_closest5day)
            MKEinterp_closest5day  = mke5day[k,latindex1[jj], lonindex1[jj]]*(x2_closest5day /xtot_closest5day ) + mke5day[k, latindex1[jj], lonindex2[jj]]*(x1_closest5day /xtot_closest5day)
            #
            # Interpolate Linearly Along Furthest Latitude (latvalue2)
            x1_furthest5day  = abs(pos2dist([lonvalue1[jj],latvalue2[jj],floatlon[jj],latvalue2[jj]]))
            x2_furthest5day  = abs(pos2dist([lonvalue2[jj],latvalue2[jj],floatlon[jj],latvalue2[jj]]))
            xtot_furthest5day  = abs(pos2dist([lonvalue1[jj],latvalue2[jj],lonvalue2[jj],latvalue2[jj]]))
            Uinterp_furthest5day  = udaily[k,latindex2[jj], lonindex1[jj]]*(x2_furthest5day /xtot_furthest5day ) + u[k, latindex2[jj], lonindex2[jj]]*(x1_furthest5day /xtot_furthest5day)
            Vinterp_furthest5day  = vdaily[k,latindex2[jj], lonindex1[jj]]*(x2_furthest5day /xtot_furthest5day ) + v[k, latindex2[jj], lonindex2[jj]]*(x1_furthest5day /xtot_furthest5day)
            MKEinterp_furthest5day  = mke5day[k,latindex2[jj], lonindex1[jj]]*(x2_furthest5day /xtot_furthest5day ) + mke5day[k, latindex2[jj], lonindex2[jj]]*(x1_furthest5day /xtot_furthest5day)
            #
            # Tnterpolate Along Londitude
            y_closest5day  = abs(pos2dist([floatlon[jj],latvalue1[jj],floatlon[jj],floatlat[jj]]))
            y_furthest5day  = abs(pos2dist([floatlon[jj],latvalue2[jj],floatlon[jj],floatlat[jj]]))
            ytot5day  = abs(pos2dist([floatlon[jj],latvalue1[jj],floatlon[jj],latvalue2[jj]]))
            floatU5day [jj]= Uinterp_furthest5day *(y_closest5day /ytot5day ) + Uinterp_closest5day *(y_furthest5day /ytot5day )            
            floatV5day [jj]= Vinterp_furthest5day *(y_closest5day /ytot5day ) + Vinterp_closest5day *(y_furthest5day /ytot5day ) 
            floatMKE5day [jj]= MKEinterp_furthest5day *(y_closest5day /ytot5day ) + MKEinterp_closest5day *(y_furthest5day /ytot5day )
# ===================================================================================
# Daily Bi-Linear Interpolation
# ===================================================================================
#
# Initializing Arrays for Interpolation
floatU = np.zeros([arraylen])
floatV = np.zeros([arraylen])
floatMKE = np.zeros([arraylen])
#
# Loop through the float data
for jj in range(0,arraylen):
    #
    # Loop though the daily oscar data
    for k in range(0,len(oscardates)-5):
        #
        # Loop through the daily oscar data ates to find where it equals the float data dates
        if (oscardates[k] == floatdate[jj]):
            #
            # Interpolate Linearly Along Closest Latitude (latvalue1)
            x1_closest = abs(pos2dist([lonvalue1[jj],latvalue1[jj],floatlon[jj],latvalue1[jj]]))
            x2_closest = abs(pos2dist([lonvalue2[jj],latvalue1[jj],floatlon[jj],latvalue1[jj]]))

            xtot_closest = abs(pos2dist([lonvalue1[jj],latvalue1[jj],lonvalue2[jj],latvalue1[jj]]))
            Uinterp_closest = udaily[k,latindex1[jj], lonindex1[jj]]*(x2_closest/xtot_closest) + udaily[k, latindex1[jj], lonindex2[jj]]*(x1_closest/xtot_closest)
            Vinterp_closest = vdaily[k,latindex1[jj], lonindex1[jj]]*(x2_closest/xtot_closest) + vdaily[k, latindex1[jj], lonindex2[jj]]*(x1_closest/xtot_closest)
            MKEinterp_closest = mke[k,latindex1[jj], lonindex1[jj]]*(x2_closest/xtot_closest) + mke[k, latindex1[jj], lonindex2[jj]]*(x1_closest/xtot_closest)
            #
            # Interpolate Linearly Along Furthest Latitude (latvalue2)
            x1_furthest = abs(pos2dist([lonvalue1[jj],latvalue2[jj],floatlon[jj],latvalue2[jj]]))
            x2_furthest = abs(pos2dist([lonvalue2[jj],latvalue2[jj],floatlon[jj],latvalue2[jj]]))
            xtot_furthest = abs(pos2dist([lonvalue1[jj],latvalue2[jj],lonvalue2[jj],latvalue2[jj]]))
            Uinterp_furthest = udaily[k,latindex2[jj], lonindex1[jj]]*(x2_furthest/xtot_furthest) + udaily[k, latindex2[jj], lonindex2[jj]]*(x1_furthest/xtot_furthest)
            Vinterp_furthest = vdaily[k,latindex2[jj], lonindex1[jj]]*(x2_furthest/xtot_furthest) + vdaily[k, latindex2[jj], lonindex2[jj]]*(x1_furthest/xtot_furthest)
            MKEinterp_furthest = mke[k,latindex2[jj], lonindex1[jj]]*(x2_furthest/xtot_furthest) + mke[k, latindex2[jj], lonindex2[jj]]*(x1_furthest/xtot_furthest)
            #
            # Tnterpolate Along Londitude
            y_closest = abs(pos2dist([floatlon[jj],latvalue1[jj],floatlon[jj],floatlat[jj]]))
            y_furthest = abs(pos2dist([floatlon[jj],latvalue2[jj],floatlon[jj],floatlat[jj]]))
            ytot = abs(pos2dist([floatlon[jj],latvalue1[jj],floatlon[jj],latvalue2[jj]]))
            floatU[jj]= Uinterp_furthest*(y_closest/ytot) + Uinterp_closest*(y_furthest/ytot)            
            floatV[jj]= Vinterp_furthest*(y_closest/ytot) + Vinterp_closest*(y_furthest/ytot) 
            floatMKE[jj]= MKEinterp_furthest*(y_closest/ytot) + MKEinterp_closest*(y_furthest/ytot) 
#
# ===================================================================================
# The Tracer U and V components for each float surfacing
# ===================================================================================
#
# Initialize Arrays 
tracerU_kmday = np.zeros([arraylen])
tracerV_kmday = np.zeros([arraylen])
tracerX_km = np.zeros([arraylen])
tracerY_km = np.zeros([arraylen])
for jj in range(0,arraylen):
    #
    # U and V Velocitues in km/day
    tracerU_kmday[jj] = floatU[jj]*86400/1000
    tracerV_kmday[jj] = floatV[jj]*86400/1000
    #
    # Distance Tracer will go in 2 days
    tracerX_km[jj] = tracerU_kmday[jj]*2
    tracerY_km[jj] = tracerV_kmday[jj]*2
#
# ===========================================================
# Tracer bearing and magnitude from float point
# ===========================================================
#
# Initialize Arrays
bearing_U_V = np.zeros([arraylen])
bearing_rads_U_V = np.zeros([arraylen])
magnitude_U_V = np.zeros([arraylen])
tracer_newlat_U_V = np.zeros([arraylen])
tracer_newlon_U_V = np.zeros([arraylen])
#
for jj in range(0,arraylen):
    #
    # bearing and magnitude U,V
    bearing_rads_U_V[jj], bearing_U_V[jj], magnitude_U_V[jj] = find_bearing_magnitude(tracerX_km[jj], tracerY_km[jj])
#
for jj in range(0,arraylen):
    #
    # Call Function to Calculate Tracer Locations
    tracer_newlat_U_V[jj],              tracer_newlon_U_V[jj] = destpoint(ma.radians(floatlat[jj]),ma.radians(floatlon[jj]),ma.radians(bearing_U_V[jj]), magnitude_U_V[jj])
# ===========================================================
# Tracer bearing and magnitude
# ===========================================================
#
origtotracer = np.zeros([arraylen-1])
origtofloat = np.zeros([arraylen-1])
tracertofloat = np.zeros([arraylen-1])
tracerfloatratio = np.zeros([arraylen-1])
#
for jj in range(0,arraylen-1):
    origtotracer[jj] = pos2dist([floatlon[jj],floatlat[jj], tracer_newlon_U_V[jj], tracer_newlat_U_V[jj]])
    origtofloat[jj] = pos2dist([floatlon[jj],floatlat[jj], floatlon[jj+1],floatlat[jj+1]])
    tracertofloat[jj] = pos2dist([floatlon[jj+1],floatlat[jj+1], tracer_newlon_U_V[jj], tracer_newlat_U_V[jj]])
    tracerfloatratio[jj] = abs(tracertofloat[jj])/origtotracer[jj]
#
# Do you want to plot Gridded Ws?  Yes: gridWplt = 1, No: else
gridWdat = 1
if gridWdat == 1:
    #
    # ===========================================================
    # Read the Float Data and Interpolate for Gridding
    # ===========================================================
    #
    os.chdir('/data/orbprocess_mail/alex/Oct02_Dec25_2013/data/type/')
    Ygrid = np.linspace(-1800,-10,1791*2)
    interplength = len(Ygrid)
    GridData = np.zeros([interplength,7,arraylen])
    ChlyDumby = np.zeros([interplength,arraylen])
    for i in range(0,arraylen):
        dataout, rows, cols = csvread(fullnames[i])
        #
        TemInterp = interpolate.interp1d(dataout[0,:], dataout[1,:],kind='linear')
        DenInterp = interpolate.interp1d(dataout[0,:], dataout[2,:],kind='linear')
        SalInterp = interpolate.interp1d(dataout[0,:], dataout[3,:],kind='linear')
        PreInterp = interpolate.interp1d(dataout[0,:], dataout[4,:],kind='linear')
        OxyInterp = interpolate.interp1d(dataout[0,:], dataout[7,:],kind='linear')
        ChlInterp = interpolate.interp1d(dataout[0,:], dataout[11,:],kind='linear')
        BskInterp = interpolate.interp1d(dataout[0,:], dataout[12,:],kind='linear')
        CDMInterp = interpolate.interp1d(dataout[0,:], dataout[13,:],kind='linear')
        #
        for j in range(0,interplength):
            GridData[j,0,i] = TemInterp(Ygrid[j])
            GridData[j,1,i] = DenInterp(Ygrid[j])
            GridData[j,2,i] = SalInterp(Ygrid[j])
            GridData[j,3,i] = OxyInterp(Ygrid[j])            
            ChlyDumby[j,i] = ChlInterp(Ygrid[j])    
            GridData[j,5,i] = BskInterp(Ygrid[j])   
            GridData[j,6,i] = CDMInterp(Ygrid[j])   
            if (ChlInterp(Ygrid[j]) <= 0.0):
                GridData[j,4,i] = GridData[j-1,4,i]
            else: 
                GridData[j,4,i] = ChlyDumby[j,i]    
    #
    # ===========================================================
    # Find Gridded Change Velocities Based on Nearest Values
    # ===========================================================
    #
    GridWData = np.zeros([interplength,7,arraylen-1])
    GridWDepths = np.zeros([interplength,7,arraylen-1])
    GridWDates = np.zeros([interplength,7,arraylen-1])
    for k in range(0,7):
        for i in range (0,arraylen-1):
            for j in range(0,interplength):
                value, index = find_nearest(GridData[:,k,i+1], GridData[j,k,i])
                startdepth = Ygrid[j]
                newdepth = Ygrid[index]
                GridWDepths[j,k,i] = (newdepth + startdepth)/2
                dz = newdepth-startdepth
                GridWData[j,k,i] = dz/172800
                GridWDates[j,k,i] = (floatdate[i] + floatdate[i+1])/2
                # checks.
    #
    # ===========================================================
    # Find Weekly Avg Gridded Vertical Velocities
    # ===========================================================
    #
    GridWDataAvg = np.zeros([interplength,7,arraylen-3])
    GridWDepthsAvg = np.zeros([interplength,7,arraylen-3])
    GridWDatesAvg = np.zeros([interplength,7,arraylen-3])
    for k in range(0,7):
        for i in range (0,arraylen-3):
            for j in range(0,interplength):
                GridWDataAvg[j,k,i] = (GridWData[j,k,i] + GridWData[j,k,i+1] + GridWData[j,k,i+2])/3
                GridWDepthsAvg[j,k,i] = (GridWDepths[j,k,i] + GridWDepths[j,k,i+1] + GridWDepths[j,k,i+2])/3
                GridWDatesAvg[j,k,i] = (GridWDates[j,k,i] + GridWDates[j,k,i+1] + GridWDates[j,k,i+2])/3
    #
    # ===========================================================
    # Find Depth of Isochly Contours
    # ===========================================================
    #
    isobiolines = np.linspace(.1,.40,20)
    numbiolines = len(isobiolines)
    Biodepths = np.zeros([numbiolines,arraylen])
    for i in range(0,arraylen):
        #print i
        for j in range(0,numbiolines):
            #print j
            count = 0
            while (GridData[count,4,i] < isobiolines[j]):
                    count = count + 1 
            Biodepths[j,i] = Ygrid[count]
    #
    # ===========================================================
    # Find Depth of Isocdom Contours
    # ===========================================================
    #
    isocdomlines = np.linspace(1.7,2.2,21)
    numcdomlines = len(isocdomlines)
    cdomdepths = np.zeros([numcdomlines,arraylen])
    for i in range(0,arraylen):
        #print i
        for j in range(0,numcdomlines):
            #print j
            count = 0
            while (GridData[count,6,i] > isocdomlines[j]):
                    count = count + 1 
            cdomdepths[j,i] = Ygrid[count]

    #
    # ===========================================================
    # Find Isochly Contour Change Velocities
    # ===========================================================
    #
    IsoBioWData = np.zeros([numbiolines,arraylen-1])
    IsoBioWDepths = np.zeros([numbiolines,arraylen-1])
    IsoBioWDates = np.zeros([numbiolines,arraylen-1])
    IsoBioFake = np.zeros([numbiolines,arraylen-1])
    for i in range(0,arraylen-1):
        for k in range(0,numbiolines):
            j = i
            IsoBioWData[k,i] = (Biodepths[k,j+1] - Biodepths[k,j])/172800
            IsoBioFake[k,i] = 0
            IsoBioWDepths[k,i] = (Biodepths[k,j+1] + Biodepths[k,j])/2
            IsoBioWDates[k,i] = (floatdate[j+1] + floatdate[j])/2
    #
    # ===========================================================
    # Find Isochly Contour Change Velocities
    # ===========================================================
    #
    IsoCdomWData = np.zeros([numcdomlines,arraylen-1])
    IsoCdomWDepths = np.zeros([numcdomlines,arraylen-1])
    IsoCdomWDates = np.zeros([numcdomlines,arraylen-1])
    IsoCdomFake = np.zeros([numcdomlines,arraylen-1])
    for i in range(0,arraylen-1):
        for k in range(0,numcdomlines):
            j = i
            IsoCdomWData[k,i] = (cdomdepths[k,j+1] - cdomdepths[k,j])/172800
            IsoCdomFake[k,i] = 0
            IsoCdomWDepths[k,i] = (cdomdepths[k,j+1] + cdomdepths[k,j])/2
            IsoCdomWDates[k,i] = (floatdate[j+1] + floatdate[j])/2
    #
    # ===========================================================
    # Find Avg Isoline Change Velocities
    # ===========================================================
    #
    IsoBioWDataAvg = np.zeros([numbiolines,arraylen-3])
    IsoBioWDepthsAvg = np.zeros([numbiolines,arraylen-3])
    IsoBioWDatesAvg = np.zeros([numbiolines,arraylen-3])
    IsoBioFakeAvg = np.zeros([numbiolines,arraylen-3])
    for i in range (0,arraylen-3):
        for j in range(0,numbiolines):
                IsoBioWDataAvg[j,i] = (IsoBioWData[j,i] + IsoBioWData[j,i+1] + IsoBioWData[j,i+2])/3
                IsoBioWDepthsAvg[j,i] = (IsoBioWDepths[j,i] + IsoBioWDepths[j,i+1] + IsoBioWDepths[j,i+2])/3
                IsoBioWDatesAvg[j,i] = (IsoBioWDates[j,i] + IsoBioWDates[j,i+1] + IsoBioWDates[j,i+2])/3
                IsoBioFakeAvg[j,i] = (IsoBioFake[j,i] + IsoBioFake[j,i+1] + IsoBioFake[j,i+2])/3
    #
    # ===========================================================
    # Find Avg IsoCDOM Change Velocities
    # ===========================================================
    #
    IsoCdomWDataAvg = np.zeros([numcdomlines,arraylen-3])
    IsoCdomWDepthsAvg = np.zeros([numcdomlines,arraylen-3])
    IsoCdomWDatesAvg = np.zeros([numcdomlines,arraylen-3])
    IsoCdomFakeAvg = np.zeros([numcdomlines,arraylen-3])
    for i in range (0,arraylen-3):
        for j in range(0,numcdomlines):
                IsoCdomWDataAvg[j,i] = (IsoCdomWData[j,i] + IsoCdomWData[j,i+1] + IsoCdomWData[j,i+2])/3
                IsoCdomWDepthsAvg[j,i] = (IsoCdomWDepths[j,i] + IsoCdomWDepths[j,i+1] + IsoCdomWDepths[j,i+2])/3
                IsoCdomWDatesAvg[j,i] = (IsoCdomWDates[j,i] + IsoCdomWDates[j,i+1] + IsoCdomWDates[j,i+2])/3
                IsoCdomFakeAvg[j,i] = (IsoCdomFake[j,i] + IsoCdomFake[j,i+1] + IsoCdomFake[j,i+2])/3
    #
    # ===========================================================
    # Find Depth Averaged Avg Isoline Change Velocities
    # ===========================================================
    #
    IsoBioWDataDepthAvgAvg = np.zeros(arraylen-3)
    IsoBioWDepthsDepthAvgAvg = np.zeros(arraylen-3)
    IsoBioWDatesDepthAvgAvg = np.zeros(arraylen-3)
    IsoBioFakeDepthAvgAvg = np.zeros(arraylen-3)
    for i in range (0,arraylen-3):
        IsoBioWDataDepthAvgAvg[i] = np.mean(IsoBioWDataAvg[:,i])
        IsoBioWDepthsDepthAvgAvg[i] = np.mean(IsoBioWDepthsAvg[:,i])
        IsoBioWDatesDepthAvgAvg[i] = np.mean(IsoBioWDatesAvg[:,i])
        IsoBioFakeDepthAvgAvg[i] = np.mean(IsoBioFakeAvg[:,i])
    #
    # ===========================================================
    # Find Depth Averaged Avg Isoline Change Velocities
    # ===========================================================
    #
    IsoCdomWDataDepthAvgAvg = np.zeros(arraylen-3)
    IsoCdomWDepthsDepthAvgAvg = np.zeros(arraylen-3)
    IsoCdomWDatesDepthAvgAvg = np.zeros(arraylen-3)
    IsoCdomFakeDepthAvgAvg = np.zeros(arraylen-3)
    for i in range (0,arraylen-3):
        IsoCdomWDataDepthAvgAvg[i] = np.mean(IsoCdomWDataAvg[:,i])
        IsoCdomWDepthsDepthAvgAvg[i] = np.mean(IsoCdomWDepthsAvg[:,i])
        IsoCdomWDatesDepthAvgAvg[i] = np.mean(IsoCdomWDatesAvg[:,i])
        IsoCdomFakeDepthAvgAvg[i] = np.mean(IsoCdomFakeAvg[:,i])
    #
    # ===========================================================
    # Find Sink Vel
    # ===========================================================
    #
    NearestWGridded  = np.zeros([numbiolines,arraylen-1])
    SinkWData = np.zeros([numbiolines,arraylen-1])
    for i in range (0,arraylen-1):
        for k in range(0,numbiolines):
            value, index = find_nearest(GridWDepths[:,1,i], IsoBioWDepths[k,i])
            SinkWData[k,i] = IsoBioWData[k,i] - GridWData[index,1,i]
            NearestWGridded[k,i] = GridWData[index,1,i]
    #
    # ===========================================================
    # Find Avg Sink Vel
    # ===========================================================
    #
    SinkWDataAvg = np.zeros([numbiolines,arraylen-3])
    NearestWGriddedAvg  = np.zeros([numbiolines,arraylen-3])
    for i in range (0,arraylen-3):
        for k in range(0,numbiolines):
            value, index = find_nearest(GridWDepthsAvg[:,1,i], IsoBioWDepthsAvg[k,i])
            SinkWDataAvg[k,i] = IsoBioWDataAvg[k,i] - GridWDataAvg[index,1,i]
            NearestWGriddedAvg[k,i] = GridWDataAvg[index,1,i]
    #
    # ===========================================================
    # Find Depth Averaged Avg Sink Vel
    # ===========================================================
    #
    SinkWDataDepthAvgAvg = np.zeros(arraylen-3)
    NearestWGriddedAvgAvg  = np.zeros(arraylen-3)
    for i in range (0,arraylen-3):
        SinkWDataDepthAvgAvg[i] = np.mean(SinkWDataAvg[:,i])
        NearestWGriddedAvgAvg[i] = np.mean(NearestWGriddedAvg[:,i])
    #
    # ===========================================================
    # Find Average Chly Concentration in the upper 30 m
    # ===========================================================
    #
    chlyavg30meter = np.zeros(arraylen-6)
    for i in range(0,arraylen-6):
            dataout2, rows2, cols2 = csvread(fullnames[i])
            gdataout2, grows2, gcols = csvread(gradnames[i])
            depthlength = len (dataout2[0,:])
            counter = 0
            chlytotal = 0
            for ii in range(0,depthlength):
                if (abs(dataout2[0,ii]) < 30.0):
                    counter = counter + 1
                    chlytotal = chlytotal + dataout2[11,ii]
            chlyavg30meter[i] = chlytotal/counter   
    #
    # ===========================================================
    # Demean Average Chly Concentration in the upper 30 m; MKE
    # ===========================================================
    #   
    # Mean chl and MKE over the entire time series
    meanchlyavg30meter = np.mean(chlyavg30meter)
    meanfloatMKE = np.mean(floatMKE)
    #
    # Initialization for demean data
    demeaned_chlyavg30meter = np.zeros(arraylen-6)
    demeaned_floatMKE = np.zeros(arraylen-6)
    #
    # Initialization for Abreviated MKE, Chly time series
    floatdateabrv = np.zeros(arraylen-6)    
    floatMKEabrv = np.zeros(arraylen-6)
    chlyavg30meterabr = np.zeros(arraylen-6)
    #
    # Demean loop
    for i in range(0,arraylen-6):
        demeaned_chlyavg30meter[i] = chlyavg30meter[i]-meanchlyavg30meter
        demeaned_floatMKE[i] = floatMKE[i] - meanfloatMKE
        floatdateabrv[i] = floatdate[i]
        floatMKEabrv[i] = floatMKE[i]
        chlyavg30meterabr[i] = chlyavg30meter[i]
    #
    # ===========================================================
    # Relative Change in Demean, Normal Average Chly Concentration in the upper 30 m; MKE
    # ===========================================================
    #  
    # Initialization Demean Change
    delta_demeaned_chlyavg30meter = np.zeros(arraylen-7) 
    delta_demeaned_floatMKE = np.zeros(arraylen-7)
    #
    # initialization Change
    delta_floatdate = np.zeros(arraylen-7) 
    delta_chlyavg30meter = np.zeros(arraylen-7)  
    delta_floatMKE = np.zeros(arraylen-7)
    relative_chlychange = np.zeros(arraylen-7)
    relative_MKEchange = np.zeros(arraylen-7)
    for i in range(0,arraylen-7):
        #
        # Demean Relative Change
        delta_demeaned_chlyavg30meter[i] = (demeaned_chlyavg30meter[i+1] - demeaned_chlyavg30meter[i])/(floatdate[i+1]-floatdate[i])
        delta_demeaned_floatMKE[i] = (demeaned_floatMKE[i+1] - demeaned_floatMKE[i])/(floatdate[i+1]-floatdate[i])
        delta_floatdate[i] = (floatdate[i] + floatdate[i+1])/2
        #
        # Relative change
        delta_chlyavg30meter[i] = (chlyavg30meter[i+1] - chlyavg30meter[i])
        delta_floatMKE[i] = (floatMKE[i+1] - floatMKE[i])
        relative_chlychange[i] = delta_chlyavg30meter[i]/chlyavg30meter[i]
        relative_MKEchange[i] = delta_floatMKE[i]/floatMKE[i]
    #
    # ===========================================================
    # drho/dz raw data max
    # ===========================================================
    #  
    dengradmaxes = np.zeros(arraylen)
    dengradmaxes2 = np.zeros(arraylen)

    dengraddepths = np.zeros(arraylen)
    for i in range(0,arraylen):
        gdataout, grows, gcols = csvread('/data/orbprocess_mail/alex/Oct02_Dec25_2013/data/type/' + gradnames[i])
        #gradarray_rowlengths = len()
        #dengradmaxes2[i] = findmaxs(gdataout[2,:], gdataout[1,:], arraylen, grows)
        mx = np.amax(abs(gdataout[2,:]))
        r33 = 0
        for r3 in range(0, grows):
                r33 = 0 + r3
                if (abs(gdataout[2,r3]) == mx):
                        break
        dengradmaxes[i] = gdataout[2,r33]
        dengraddepths[i] = gdataout[0,r33]
        # dengradmaxes[i] = np.amin(gdataout[2,:])
        # value, index = find_nearest(gdataout[2,:],dengradmaxes[i])
        # dengraddepths[i] = gdataout[0,index]


#print dengraddepths
#print dengradmaxes
#print dengradmaxes2



os.chdir('/home/ardavies/satdata/OSCAR/pdfoutput')



pp1 = PdfPages('linear_density.pdf')
for i in range(0,arraylen):
    dataout, rows, cols = csvread('/data/orbprocess_mail/alex/Oct02_Dec25_2013/data/type/' + fullnames[i])

    # ===========================================================
    # drho/dz raw data max
    # ===========================================================
    #
    denavg3 = np.zeros(rows-2)
    depthavg3 = np.zeros(rows-2)
    denavg5 = np.zeros(rows-4)
    depthavg5 = np.zeros(rows-4)
    denavg7 = np.zeros(rows-6)
    depthavg7 = np.zeros(rows-6)
    for ii in range(0,rows-2):
        denavg3[ii] = (dataout[2,ii] + dataout[2,ii + 1] + dataout[2,ii + 2])/3
        depthavg3[ii] = (dataout[0,ii] + dataout[0,ii + 1] + dataout[0,ii + 2])/3
    for iii in range(0,rows-4):   
        denavg5[iii] = (dataout[2,iii] + dataout[2,iii + 1] + dataout[2,iii + 2]+ dataout[2,iii + 3]+ dataout[2,iii + 4])/5
        depthavg5[iii] = (dataout[0,iii] + dataout[0,iii + 1] + dataout[0,iii + 2]+ dataout[0,iii + 3]+ dataout[0,iii + 4])/5
    for iiii in range(0,rows-6):
        denavg7[iiii] = (dataout[2,iiii] + dataout[2,iiii + 1] + dataout[2,iiii + 2]+ dataout[2,iiii + 3] + dataout[2,iiii + 4] + dataout[2,iiii + 5] + dataout[2,iiii + 6])/7
        depthavg7[iiii] = (dataout[0,iiii] + dataout[0,iiii + 1] + dataout[0,iiii + 2]+ dataout[0,iiii + 3] + dataout[0,iiii + 4] + dataout[0,iiii + 5] + dataout[0,iiii + 6])/7


    drhodz = np.zeros(rows-1)
    drhodz_depth = np.zeros(rows-1)
    drhodz3 = np.zeros(rows-3)
    drhodz3_depth = np.zeros(rows-3)
    drhodz5 = np.zeros(rows-5)
    drhodz5_depth = np.zeros(rows-5)
    drhodz7 = np.zeros(rows-7)
    drhodz7_depth = np.zeros(rows-7)




    for ii in range(0,rows-1):
        drhodz[ii] = (dataout[2,ii+1]-dataout[2,ii])/(dataout[0,ii+1]-dataout[0,ii])
        drhodz_depth[ii] = (dataout[0,ii]+dataout[0,ii+1])/2
    #if i == 0:
        #print dataout[0,rows-1]
        #print drhodz
    for ii in range(0,rows-3):
        drhodz3[ii] = (denavg3[ii+1]-denavg3[ii])/(depthavg3[ii+1] -depthavg3[ii])
        drhodz3_depth[ii] = (depthavg3[ii] + depthavg3[ii+1])/2
    for ii in range(0,rows-5):
        drhodz5[ii] = (denavg5[ii+1]-denavg5[ii])/(depthavg5[ii+1] -depthavg5[ii])
        drhodz5_depth[ii] = (depthavg5[ii] + depthavg5[ii+1])/2
    for ii in range(0,rows-7):
        drhodz7[ii] = (denavg7[ii+1]-denavg7[ii])/(depthavg7[ii+1] -depthavg7[ii])
        drhodz7_depth[ii] = (depthavg7[ii] + depthavg7[ii+1])/2
    
    # Array Max
    mx1 = np.amin(drhodz)
    mx3 = np.amin(drhodz3)
    mx5 = np.amin(drhodz5)
    mx7 = np.amin(drhodz7)
    # Index of Max
    idx1 = (np.abs(drhodz-mx1)).argmin()
    idx3 = (np.abs(drhodz3-mx3)).argmin()
    idx5 = (np.abs(drhodz5-mx5)).argmin()
    idx7 = (np.abs(drhodz7-mx7)).argmin()
    # Depth of Max
    thedepthofmax1 = drhodz_depth[idx1]
    thedepthofmax3 = drhodz3_depth[idx3]
    thedepthofmax5 = drhodz5_depth[idx5]
    thedepthofmax7 = drhodz7_depth[idx7]
    #
    #
    d2rhodz = np.zeros(rows-2)
    d2rhodz_depth = np.zeros(rows-2)
    d2rhodz3 = np.zeros(rows-4)
    d2rhodz3_depth = np.zeros(rows-4)
    d2rhodz5 = np.zeros(rows-6)
    d2rhodz5_depth = np.zeros(rows-6)
    d2rhodz7 = np.zeros(rows-8)
    d2rhodz7_depth = np.zeros(rows-8)




    for ii in range(0,rows-2):
        d2rhodz[ii] = (drhodz[ii+1]-drhodz[ii])/(drhodz_depth[ii+1]-drhodz_depth[ii])
        d2rhodz_depth[ii] = (drhodz_depth[ii+1]+drhodz_depth[ii])/2
    for ii in range(0,rows-4):
        d2rhodz3[ii] = (drhodz3[ii+1]-drhodz3[ii])/(drhodz3_depth[ii+1] -drhodz3_depth[ii])
        d2rhodz3_depth[ii] = (drhodz3_depth[ii] + drhodz3_depth[ii+1])/2
    for ii in range(0,rows-6):
        d2rhodz5[ii] = (drhodz5[ii+1]-drhodz5[ii])/(drhodz5_depth[ii+1] -drhodz5_depth[ii])
        d2rhodz5_depth[ii] = (drhodz5_depth[ii] + drhodz5_depth[ii+1])/2
    for ii in range(0,rows-8):
        d2rhodz7[ii] = (drhodz7[ii+1]-drhodz7[ii])/(drhodz7_depth[ii+1] -drhodz7_depth[ii])
        d2rhodz7_depth[ii] = (drhodz7_depth[ii] + drhodz7_depth[ii+1])/2
    

    idx21 = (np.abs(d2rhodz-0.0)).argmin()
    idx23 = (np.abs(d2rhodz3-0.0)).argmin()
    idx25 = (np.abs(d2rhodz5-0.0)).argmin()
    idx27 = (np.abs(d2rhodz7-0.0)).argmin()


    thedepthofmax21 = d2rhodz_depth[idx1]
    thedepthofmax23 = d2rhodz3_depth[idx3]
    thedepthofmax25 = d2rhodz5_depth[idx5]
    thedepthofmax27 = d2rhodz7_depth[idx7]

    print thedepthofmax25

    # How many plots? Which axis
    fig = plt.figure()
    host = fig.add_subplot(111)
    par1 = host.twiny()

    from pylab import *
    #
    #
    # What are you plotting?
    thedepthofmax = dengraddepths[i]
    p1 = host.plot(dataout[2,:],dataout[0,:],'0.75',linewidth = 1)
    host.scatter(dataout[2,:],dataout[0,:],c='0.75')
    host.plot(denavg3,depthavg3,'m',linewidth = 1)
    host.plot(denavg5,depthavg5,'b',linewidth = 1)
    host.plot(denavg7,depthavg7,'k',linewidth = 1)

    #p11, = host.scatter(dengradmaxes[i],dengraddepths[i], c='r')
    p2 = par1.plot(dataout[11,:],dataout[0,:],'g',linewidth = 1)
    #
    # Setting axi
    host.xaxis.set_major_locator(MaxNLocator(5))
    host.set_xlim([1026.65,1027.25])
    host.set_ylim([-400,0])
    #host.axhline(xmin=0.0, xmax=0.18, y=thedepthofmax,linewidth=1, color='r')
    host.axhline(xmin=0.0, xmax=0.16,y=thedepthofmax1,linewidth=1, color='0.75')
    host.axhline(xmin=0.0, xmax=0.14,y=thedepthofmax3,linewidth=1, color='m')
    host.axhline(xmin=0.0, xmax=0.11,y=thedepthofmax5,linewidth=1, color='b')
    host.axhline(xmin=0.0, xmax=0.10,y=thedepthofmax7,linewidth=1, color='k')


    host.axhline(xmin=0.84, xmax=1,y=thedepthofmax21,linewidth=1, color='0.75')
    host.axhline(xmin=0.86, xmax=1,y=thedepthofmax23,linewidth=1, color='m')
    host.axhline(xmin=0.88, xmax=1,y=thedepthofmax25,linewidth=1, color='b')
    host.axhline(xmin=0.9, xmax=1,y=thedepthofmax27,linewidth=1, color='k')


    par1.set_xlim([0,6])
    #par1.set_yscale('log')
    #par1.yaxis.set_ylim([-.12, .12])
    #
    # axis labels
    host.set_xlabel("Density")
    host.set_ylabel("Depth")
    par1.set_ylabel("Chlorophyll Fluorescence (micro g/l)")
    #
    #
    # colors
    #host.yaxis.label.set_color(p1.get_color())
    #par1.yaxis.label.set_color(p2.get_color())
    #
    # Using the colors and setting other things about plot
    #tkw = dict(size=4, width=1.5)
    #host.tick_params(axis='y', colors=p1.get_color(), **tkw)
    #par1.tick_params(axis='y', colors=p2.get_color(), **tkw)
    #host.tick_params(axis='x2_closest5day', **tkw)
    #lines = [p1, p2]
    fig.suptitle("Date: " + str(floatdate[i])) 
    fig.savefig(pp1, format = "pdf")  
pp1.close()
os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/linear_density.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')





kjhjbhj












    #






import numpy, scipy, pylab, random  

pp1 = PdfPages('densityfft.pdf')
for i in range(0,arraylen):
    dataout, rows, cols = csvread('/data/orbprocess_mail/alex/Oct02_Dec25_2013/data/type/' + fullnames[i])

    # #
    # #chdata = np.zeros(rows)
    # for ii in range(0,rows):
    #     dataout[11,ii] = np.log(dataout[11,ii])
    # TemInterp = interpolate.interp1d(dataout[0,:], dataout[1,:],kind='linear')
    # DenInterp = interpolate.interp1d(dataout[0,:], dataout[2,:],kind='linear')
    # SalInterp = interpolate.interp1d(dataout[0,:], dataout[3,:],kind='linear')
    # PreInterp = interpolate.interp1d(dataout[0,:], dataout[4,:],kind='linear')
    # OxyInterp = interpolate.interp1d(dataout[0,:], dataout[7,:],kind='linear')
    # ChlInterp = interpolate.interp1d(dataout[0,:], dataout[11,:],kind='linear')
    # BskInterp = interpolate.interp1d(dataout[0,:], dataout[12,:],kind='linear')
    # CDMInterp = interpolate.interp1d(dataout[0,:], dataout[13,:],kind='linear')
    # #
    # for j in range(0,interplength):
    #     GridData[j,0,i] = TemInterp(Ygrid[j])
    #     GridData[j,1,i] = DenInterp(Ygrid[j])
    #     GridData[j,2,i] = SalInterp(Ygrid[j])
    #     GridData[j,3,i] = PreInterp(Ygrid[j])   
    #     GridData[j,4,i] = OxyInterp(Ygrid[j])            
    #     ChlyDumby[j,i] = ChlInterp(Ygrid[j])    
    #     GridData[j,6,i] = BskInterp(Ygrid[j])   
    #     GridData[j,7,i] = CDMInterp(Ygrid[j])   
    #     #if (ChlInterp(Ygrid[j]) < 0):
    #     #   GridData[j,5,i] = GridData[j-1,5,i]
    #     #else: 
    #     GridData[j,5,i] = ChlyDumby[j,i]    

    denfft=scipy.fft(dataout[2,:]) 

    fig = plt.figure()    
    ax = fig.add_subplot(111)
    ax.plot(denfft,'k',label='MKE',linewidth = 3)
    fig.savefig(pp1, format = "pdf")  
pp1.close()
os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/densityfft.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')

aggfh


    
#     #Verification Ploting 
#     fig = plt.figure()
#     ax1 = fig.add_subplot(1,1,1)
#     #p11 = ax1.plot(Ygrid,GridData[:,5,i],'k')
#     #p12 = ax1.plot(Ygrid,ChlyDumby[:,i])
#     #p13 = ax1.scatter(dataout[0,:], dataout[11,:])
#     #ax1.set_xlabel("Depth(m)")
#     #ax1.set_ylabel("Chlorophyll Concentration")

#     #pp11 = ax1.plot(Ygrid,GridData[:,5,i],'k')
#     pp12 = ax1.plot(dataout[2,:],dataout[0,:])
#     #pp13 = ax1.scatter(dataout[2,:],dataout[0,:])
#     pp14 = ax1.scatter(dengradmaxes[i],dengraddepths[i], c='r')
#     ax1.set_xlabel("Depth(m)")
#     ax1.set_xlim()
#     ax1.set_ylim()
#     ax1.set_ylabel("Density")
#     # ax1.set_yscale('log')
#     fig.suptitle("Linear Interpolation Pre Log. Date: " + str(floatdate[i])) 
#     fig.savefig(pp1, format = "pdf")  
# pp1.close()
# os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/linear_density.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')

# dasfa



pp = PdfPages('floatmke-chly-daily-delta.pdf')

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(delta_floatdate,delta_demeaned_floatMKE,'k',label='MKE',linewidth = 3)
ax2 = ax.twinx()
ax2.plot(delta_floatdate,delta_demeaned_chlyavg30meter,'b',label='30 meter Avg Chl Concentration',linewidth = 3)
ax.legend(loc=2)
ax2.legend(loc=0)
ax.grid()
ax.set_xlabel("DOY")
ax.set_ylabel("d(MKE)/dt")
ax2.set_ylabel("d(chl)/dt")
ax.set_ylim([-1200,1200])
ax2.set_ylim([-1.2,1.2])
plt.savefig(pp, format = "pdf")
pp.close()
os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/floatmke-chly-daily-delta.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')


#
# ===========================================================
# Function to find depth of max
# ===========================================================
#
def linearregress(x, y):
    from scipy import stats
    #
    # Linearly Regression
    slope, intercept, rvalue, pvalue, stderr = stats.linregress(x,y)
    #
    # Regression Line
    regressline = slope*x+intercept
    #
    return regressline, rvalue, pvalue, stderr


# Linear Regression
line, r_value, p_value, std_err = linearregress(delta_demeaned_floatMKE[15:],delta_demeaned_chlyavg30meter[15:])

# Regression Plot
pp = PdfPages('floatmke-chly-daily-delta-demean-regression-cut.pdf')
fig = plt.figure()
ax = fig.add_subplot(111)
line1,=ax.plot(delta_demeaned_floatMKE[15:],delta_demeaned_chlyavg30meter[15:],'o', )
ax.set_xlabel("d(MKE)/dt")
ax.set_ylabel("d(30m Avg Chl)/dt")
line2,=ax.plot(delta_demeaned_floatMKE[15:],line,'k-')

fig.suptitle('R value: ' + str(r_value) + ',  P value: ' + str(p_value))
plt.savefig(pp, format = "pdf")
pp.close()
os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/floatmke-chly-daily-delta-demean-regression-cut.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')

# Linear Regression
line, r_value, p_value, std_err = linearregress(relative_MKEchange[:15],relative_chlychange[:15])

# Regression Plot
pp = PdfPages('floatmke-chly-daily-delta-demean-regression-change-cut.pdf')
fig = plt.figure()
ax = fig.add_subplot(111)
line1,=ax.plot(relative_MKEchange[:15],relative_chlychange[:15],'o', )
ax.set_xlabel("d(MKE)/dt")
ax.set_ylabel("d(30m Avg Chl)/dt")
line2,=ax.plot(relative_MKEchange[:15],line,'k-')

fig.suptitle('R value: ' + str(r_value) + ',  P value: ' + str(p_value))
plt.savefig(pp, format = "pdf")
pp.close()
os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/floatmke-chly-daily-delta-demean-regression-change-cut.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')
fghfhjf

gggdfs




pp = PdfPages('floatmke-chly-daily2_secondmission.pdf')
fig = plt.figure()
#fig.subplots_adjust(bottom=0.22)
#
# How many plots? Which axis
host = fig.add_subplot(111)
par1 = host.twinx()

from pylab import *
#
#
# What are you plotting?
p1, = host.plot(floatdate,floatMKE,'k',label='Using Daily Averaged OSCAR Data',linewidth = 3)
p2, = par1.plot(floatdateabrv,chlyavg30meter,'b',label='30 meter Avg Chl Concentration',linewidth = 3)

l = legend(loc = 2)

#
# Setting axi
host.xaxis.set_major_locator(MaxNLocator(5))
host.yaxis.set_major_locator(MaxNLocator(5))
par1.yaxis.set_major_locator(MaxNLocator(5))
#
# axis labels
host.set_xlabel("DOY")
host.set_ylabel("MKE (cm^2/s^2)")
par1.set_ylabel("Chlorophyll Fluorescence (micro g/l)")
#
host.set_yscale('log')
par1.set_yscale('log')
#
# colors
host.yaxis.label.set_color(p1.get_color())
par1.yaxis.label.set_color(p2.get_color())
#
# Using the colors and setting other things about plot
tkw = dict(size=4, width=1.5)
host.tick_params(axis='y', colors=p1.get_color(), **tkw)
par1.tick_params(axis='y', colors=p2.get_color(), **tkw)
host.tick_params(axis='x2_closest5day', **tkw)
lines = [p1, p2]
plt.savefig(pp, format = "pdf")
pp.close()
os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/floatmke-chly-daily2_secondmission.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')





fig = plt.figure()
#fig.subplots_adjust(bottom=0.22)
#
# How many plots? Which axis
host = fig.add_subplot(111)
par1 = host.twinx()

from pylab import *
#
#
# What are you plotting?
p1, = host.plot(floatdateabrv,demeaned_floatMKE,'k',label='Using Daily Averaged OSCAR Data',linewidth = 3)
p2, = par1.plot(floatdateabrv,demeaned_chlyavg30meter,'b',label='30 meter Avg Chl Concentration',linewidth = 3)

l = legend(loc = 2)

#
# Setting axi
host.xaxis.set_major_locator(MaxNLocator(5))
host.yaxis.set_ylim([-1200, 1200])
par1.yaxis.set_ylim([-.12, .12])
#
# axis labels
host.set_xlabel("DOY")
host.set_ylabel("MKE (cm^2/s^2)")
par1.set_ylabel("Chlorophyll Fluorescence (micro g/l)")
#
#
# colors
host.yaxis.label.set_color(p1.get_color())
par1.yaxis.label.set_color(p2.get_color())
#
# Using the colors and setting other things about plot
tkw = dict(size=4, width=1.5)
host.tick_params(axis='y', colors=p1.get_color(), **tkw)
par1.tick_params(axis='y', colors=p2.get_color(), **tkw)
host.tick_params(axis='x2_closest5day', **tkw)
lines = [p1, p2]
plt.savefig(pp, format = "pdf")
pp.close()
os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/floatmke-chly-daily3_secondmission.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')



#
# ===========================================================
# Function to find depth of max
# ===========================================================
#
def linearregress(x, y):
    from scipy import stats
    #
    # Linearly Regression
    slope, intercept, rvalue, pvalue, stderr = stats.linregress(x,y)
    #
    # Regression Line
    regressline = slope*x+intercept
    #
    return regressline, rvalue, pvalue, stderr


# Linear Regression
line, r_value, p_value, std_err = linearregress(demeaned_floatMKE,demeaned_chlyavg30meter)

# Regression Plot
pp = PdfPages('floatmke-chly-daily-regression.pdf')
fig = plt.figure()
ax = fig.add_subplot(111)
line1,=ax.plot(demeaned_floatMKE,demeaned_chlyavg30meter,'o', )
ax.set_xlabel("MKE")
ax.set_ylabel("30m Avg Chl")
line2,=ax.plot(demeaned_floatMKE,line,'k-')

fig.suptitle('R value: ' + str(r_value) + ',  P value: ' + str(p_value))
plt.savefig(pp, format = "pdf")
pp.close()
os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/floatmke-chly-daily-regression.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')
fghfhjf










































gridWplt = 2
if gridWplt == 1:    
    #
    # ===========================================================
    # Contour Plotting Grided Density W Data
    # ===========================================================
    #
    # Changing Directory to plotting output
    os.chdir('/home/ardavies/satdata/OSCAR/pdfoutput')
    #
    # Plot Set-up
    pp = PdfPages('GriddedW_Density.pdf')
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    #
    # Make Blue-Red Color Scheme
    from matplotlib.colors import LinearSegmentedColormap
    cdict3 = {'red':  ((0.0, 0.0, 0.0),
                       (0.25,0.0, 0.0),
                       (0.5, 0.8, 1.0),
                       (0.75,1.0, 1.0),
                       (1.0, 0.4, 1.0)),

             'green': ((0.0, 0.0, 0.0),
                       (0.25,0.0, 0.0),
                       (0.5, 0.9, 0.9),
                       (0.75,0.0, 0.0),
                       (1.0, 0.0, 0.0)),

             'blue':  ((0.0, 0.0, 0.4),
                       (0.25,1.0, 1.0),
                       (0.5, 1.0, 0.8),
                       (0.75,0.0, 0.0),
                       (1.0, 0.0, 0.0))
            }
    #
    # Make a modified version of cdict3 with some transparency
    # in the middle of the range.
    cdict4 = cdict3.copy()
    cdict4['alpha'] = ((0.0, 1.0, 1.0),
                    #   (0.25,1.0, 1.0),
                       (0.5, 0.3, 0.3),
                    #   (0.75,1.0, 1.0),
                       (1.0, 1.0, 1.0))
    plt.register_cmap(name='BlueRedAlpha', data=cdict4)
    #
    # Text/Front Set up
    from matplotlib import rcParams
    rcParams['axes.labelsize'] = 14
    rcParams['xtick.labelsize'] = 12
    rcParams['ytick.labelsize'] = 12
    rcParams['legend.fontsize'] = 10
    rcParams['font.family'] = 'serif'
    rcParams['font.serif'] = ['Computer Modern Roman']
    rcParams['text.usetex'] = True
    #
    # Contour Plotting & Color Bar
    from pylab import *
    cs = plt.contourf(GridWDates[:,1,:],GridWDepths[:,1,:],GridWData[:,1,:]*10**(4), levels = np.linspace(-6,6,201))
    plt.set_cmap('BlueRedAlpha')
    cbar = plt.colorbar(cs,spacing='proportional')
    cbar.set_label(r"\textbf{Isopycnal Vertical Velocity $\times 10**4$ (ms$^{-1}$)}")
    cbar.set_ticks([-5,-2.5, 0,2.5, 5])
    cbar.set_ticklabels([-5,-2.5, 0,2.5, 5])
    plt.plot(IsoBioWDatesDepthAvgAvg, IsoBioWDepthsDepthAvgAvg, 'k',label=r'Averaged Depth of 3 Profile Averaged Chlorophyll and Sinking Vertical Velocties',linewidth = 2)
    l = legend(loc = 2)
    ax.set_xlabel(r"\textbf{Day of Year}")
    ax.set_xlim([70,110])
    ax.set_ylim([-900,0])
    ax.set_ylabel(r"\textbf{Depth (m)}")
    #
    # Save and spc to storm
    plt.savefig(pp, format = "pdf")
    pp.close()
    os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/GriddedW_Density.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/Jan01_Jun04_2013/verticalvelocities/')
#
# ===================================================================================
# PLOTTING MKE LINES
# ===================================================================================
#
# Changing Directory to plotting output
os.chdir('/home/ardavies/satdata/OSCAR/pdfoutput')

#
# Float MKE Line Plots

pp = PdfPages('floatmke-chly-daily_secondmission.pdf')
fig = plt.figure()
from matplotlib import rcParams
rcParams['axes.labelsize'] = 14
rcParams['xtick.labelsize'] = 12
rcParams['ytick.labelsize'] = 12
rcParams['legend.fontsize'] = 10
rcParams['font.family'] = 'serif'
rcParams['font.serif'] = ['Computer Modern Roman']
rcParams['text.usetex'] = True

ax = fig.add_subplot(1, 1, 1)
from pylab import *

ax.plot(floatdate,chlyavg30meter,'b',label=r'30 meter Avg Chl Concentration',linewidth = 3)
ax.plot(floatdate,floatMKE,'k',label=r'Using Daily Averaged OSCAR Data',linewidth = 3)
l = legend(loc = 2)
ax.set_yscale('log')
ax.set_xlabel(r"\textbf{Day of Year}")
ax.set_ylabel(r"\textbf{Mean Kinetic Energy (cm$^2$s$^{-2}$)}")
plt.savefig(pp, format = "pdf")
pp.close()
os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/floatmke-chly-daily_secondmission.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')


lineplot = 2
if lineplot == 1:
    #
    # Plot Distances between float origin, next float location, origin, and so on..
    pp = PdfPages('allthree_secondmission.pdf')
    fig = plt.figure()
    from matplotlib import rcParams
    rcParams['axes.labelsize'] = 18
    rcParams['xtick.labelsize'] = 18
    rcParams['ytick.labelsize'] = 18
    rcParams['legend.fontsize'] = 16
    rcParams['font.family'] = 'serif'
    rcParams['font.serif'] = ['Computer Modern Roman']
    rcParams['text.usetex'] = True
    #
    ax = fig.add_subplot(1, 1, 1)
    ax.set_yscale('linear')
    from pylab import *
    ax.plot(floatdate[0:arraylen-1],origtofloat,'k',label=r'Current Float to Next Float Distance',linewidth = 3)
    ax.plot(floatdate[0:arraylen-1],tracertofloat,'b',label=r'Next Float to Tracer Distance',linewidth = 3)
    ax.plot(floatdate[0:arraylen-1],origtotracer,'0.75',label=r'Current Float to Tracer Distance',linewidth = 3)
    l = legend(loc = 2)
    ax.set_xlabel(r"\textbf{Day of Year}")
    ax.set_xlim([275,350])
    #ax.set_ylim([0,110])
    ax.set_ylabel(r"\textbf{Distance (km)}")
    plt.savefig(pp, format = "pdf")
    pp.close()
    os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/allthree_secondmission.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')
    #
    # Plot Lagrangian Ratio
    pp = PdfPages('ratio_secondmission.pdf')
    fig = plt.figure()
    from matplotlib import rcParams
    rcParams['axes.labelsize'] = 14
    rcParams['xtick.labelsize'] = 12
    rcParams['ytick.labelsize'] = 12
    rcParams['legend.fontsize'] = 10
    rcParams['font.family'] = 'serif'
    rcParams['font.serif'] = ['Computer Modern Roman']
    rcParams['text.usetex'] = True
    #    
    ax = fig.add_subplot(1, 1, 1)
    from pylab import *
    from matplotlib.font_manager import FontProperties
    fontP = FontProperties()
    fontP.set_size('x-small')
    ax.set_yscale('linear')
    ax.plot(floatdate[0:arraylen-1],tracerfloatratio,'k',linewidth = 3)
    ax.set_xlabel(r"\textbf{Day of Year}")
    ax.set_xlim([275,350])
    #ax.set_ylim([0,2])
    ax.set_ylabel(r"\textbf{Lagrangian Ratio}")
    plt.savefig(pp, format = "pdf")
    pp.close()
    os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/ratio_secondmission.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')
    #
    # Float MKE Line Plots
    pp = PdfPages('floatmke-daily_secondmission.pdf')
    fig = plt.figure()
    from matplotlib import rcParams
    rcParams['axes.labelsize'] = 14
    rcParams['xtick.labelsize'] = 12
    rcParams['ytick.labelsize'] = 12
    rcParams['legend.fontsize'] = 10
    rcParams['font.family'] = 'serif'
    rcParams['font.serif'] = ['Computer Modern Roman']
    rcParams['text.usetex'] = True

    ax = fig.add_subplot(1, 1, 1)
    from pylab import *

    ax.plot(floatdate,floatMKE5day,'0.75',label=r'From Original, 5 Day Averaged Data',linewidth = 3)
    ax.plot(floatdate,floatMKE,'k',label=r'Using Daily Averaged OSCAR Data',linewidth = 3)
    l = legend(loc = 2)
    ax.set_yscale('log')
    ax.set_xlabel(r"\textbf{Day of Year}")
    ax.set_ylabel(r"\textbf{Mean Kinetic Energy (cm$^2$s$^{-2}$)}")
    plt.savefig(pp, format = "pdf")
    pp.close()
    os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/floatmke-daily_secondmission.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')




    #
    # Float MKE Line Plots-Cropped
    pp = PdfPages('floatmke-daily-compare-crop_secondmission.pdf')
    fig = plt.figure()
    from matplotlib import rcParams
    rcParams['axes.labelsize'] = 14
    rcParams['xtick.labelsize'] = 12
    rcParams['ytick.labelsize'] = 12
    rcParams['legend.fontsize'] = 10
    rcParams['font.family'] = 'serif'
    rcParams['font.serif'] = ['Computer Modern Roman']
    rcParams['text.usetex'] = True
    #   
    ax = fig.add_subplot(1, 1, 1)
    ax.set_yscale('log')
    ax.plot(floatdate,floatMKE5day,'0.75',label=r'From Original, 5 Day Averaged Data',linewidth = 3)
    ax.plot(floatdate,floatMKE,'k',label=r'Using Daily Averaged OSCAR Data',linewidth = 3)
    l = legend(loc = 2)
    ax.set_xlabel(r"\textbf{Day of Year}")
    ax.set_ylabel(r"\textbf{Mean Kinetic Energy (cm$^2$s$^{-2}$)}")
    ax.set_xlim([275,350])
    #ax.set_ylim([30,2000])
    plt.savefig(pp, format = "pdf")
    pp.close()
    os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/floatmke-daily-compare-crop_secondmission.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')
    # #
    # # Float MKE Line Plots-Cropped
    # pp = PdfPages('floatmke-daily-crop_secondmission.pdf')
    # fig = plt.figure()

    # from matplotlib import rcParams
    # rcParams['axes.labelsize'] = 14
    # rcParams['xtick.labelsize'] = 12
    # rcParams['ytick.labelsize'] = 12
    # rcParams['legend.fontsize'] = 10
    # rcParams['font.family'] = 'serif'
    # rcParams['font.serif'] = ['Computer Modern Roman']
    # rcParams['text.usetex'] = True

    # ax = fig.add_subplot(1, 1, 1)
    # ax.set_yscale('log')
    # ax.plot(floatdate,floatMKE,'k',label=r'Using Daily Averaged OSCAR Data',linewidth = 3)
    # ax.set_xlabel(r"\textbf{Day of Year}")
    # ax.set_ylabel(r"\textbf{Mean Kinetic Energy (cm$^2$s$^{-2}$)}")
    # ax.set_xlim([275,350])
    # #ax.set_ylim([30,2000])
    # plt.savefig(pp, format = "pdf")
    # pp.close()
    # os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/floatmke-daily-crop_secondmission.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')
#
# MKE Zoomed in Plotting
pp = PdfPages('currents_secondmission.pdf')
figtit1 = 'OSCAR 5 Day Avg '
figtit2 = ' Centered About DOY '
from matplotlib.colors import LogNorm

for l in range(0,arraylen-7):
    #
    # Find the correct OSCAR data plot
    k = 0
    for kk in range(0,len(oscardates)-5):
        if (oscardates[k] == floatdate[l]):
            break
        else:
            k = k + 1
    #
    usefloatlon= floatlon[l]
    usefloatlat = floatlat[l]
    usefloatlon_new = floatlon[l+1]
    usefloatlat_new = floatlat[l+1]
    usetracerlat_U_V = tracer_newlat_U_V[l]
    usetracerlon_U_V = tracer_newlon_U_V[l]
    #
    import matplotlib.path as mpath
    import matplotlib.patches as mpatches
    import matplotlib as mpl
    from matplotlib.collections import PatchCollection
    Path = mpath.Path
    fig=plt.figure(figsize=(8,4.5))
    m = Basemap(llcrnrlon=309,llcrnrlat=-53,urcrnrlon=337,urcrnrlat=-43,projection='mill',resolution='i')
    ax = fig.add_axes([0.05,0.05,0.9,0.85])
    m.drawcoastlines(linewidth=0.25)
    m.drawcountries(linewidth=0.25)
    m.fillcontinents(color='gray',lake_color='white')
    m.drawmapboundary(fill_color='white')
    m.drawparallels(np.arange(-90.,90,5.),labels=[1,0,0,0])
    m.drawmeridians(np.arange(180.,360.,5. ),labels=[0,0,0,1])
    x, y = m(lon, lat)
    floatx,floaty =m(usefloatlon, usefloatlat)
    #floatx_new,floaty_new =m(usefloatlon_new, usefloatlat_new)
    #tracerx_U_V,tracery_U_V = m(usetracerlon_U_V, usetracerlat_U_V)
    #cs = m.scatter(floatx_new,floaty_new, s = 10, marker=(5,0), c='b')    
    cs1 = m.scatter(floatx,floaty, s = 20, marker=(5,0), c='r')
    #cs2 = m.scatter(tracerx_U_V,tracery_U_V, s = 10, marker=(5,0), c='r')
    Q = m.quiver(x,y,udaily[k,:,:],vdaily[k,:,:])
    qk = plt.quiverkey(Q, .02, .95, 1, '1 m/s', labelpos='W',color = 'r')
    plt.title('OSCAR Currents, Float and Tracer on DOY: '  + str(int(oscardates[k])))
    plt.savefig(pp, format = "pdf")

pp.close()
os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/currents_secondmission.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')




pp = PdfPages('currents-daily_secondmissoin.pdf')
figtit1 = 'OSCAR 5 Day Avg Daily Plot on DOY: '
from matplotlib.colors import LogNorm

for ii in range(0,len(oscardates)-5):    #
    # 
    import matplotlib.path as mpath
    import matplotlib.patches as mpatches
    import matplotlib as mpl
    from matplotlib.collections import PatchCollection
    Path = mpath.Path
    fig=plt.figure(figsize=(8,4.5))
    m = Basemap(llcrnrlon=309,llcrnrlat=-53,urcrnrlon=337,urcrnrlat=-43,projection='mill',resolution='i')
    ax = fig.add_axes([0.05,0.05,0.9,0.85])
    m.drawcoastlines(linewidth=0.25)
    m.drawcountries(linewidth=0.25)
    m.fillcontinents(color='gray',lake_color='white')
    m.drawmapboundary(fill_color='white')
    m.drawparallels(np.arange(-90.,90,5.),labels=[1,0,0,0])
    m.drawmeridians(np.arange(180.,360.,5. ),labels=[0,0,0,1])
    x, y = m(lon, lat)
    Q = m.quiver(x,y,udaily[ii,:,:],vdaily[ii,:,:])
    qk = plt.quiverkey(Q,.02, .95, 1, '1 m/s', labelpos='W',color = 'r')
    for k in range(0,numdates-1):
        if (centerdates[k] == oscardates[ii]):
                Q = m.quiver(x,y,u[k,:,:],v[k,:,:], color = 'b')
    plt.title(figtit1 + str(int(oscardates[ii])))
    plt.savefig(pp, format = "pdf")
pp.close()

os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/currents-daily_secondmissoin.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')



#
# ===================================================================================
# Plotting FLoat Locatoins, Trajectories, and MKE and Each Space.
# ===================================================================================
#
toplot = 2
if toplot == 1:
    pp = PdfPages('floatlocations.pdf')
    from matplotlib.colors import LogNorm
    from mpl_toolkits.basemap import Basemap
    import matplotlib.pyplot as plt
    import numpy as np
    fig = plt.figure() 
    # ax = fig.add_axes()
    ax = fig.add_axes([0.07,0.06,0.83,0.97])

    #f, ax = plt.subplots()
    #
    # Setting Fonts
    from matplotlib import rcParams
    rcParams['axes.labelsize'] = 16
    rcParams['xtick.labelsize'] = 16
    rcParams['ytick.labelsize'] = 16
    rcParams['legend.fontsize'] = 16
    rcParams['font.family'] = 'serif'
    rcParams['font.serif'] = ['Computer Modern Roman']
    rcParams['text.usetex'] = True
    #
    # Setting the Larger Map
    m = Basemap(llcrnrlon=300.8,llcrnrlat=-59.5,urcrnrlon=305.9,urcrnrlat=-55.5,projection='cyl',resolution='i')
    m.drawcoastlines(linewidth=0.25)
    m.drawcountries(linewidth=0.25)
    m.fillcontinents(color='gray',lake_color='white')
    m.drawmapboundary(fill_color='white')
    m.drawparallels(np.arange(-90.,90,2.),labels=[1,0,0,0],linewidth=0.0,fontsize=16)
    m.drawmeridians(np.arange(180.,360.,2. ),labels=[0,0,0,1],linewidth=0.0,fontsize=16)
    #
    # Plotting Tracer and Float locations over the correct dates
    x,y = m(lon,lat)
    for l in range(0,21):
        ll = l + 30
        #
        # Find the correct OSCAR data plot
        k = 0
        for kk in range(0,len(oscardates)-5):
            if (oscardates[k] == floatdate[ll]):
                break
            else:
                k = k + 1
        #
        usefloatlon= floatlon[ll]
        usefloatlat = floatlat[ll]
        usefloatlon_new = floatlon[ll+1]
        usefloatlat_new = floatlat[ll+1]
        usetracerlat_U_V = tracer_newlat_U_V[ll]
        usetracerlon_U_V = tracer_newlon_U_V[ll]
        #
        floatx,floaty =m(usefloatlon, usefloatlat)
        #
        # Plotting Float Locations And Tracers
        tracerx_U_V,tracery_U_V =m(usetracerlon_U_V, usetracerlat_U_V)
        cs3 = m.plot([floatx,tracerx_U_V],[floaty,tracery_U_V],linewidth = 0.5, color = '#696969')
        cs1 = m.scatter(tracerx_U_V,tracery_U_V, s = 5, c='#696969',edgecolors='none')
    #
    # Plotting MKE Color for each float date and the DOY
    floatx2 = np.zeros([21])
    floaty2 = np.zeros([21])
    floatx22 = np.zeros([11])
    floaty22 = np.zeros([11])
    floatmke2 = np.zeros([21])
    floatdates2 = [None]*11
    offsetx = np.zeros([11])
    offsety = np.zeros([11])
    offsetx = ([  0,    0, -.1, -.13, -.15, -.13,  -.13, -.15,  .15, -.15,  0])
    offsety = ([.14,  .13,  .1, 0.03,    0,    0,  -0.01,    0,    0,   .1, .15])
    lll = 0
    for l in range(0,21):
        ll = l + 30
        floatx2[l],floaty2[l] = m(floatlon[ll], floatlat[ll])
        floatmke2[l] = floatMKE[ll]
        if l%2==0:
            floatx22[lll] = floatx2[l] + offsetx[lll]
            floaty22[lll] = floaty2[l] + offsety[lll]
            if lll == 0:
                floatdates2[lll] = "DOY: " + str(int(floatdate[ll]))
            else:
                floatdates2[lll] = str(int(floatdate[ll]))
            lll = lll + 1
    for label, xpt, ypt in zip(floatdates2, floatx22, floaty22):
        plt.text(xpt, ypt, label,ha='center',va='center',fontsize=12, family='Computer Modern Roman')
   
    import math as ma
    from mpl_toolkits.axes_grid1 import make_axes_locatable

    xponent1 = ma.log(30)/ma.log(10)
    xponent2 = ma.log(2000)/ma.log(10)
    # cs2 = m.scatter(floatx2,floaty2,s = 35,c=floatmke2, norm=matplotlib.colors.LogNorm())#,levels= np.logspace(xponent1,xponent2,100))
    # cbar = plt.colorbar(cs2,spacing='proportional', norm = LogNorm())

    cs2 = m.scatter(floatx2,floaty2,s = 40,c=floatmke2, norm=LogNorm(vmin=floatmke2.min(), vmax=floatmke2.max()))# , norm=matplotlib.colors.LogNorm())#,levels= np.logspace(xponent1,xponent2,100))
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.05)
    cbar = plt.colorbar(cs2,cax=cax, norm=LogNorm(vmin=floatmke2.min(), vmax=floatmke2.max()))
    minorticks = cs2.norm(np.array([60, 70, 80, 90, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000]))
    cbar.ax.yaxis.set_ticks(minorticks, minor=True)

    print np.log(floatmke2)

    cbar.set_label(r"\textbf{Mean Kinetic Energy (cm$^2$s$^{-2}$)}")
    cbar.set_ticks([100, 1000])
    cbar.set_ticklabels([r'10$^{2}$', r'10$^{3}$'])
    #
    # Base Map set-up for the insert map
    m2 = Basemap(llcrnrlon=292.2,llcrnrlat=-64.2,urcrnrlon=308.5,urcrnrlat=-50.9,projection='cyl',resolution='i')
    ax2 = fig.add_axes([0.061,0.545,0.4,0.4])
    m2.drawcoastlines(linewidth=0.25)
    m2.drawcountries(linewidth=0.25)
    m2.fillcontinents(color='gray',lake_color='white')
    m2.drawmapboundary(fill_color='white')
    #
    # Plotting all the float locations in gray
    floatx3,floaty3 =m2(floatlon, floatlat)
    cs7 = m2.plot(floatx3,floaty3,linewidth = 0.2, color = '#696969')
    cs8 = m2.scatter(floatx3,floaty3, s = 3, c='#696969',edgecolors='none')
    #
    # Overplotting the ones used in larger figure in red
    floatx4 = np.zeros([21])
    floaty4 = np.zeros([21])
    lll = 0
    for l in range(0,21):
        ll = l + 30
        floatx4[l],floaty4[l] = m2(floatlon[ll], floatlat[ll])
        lll = lll + 1
    cs9 = m2.plot(floatx4,floaty4,linewidth = 0.2, color = 'b')
    cs10 = m2.scatter(floatx4,floaty4, s = 6, c='b',edgecolors='b')
    #
    # Plotting Labels
    places = [r'\textbf{Argentina}', r'\textbf{Falkland Islands}', r'\textbf{Antarctica}']
    placeslons = ([295.6, 305.0, 295.5])
    placeslats = ([-53.6, -52.465, -62.3])
    placesx, placesy = m2(placeslons, placeslats)
    for place, xpt2, ypt2 in zip(places, placesx, placesy):
        plt.text(xpt2, ypt2, place,ha='center',va='center',fontsize=11, family = 'Computer Modern Roman')
    #
    # Saving
    plt.savefig(pp, format = "pdf")
    pp.close()
    os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/floatlocations.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/grlpaperplots/')
#
# ===================================================================================
# Plotting 
# ===================================================================================
#
sideplot = 2
if sideplot == 1:
    from matplotlib import rc
    from matplotlib.numerix import arange, cos, pi
    rc('text', usetex=True)
    from pylab import *
    numbiolinesavg = len(IsoBioWDataAvg)
    pp = PdfPages('figure3a.pdf')
    fig = plt.figure()
    from matplotlib.font_manager import FontProperties
    legendfont = FontProperties()
    legendfont.set_name('Computer Modern Roman')
    legendfont.set_size('x-small')
    rcParams['axes.labelsize'] = 18
    rcParams['xtick.labelsize'] = 18
    rcParams['ytick.labelsize'] = 18
    rcParams['legend.fontsize'] = 16

    from matplotlib import rcParams
    rcParams['font.family'] = 'serif'
    rcParams['font.serif'] = ['Computer Modern Roman']
    rcParams['text.usetex'] = True

    #ax = fig.add_subplot(1, 1, 1)
    ax = fig.add_axes([0.15,0.1,0.81,0.68])

    #
    # Which Plotting Information?
    # axhline(0, color='k')
    for j in range(0,numbiolinesavg):
        if j == 0:
            plt.plot(IsoBioWDatesAvg[j], IsoBioWDataAvg[j]*10**4, '0.75',label=r'Chl Velocities')
        else:
            plt.plot(IsoBioWDatesAvg[j], IsoBioWDataAvg[j]*10**4, '0.75')
    plt.plot(IsoBioWDatesDepthAvgAvg, NearestWGriddedAvgAvg*10**4, 'b',label=r'Environmental Velocity',linewidth = 2)
    plt.plot(IsoBioWDatesDepthAvgAvg, SinkWDataDepthAvgAvg*10**4, 'r',label=r'Depth Avg Sinking Velocity',linewidth = 2)
    plt.plot(IsoBioWDatesDepthAvgAvg, IsoBioWDataDepthAvgAvg*10**4, 'k',label=r'Depth Avg Chl Velocity',linewidth = 2)
    l = legend(loc = 3)
    #ax.set_xlabel("Days since 01/01/2013")
    ax.set_ylabel(r'\textbf{Veloctiy $\times 10^4$ (ms$^{-1}$)}' )
    ax.get_yaxis().set_label_coords(-0.1,0.5)
    ax.set_xlim([70,110])
    ax.set_ylim([-13,4])
    plt.setp(ax.get_xticklabels(), visible=False)
    plt.savefig(pp, format = "pdf")
    pp.close()
    os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/figure3a.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/grlpaperplots/')

    #ax.set(adjustable='box-forced', aspect=1)
    #
    #
sideplot = 2
if sideplot == 1:
    from matplotlib import rc
    from matplotlib.numerix import arange, cos, pi
    rc('text', usetex=True)
    from pylab import *
    numbiolinesavg = len(IsoBioWDataAvg)
    pp = PdfPages('figure3b.pdf')
    fig = plt.figure()
    from matplotlib.font_manager import FontProperties
    legendfont = FontProperties()
    legendfont.set_name('Computer Modern Roman')
    legendfont.set_size('x-small')
    rcParams['axes.labelsize'] = 18
    rcParams['xtick.labelsize'] = 18
    rcParams['ytick.labelsize'] = 18
    rcParams['legend.fontsize'] = 16

    from matplotlib import rcParams
    rcParams['font.family'] = 'serif'
    rcParams['font.serif'] = ['Computer Modern Roman']
    rcParams['text.usetex'] = True

    ax2 = fig.add_axes([0.15,0.1,0.81,0.68])
    #
    # Which Plotting Information?
    axhline(0, color='k')
    for j in range(0,numbiolinesavg):
        if j == 0:
            plt.plot(IsoBioWDatesAvg[j], IsoBioWDepthsAvg[j], '0.75',label=r'Depth of Chl Velocities')
        else:
            plt.plot(IsoBioWDatesAvg[j], IsoBioWDepthsAvg[j], '0.75')
    plt.plot(IsoBioWDatesDepthAvgAvg, IsoBioWDepthsDepthAvgAvg, 'k',label=r'Depth of Averaged Velocties',linewidth = 2)
    #
    #
    l2 = legend(loc = 3)
    #ax2.set_xlabel("Days since 01/01/2013")
    ax2.set_ylabel(r'\textbf{Depth (m)}')
    ax2.get_yaxis().set_label_coords(-0.1,0.5)
    ax2.set_xlim([70,110])
    ax2.set_ylim([-700, -90])
    ax2.set_yticks([0, -200, -400, -600])
    plt.setp(ax2.get_xticklabels(), visible=False)
    plt.savefig(pp, format = "pdf")
    pp.close()
    os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/figure3b.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/grlpaperplots/')
    #
sideplot = 2
if sideplot == 1:
    from matplotlib import rc
    from matplotlib.numerix import arange, cos, pi
    rc('text', usetex=True)
    from pylab import *
    numbiolinesavg = len(IsoBioWDataAvg)
    pp = PdfPages('figure3c.pdf')
    fig = plt.figure()
    from matplotlib.font_manager import FontProperties
    legendfont = FontProperties()
    legendfont.set_name('Computer Modern Roman')
    legendfont.set_size('x-small')
    rcParams['axes.labelsize'] = 18
    rcParams['xtick.labelsize'] = 18
    rcParams['ytick.labelsize'] = 18
    rcParams['legend.fontsize'] = 16

    from matplotlib import rcParams
    rcParams['font.family'] = 'serif'
    rcParams['font.serif'] = ['Computer Modern Roman']
    rcParams['text.usetex'] = True

    ax3 = fig.add_axes([0.15,0.1,0.81,0.68])
    #
    #
    #ax3 = fig.add_subplot(1, 1, 1)
    ax3.set_yscale('log')
    ax3.plot(floatdate,floatMKE,'k',label=r'Float MKE',linewidth = 3)
    l3 = legend(loc = 3)
    #l3 = legend(loc = 3,prop = legendfont)
    ax3.set_xlabel(r'\textbf{Days since 01/01/2013}')
    ax3.set_ylabel(r'\textbf{MKE (cm$^2$ s$^{-2}$)}')
    ax3.get_yaxis().set_label_coords(-0.1,0.5)
    ax3.set_xlim([70,110])
    ax3.set_ylim([40,1500])
    plt.savefig(pp, format = "pdf")
    pp.close()
    os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/figure3c.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/grlpaperplots/')





contplt = 2
if contplt == 1:
    from matplotlib import rc
    from matplotlib.numerix import arange, cos, pi
    rc('text', usetex=True)
    from pylab import *
    numbiolinesavg = len(IsoBioWDataAvg)
    pp = PdfPages('figure2a.pdf')
    fig = plt.figure()
    from matplotlib.font_manager import FontProperties
    legendfont = FontProperties()
    legendfont.set_name('Computer Modern Roman')
    legendfont.set_size('x-small')
    rcParams['axes.labelsize'] = 18
    rcParams['xtick.labelsize'] = 18
    rcParams['ytick.labelsize'] = 18
    rcParams['legend.fontsize'] = 16

    from matplotlib import rcParams
    rcParams['font.family'] = 'serif'
    rcParams['font.serif'] = ['Computer Modern Roman']
    rcParams['text.usetex'] = True
    #
    # Plotting correct data
    contourdat = np.zeros([interplength,arraylen])
    for ci in range(0,interplength):
        for cj in range(0,arraylen):
            #cjj = cj + 1
            if GridData[ci,4,cj] < 10**-2.5:
                contourdat[ci,cj] = 10**-2.5
            else:
                contourdat[ci,cj] = GridData[ci,4,cj]
    #
    # Plot Set-up
    import math as ma
    from mpl_toolkits.axes_grid1 import make_axes_locatable
    fig = plt.figure()
    ax = fig.add_axes([0.15,0.1,0.70,0.85])
    #ax = fig.add_subplot(1, 1, 1)

    contlevels= np.logspace(-2.5,np.log(contourdat.max())/np.log(10),250)
    conticks = [0.01,0.1,1.00]
    cs = plt.contourf(floatdate,Ygrid,contourdat, levels= contlevels, norm=LogNorm())
    cs = plt.contourf(floatdate,Ygrid,contourdat, levels= contlevels, norm=LogNorm())
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.05)

    cbar = plt.colorbar(cs,cax=cax, norm=LogNorm())
    minorticks = cs.norm(np.array([0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9, 1.00, 2.00, 3.00, 4.00, 5.00, 6,00]))
    cbar.ax.yaxis.set_ticks(minorticks, minor=True)

    cbar.set_label(r"\textbf{Chlorophyll Concentration ($\mu$gl$^{-1}$)}")
    cbar.set_ticks([0.01, 0.1,1.00])
    cbar.set_ticklabels([r'10$^{-2}$',r'10$^{-1}$', r'10$^{0}$'])
    axvline(70, linewidth=0.5, color='0.75')
    axvline(110, linewidth=0.5, color='0.75')
    # Plotting
    #figtit2 = ' and Chlorophyll Quivers'
    #savetit2 = "_IsoChlyW_Quivers"
    #from pylab import *
    #Q2 = plt.quiver(IsoBioWDates[1::10], IsoBioWDepths[1::10], IsoBioFake[1::10], IsoBioWData[1::10], color='k')
    #plt.quiverkey(Q2, 0.22, 0.96, 0.0005, r'$ 5 \times 10^{-4} m/s$', labelpos='W')
    #ax = plt.gca()
    ax.set_yticks([0,-200, -400, -600, -800])
    #ax.yaxis.set_yticks([-200,-400, -600, -800])
    ax.set_yticklabels([r'0', r'200',r'400',r'600',r'800'])
    #ax.set_yticks([-100,-200,-300,-400,-500,-600,-700,-800,-900], minor=true)


    ax.set_xlabel(r'\textbf{Days since 01/01/2013}')
    ax.set_ylabel(r'\textbf{Depth (m)}')
    ax.set_ylim([-900,0])
    # plt.title(figtit)
    plt.savefig(pp, format = "pdf")
    pp.close()
    os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/figure2a.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/grlpaperplots/')


contplt = 2
if contplt == 1:
    from matplotlib import rc
    from matplotlib.numerix import arange, cos, pi
    rc('text', usetex=True)
    from pylab import *
    numbiolinesavg = len(IsoBioWDataAvg)
    pp = PdfPages('figure2b.pdf')
    fig = plt.figure()
    from matplotlib.font_manager import FontProperties
    legendfont = FontProperties()
    legendfont.set_name('Computer Modern Roman')
    legendfont.set_size('x-small')
    rcParams['axes.labelsize'] = 18
    rcParams['xtick.labelsize'] = 18
    rcParams['ytick.labelsize'] = 18
    rcParams['legend.fontsize'] = 16

    from matplotlib import rcParams
    rcParams['font.family'] = 'serif'
    rcParams['font.serif'] = ['Computer Modern Roman']
    rcParams['text.usetex'] = True
    #
    # Plotting correct data
    contourdat = np.zeros([interplength,arraylen])
    for ci in range(0,interplength):
        for cj in range(0,arraylen):
            #cjj = cj + 1
            contourdat[ci,cj] = GridData[ci,1,cj]
    print contourdat.max()
    print contourdat.min()
    #
    # Plot Set-up
    import math as ma
    from mpl_toolkits.axes_grid1 import make_axes_locatable
    fig = plt.figure()
    ax = fig.add_axes([0.15,0.1,0.70,0.85])

    contlevels= np.linspace(contourdat.min(),contourdat.max(),250)
    conticks = [1026.8, 1027.0, 1027.2, 1027.4, 1027.6, 1027.8]
    cs = plt.contourf(floatdate,Ygrid,contourdat, levels= contlevels)
    cs = plt.contourf(floatdate,Ygrid,contourdat, levels= contlevels)
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.05)
    cbar = plt.colorbar(cs,cax=cax, spacing='proportional')
    minorticks = cs.norm(np.array([1026.8, 1026.9, 1027.0, 1027.1, 1027.2, 1027.3, 1027.4, 1027.5, 1027.6, 1027.7, 1027.8]))
    cbar.ax.yaxis.set_ticks(minorticks, minor=True)

    cbar.set_label(r"\textbf{Density (kgm$^{-3}$)}")
    cbar.set_ticks([1026.8, 1027.0, 1027.2, 1027.4, 1027.6, 1027.8])
    #cbar.set_ticklabels([r'10$^{-1}$', r'10$^{0}$'])
    axvline(70, linewidth=0.5, color='0.75')
    axvline(110, linewidth=0.5, color='0.75')
    # Plotting
    #figtit2 = ' and Chlorophyll Quivers'
    #savetit2 = "_IsoChlyW_Quivers"
    #from pylab import *
    #Q2 = plt.quiver(IsoBioWDates[1::10], IsoBioWDepths[1::10], IsoBioFake[1::10], IsoBioWData[1::10], color='k')
    #plt.quiverkey(Q2, 0.22, 0.96, 0.0005, r'$ 5 \times 10^{-4} m/s$', labelpos='W')
    #ax = plt.gca()
    ax.set_yticks([0,-200, -400, -600, -800])
    #ax.yaxis.set_yticks([-200,-400, -600, -800])
    ax.set_yticklabels([r'0', r'200',r'400',r'600',r'800'])

    ax.set_xlabel(r'\textbf{Days since 01/01/2013}')
    ax.set_ylabel(r'\textbf{Depth (m)}')
    ax.set_ylim([-900,0])
    # plt.title(figtit)
    plt.savefig(pp, format = "pdf")
    pp.close()
    os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/figure2b.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/grlpaperplots/')








# # Let's add some earthquakes (fake here) :
 
# lon = np.random.random_integers(11,79,1000)/10.
# lat = np.random.random_integers(491,519,1000)/10.
# depth = np.random.random_integers(0,300,1000)/10.
# magnitude = np.random.random_integers(0,100,1000)/10.
 
# # I'm masking the earthquakes present in most of the regions (illustration of masks usage) :
# import numpy.ma as ma
# Mlon = ma.masked_outside(lon, 5.6, 7.5)
# Mlat = ma.masked_outside(lat,49.6,50.6)
# lat = ma.array(lat,mask=Mlon.mask+Mlat.mask).compressed()
# lon = ma.array(lon,mask=Mlon.mask+Mlat.mask).compressed()
# depth =ma.array(depth,mask=Mlon.mask+Mlat.mask).compressed()
# magnitude = ma.array(magnitude,mask=Mlon.mask+Mlat.mask).compressed()
 

 
# # adding a zoom plot :
# from mpl_toolkits.axes_grid.inset_locator import zoomed_inset_axes
# from mpl_toolkits.axes_grid.inset_locator import mark_inset
# from mpl_toolkits.axes_grid.anchored_artists import AnchoredSizeBar
 
# axins = zoomed_inset_axes(ax, 2, loc=2) # zoom = 4
# m.scatter(x,y,s=10*magnitude,c=depth)
# x,y = m(5.5,50.0)
# x2,y2 = m(6.5,50.5)
# axins.set_xlim(x,x2)
# axins.set_ylim(y,y2)
# plt.xticks(visible=False)
# plt.yticks(visible=False)
# mark_inset(ax, axins, loc1=1, loc2=3, fc="none", ec="0.5")
# asb = AnchoredSizeBar(axins.transData,
#  10000.,
#  "10 km",
#  loc=4,
#  pad=0.1, borderpad=0.25, sep=5,
#  frameon=False)
# axins.add_artist(asb)
 
# plt.show()





















# #
# # Zoomed-in Daily Currents
# pp = PdfPages('currents-daily-zoomed.pdf')
# figtit1 = 'Daily Averaged OSCAR Data on DOY: '
# from matplotlib.colors import LogNorm
# for ii in range(0,len(oscardates)-5):    #
#     # 
#     import matplotlib.path as mpath
#     import matplotlib.patches as mpatches
#     import matplotlib as mpl
#     from matplotlib.collections import PatchCollection
#     Path = mpath.Path
#     fig=plt.figure(figsize=(8,4.5))
#     m = Basemap(llcrnrlon=297,llcrnrlat=-61,urcrnrlon=307,urcrnrlat=-54,projection='cyl',resolution='i')
#     ax = fig.add_axes([0.05,0.05,0.9,0.85])
#     m.drawcoastlines(linewidth=0.25)
#     m.drawcountries(linewidth=0.25)
#     m.fillcontinents(color='gray',lake_color='white')
#     m.drawmapboundary(fill_color='white')
#     m.drawparallels(np.arange(-90.,90,5.),labels=[1,0,0,0])
#     m.drawmeridians(np.arange(180.,360.,5. ),labels=[0,0,0,1])
#     x, y = m(lon, lat)
#     Q = m.quiver(x,y,udaily[ii,:,:],vdaily[ii,:,:])
#     qk = plt.quiverkey(Q, -.08, .8, 1, '1 m/s', labelpos='W',color = 'r')
#     for k in range(0,numdates-1):
#         if (centerdates[k] == oscardates[ii]):
#                 Q = m.quiver(x,y,u[k,:,:],v[k,:,:], color = 'b')
#     plt.title(figtit1 + str(int(oscardates[ii])))
#     plt.savefig(pp, format = "pdf")
# pp.close()
# os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/currents-daily-zoomed.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')
# #
# # Zoomed-in Daily Currents
# pp = PdfPages('currents-daily.pdf')
# figtit1 = 'Daily Averaged OSCAR Data on DOY: '
# from matplotlib.colors import LogNorm
# for ii in range(0,len(oscardates)-5):    #
#     # 
#     import matplotlib.path as mpath
#     import matplotlib.patches as mpatches
#     import matplotlib as mpl
#     from matplotlib.collections import PatchCollection
#     Path = mpath.Path
#     fig=plt.figure(figsize=(8,4.5))
#     m = Basemap(llcrnrlon=285,llcrnrlat=-65,urcrnrlon=310,urcrnrlat=-50,projection='cyl',resolution='i')
#     ax = fig.add_axes([0.05,0.05,0.9,0.85])
#     m.drawcoastlines(linewidth=0.25)
#     m.drawcountries(linewidth=0.25)
#     m.fillcontinents(color='gray',lake_color='white')
#     m.drawmapboundary(fill_color='white')
#     m.drawparallels(np.arange(-90.,90,5.),labels=[1,0,0,0])
#     m.drawmeridians(np.arange(180.,360.,5. ),labels=[0,0,0,1])
#     x, y = m(lon, lat)
#     Q = m.quiver(x,y,udaily[ii,:,:],vdaily[ii,:,:])
#     qk = plt.quiverkey(Q, -.08, .8, 1, '1 m/s', labelpos='W',color = 'r')
#     for k in range(0,numdates-1):
#         if (centerdates[k] == oscardates[ii]):
#                 Q = m.quiver(x,y,u[k,:,:],v[k,:,:], color = 'b')
#     plt.title(figtit1 + str(int(oscardates[ii])))
#     plt.savefig(pp, format = "pdf")
# pp.close()
# os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/currents-daily.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')
# #
# # Zoomed-in Tracer Locations
# pp = PdfPages('currents-tracer-zoom.pdf')
# figtit1 = 'OSCAR 5 Day Avg '
# figtit2 = ' Centered About DOY '
# from matplotlib.colors import LogNorm

# for l in range(0,arraylen-1):
#     #
#     # Find the correct OSCAR data plot
#     k = 0
#     for kk in range(0,len(oscardates)-5):
#         if (oscardates[k] == floatdate[l]):
#             break
#         else:
#             k = k + 1
#     #
#     usefloatlon= floatlon[l]
#     usefloatlat = floatlat[l]
#     usefloatlon_new = floatlon[l+1]
#     usefloatlat_new = floatlat[l+1]
#     usetracerlat_U_V = tracer_newlat_U_V[l]
#     usetracerlon_U_V = tracer_newlon_U_V[l]
#     #
#     import matplotlib.path as mpath
#     import matplotlib.patches as mpatches
#     import matplotlib as mpl
#     from matplotlib.collections import PatchCollection
#     Path = mpath.Path
#     fig=plt.figure(figsize=(8,4.5))
#     m = Basemap(llcrnrlon=297,llcrnrlat=-61,urcrnrlon=307,urcrnrlat=-54,projection='mill',resolution='i')
#     ax = fig.add_axes([0.05,0.05,0.9,0.85])
#     m.drawcoastlines(linewidth=0.25)
#     m.drawcountries(linewidth=0.25)
#     m.fillcontinents(color='gray',lake_color='white')
#     m.drawmapboundary(fill_color='white')
#     m.drawparallels(np.arange(-90.,90,5.),labels=[1,0,0,0])
#     m.drawmeridians(np.arange(180.,360.,5. ),labels=[0,0,0,1])
#     x, y = m(lon, lat)
#     floatx,floaty =m(usefloatlon, usefloatlat)
#     floatx_new,floaty_new =m(usefloatlon_new, usefloatlat_new)
#     tracerx_U_V,tracery_U_V = m(usetracerlon_U_V, usetracerlat_U_V)
#     cs = m.scatter(floatx_new,floaty_new, s = 10, marker=(5,0), c='b')    
#     cs1 = m.scatter(floatx,floaty, s = 10, marker=(5,0), c='k')
#     cs2 = m.scatter(tracerx_U_V,tracery_U_V, s = 10, marker=(5,0), c='r')
#     Q = m.quiver(x,y,udaily[k,:,:],vdaily[k,:,:])
#     qk = plt.quiverkey(Q, -.08, .8, 1, '1 m/s', labelpos='W',color = 'r')
#     plt.title('OSCAR Currents, Float and Tracer on DOY: '  + str(int(oscardates[k])))
#     plt.savefig(pp, format = "pdf")
# pp.close()


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

# #
# # Zoomed-in Daily Currents and MKE
# pp = PdfPages('currents-daily-mke-zoomed.pdf')
# figtit1 = 'OSCAR Daily Averaged MKE on DOY: '
# from matplotlib.colors import LogNorm
# for ii in range(0,len(oscardates)-5):    
#     #
#     # Import Plotting Packages 
#     import matplotlib.path as mpath
#     import matplotlib.patches as mpatches
#     import matplotlib as mpl
#     from matplotlib.collections import PatchCollection
#     Path = mpath.Path
#     #
#     # Setting-up the map
#     fig=plt.figure(figsize=(8,4.5))
#     m = Basemap(llcrnrlon=297,llcrnrlat=-61,urcrnrlon=307,urcrnrlat=-54,projection='cyl',resolution='i')
#     ax = fig.add_axes([0.05,0.05,0.9,0.85])
#     m.drawcoastlines(linewidth=0.25)
#     m.drawcountries(linewidth=0.25)
#     m.fillcontinents(color='gray',lake_color='white')
#     m.drawmapboundary(fill_color='white')
#     m.drawparallels(np.arange(-90.,90,5.),labels=[1,0,0,0])
#     m.drawmeridians(np.arange(180.,360.,5. ),labels=[0,0,0,1])
#     #
#     # Plotting
#     x, y = m(lon, lat)
#     cs = m.contourf(x,y,mke[k,:,:], levels = np.logspace(-1,3.5,500), norm = LogNorm())
#     cb = plt.colorbar(cs, ticks=[1, 10, 100, 1000])
#     cb.ax.set_yticklabels([1, 10,100,1000 ])
#     cb.set_label('Mean Kinetic Energy (cm^2/s^2)')
#     Q = m.quiver(x,y,udaily[ii,:,:],vdaily[ii,:,:])
#     qk = plt.quiverkey(Q, -.08, .8, 1, '1 m/s', labelpos='W',color = 'r')
#     plt.title(figtit1 + str(int(oscardates[ii])))
#     plt.savefig(pp, format = "pdf")
# pp.close()
# os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/currents-daily-mke-zoomed.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')
# #
# # Daily Currents and MKE
# pp = PdfPages('currents-daily-mke.pdf')
# figtit1 = 'OSCAR Daily Averaged MKE on DOY: '
# from matplotlib.colors import LogNorm
# for ii in range(0,len(oscardates)-5):    
#     #
#     # Import Plotting Packages 
#     import matplotlib.path as mpath
#     import matplotlib.patches as mpatches
#     import matplotlib as mpl
#     from matplotlib.collections import PatchCollection
#     Path = mpath.Path
#     #
#     # Setting-up the map
#     fig=plt.figure(figsize=(8,4.5))
#     m = Basemap(llcrnrlon=285,llcrnrlat=-65,urcrnrlon=310,urcrnrlat=-50,projection='cyl',resolution='i')
#     ax = fig.add_axes([0.05,0.05,0.9,0.85])
#     m.drawcoastlines(linewidth=0.25)
#     m.drawcountries(linewidth=0.25)
#     m.fillcontinents(color='gray',lake_color='white')
#     m.drawmapboundary(fill_color='white')
#     m.drawparallels(np.arange(-90.,90,5.),labels=[1,0,0,0])
#     m.drawmeridians(np.arange(180.,360.,5. ),labels=[0,0,0,1])
#     #
#     # Plotting
#     x, y = m(lon, lat)
#     cs = m.contourf(x,y,mke[k,:,:], levels = np.logspace(-1,3.5,500), norm = LogNorm())
#     cb = plt.colorbar(cs, ticks=[1, 10, 100, 1000])
#     cb.ax.set_yticklabels([1, 10,100,1000 ])
#     cb.set_label('Mean Kinetic Energy (cm^2/s^2)')
#     Q = m.quiver(x,y,udaily[ii,:,:],vdaily[ii,:,:])
#     qk = plt.quiverkey(Q, -.08, .8, 1, '1 m/s', labelpos='W',color = 'r')
#     plt.title(figtit1 + str(int(oscardates[ii])))
#     plt.savefig(pp, format = "pdf")
# pp.close()
# os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/currents-daily-mke.pdf ardavies@storm.ceoe.udel.edu:/dev/ardavies/')


# #
# # MKE Plotting
# pp = PdfPages('currents-mke-daily.pdf')
# figtit1 = 'OSCAR 5 Day Avg '
# figtit2 = ' Centered About DOY '
# from matplotlib.colors import LogNorm
# for k in range(0,numdates-3):
#     counter = 0
#     for l in range(0,arraylen-1):
#         if (centerdates[k]-2 <= floatdate[l] <= centerdates[k] +2):
#             counter = counter + 1
#     indexes = np.zeros(counter)
#     counter = 0
#     for l in range(0,arraylen-1):
#         if (centerdates[k]-2 <= floatdate[l] <= centerdates[k] +2):
#             indexes[counter] = l
#             counter = counter + 1
#     usefloatlon = np.zeros(counter)
#     usefloatlat = np.zeros(counter)
#     for n in range(0,counter):
#         nn = int(indexes[n])
#         usefloatlon[n] = floatlon[nn]
#         usefloatlat[n] = floatlat[nn]
#     m = Basemap(llcrnrlon=285,llcrnrlat=-65,urcrnrlon=310,urcrnrlat=-50,projection='mill',resolution='i')
#     fig=plt.figure(figsize=(8,4.5))
#     ax = fig.add_axes([0.05,0.05,0.9,0.85])
#     m.drawcoastlines(linewidth=0.25)
#     m.drawcountries(linewidth=0.25)
#     m.fillcontinents(color='gray',lake_color='white')
#     m.drawmapboundary(fill_color='white')
#     m.drawparallels(np.arange(-90.,90,5.),labels=[1,0,0,0])
#     m.drawmeridians(np.arange(180.,360.,5. ),labels=[0,0,0,1])
#     x, y = m(lon, lat)
#     floatx,floaty =m(usefloatlon, usefloatlat)
#     cs = m.contourf(x,y,mke[k,:,:], levels = np.logspace(-1,3.5,500), norm = LogNorm())
#     cs2 = m.scatter(floatx,floaty, c = 'k', s = 25, marker=(5,0))
#     cs3 = m.plot(floatx,floaty, color = 'k')
#     Q = m.quiver(x,y,u[k,:,:],v[k,:,:])
#     qk = plt.quiverkey(Q, -.08, .8, 1, '1 m/s', labelpos='W',color = 'r')
#     cb = plt.colorbar(cs, ticks=[1, 10, 100, 1000])
#     cb.ax.set_yticklabels([1, 10,100,1000 ])
#     cb.set_label('Mean Kinetic Energy (cm^2/s^2)')
#     plt.title(figtit1 + 'Currents & MKE' + figtit2 + str(int(centerdates[k])))
#     plt.savefig(pp, format = "pdf")
# pp.close()







# os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/currents-daily.pdf ardavies@storm.ceoe.udel.edu:/www/dev/ardavies/')


# #
# # ===================================================================================
# # Pushing Plots to modata
# # ===================================================================================
# #
# os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/currents-tracer-zoom.pdf ardavies@storm.ceoe.udel.edu:/www/dev/ardavies/')
# os.system('scp /home/ardavies/satdata/OSCAR/pdfoutput/ratio.pdf ardavies@storm.ceoe.udel.edu:/www/dev/ardavies/')

os.chdir('/home/ardavies/satdata/OSCAR')