import sys
import readParams_moreoptions as rdp1
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
from mom_plot1 import m6plot, xdegtokm
import numpy as np
from netCDF4 import MFDataset as mfdset, Dataset as dset
import time
from plot_twamomx_budget_complete_direct_newest import extract_twamomx_terms
from plot_twamomy_budget_complete_direct_newest import extract_twamomy_terms
import pyximport
pyximport.install()
from getvaratzc import getvaratzc5, getvaratzc

def getutwa(fhgeo, fh, fh2, sl):
    dt = fh.variables['average_DT'][:]
    dt = dt[:,np.newaxis,np.newaxis,np.newaxis]
    uh = (fh.variables['uh_masked'][sl].filled(0)*dt).sum(axis=0,keepdims=True)/np.sum(dt)
    h_cu = fh.variables['h_Cu'][sl].filled(0).mean(axis=0,keepdims=True)
    h_cu = np.ma.masked_array(h_cu,mask=(h_cu<1e-3))
    dycu = fhgeo.variables['dyCu'][sl[2:]]
    utwa = uh/h_cu/dycu
    return utwa, h_cu

def getvtwa(fhgeo, fh, fh2, sl):
    dt = fh.variables['average_DT'][:]
    dt = dt[:,np.newaxis,np.newaxis,np.newaxis]
    vh = (fh.variables['vh_masked'][sl]*dt).sum(axis=0,keepdims=True)/np.sum(dt)
    h_cv = fh.variables['h_Cv'][sl].mean(axis=0,keepdims=True)
    h_cv = np.ma.masked_array(h_cv,mask=(h_cv<1e-3))
    dxcv = fhgeo.variables['dxCv'][sl[2:]]
    vtwa = vh/dxcv/h_cv
    vtwa = np.concatenate((vtwa,-vtwa[:,:,:,-1:]),axis=3)
    h_cv = np.concatenate((h_cv,h_cv[:,:,:,-1:]),axis=3)
    return vtwa, h_cv

def getpv(fhgeo, fh, fh2, xs, xe, ys, ye, zs=0, ze=None):
    sl = np.s_[:,zs:ze,ys:ye,xs:xe]
    slpy = np.s_[:,zs:ze,ys:ye+1,xs:xe]

    utwa,h_cu = getutwa(fhgeo, fh, fh2, slpy)
    dybu = fhgeo.variables['dyBu'][sl[2:]]
    utway = np.diff(utwa,axis=2)/dybu

    vtwa,h_cv = getvtwa(fhgeo, fh, fh2, sl)
    dxbu = fhgeo.variables['dxBu'][sl[2:]]
    vtwax = np.diff(vtwa,axis=3)/dxbu

    h_q = 0.25*(h_cu[:,:,:-1,:] + h_cu[:,:,1:,:] +
            h_cv[:,:,:,:-1] + h_cv[:,:,:,1:])
    f = fhgeo.variables['f'][sl[2:]]
    pvhash = (f - utway + vtwax)/h_q
    return pvhash, h_q

def extract_twapv_terms(geofil,vgeofil,fil,fil2,xstart,xend,ystart,yend,zs,ze,meanax,fil3=None,
      alreadysaved=False):

    if not alreadysaved:
        keepax = ()
        for i in range(4):
            if i not in meanax:
                keepax += (i,)

        fhgeo = dset(geofil)
        fh = mfdset(fil)
        fh2 = mfdset(fil2)
        zi = rdp1.getdims(fh)[2][0]
        dbl = -np.diff(zi)*9.8/1031
        (xs,xe),(ys,ye),dimq = rdp1.getlatlonindx(fh,wlon=xstart,elon=xend,
                slat=ystart, nlat=yend,zs=zs,ze=ze,xhxq='xq',yhyq='yq')
        dxbu = rdp1.getgeombyindx(fhgeo,xs,xe,ys,ye)[2][4]
        dybu = rdp1.getgeombyindx(fhgeo,xs,xe,ys,ye)[2][5]
        aq = rdp1.getgeombyindx(fhgeo,xs,xe,ys,ye)[1][1]
        f = rdp1.getgeombyindx(fhgeo,xs,xe,ys,ye)[-1]
        nt_const = dimq[0].size

        pvhash,hq = getpv(fhgeo, fh, fh2, xs, xe, ys, ye)
        sl = np.s_[:,zs:ze,ys:ye,xs:xe]
        slpy = np.s_[:,zs:ze,ys:ye+1,xs:xe]
        dxcu = fhgeo.variables['dxCu'][slpy[2:]]
        dycv = fhgeo.variables['dyCv'][sl[2:]]
        dycv = np.concatenate((dycv,dycv[:,-1:]),axis=1)

        if fil3:
            fh3 = mfdset(fil3)
            sltn = np.s_[-1:,zs:ze,ys:ye,xs:xe]
            islayerdeep0 = fh3.variables['islayerdeep'][-1:,0,0,0]
            islayerdeep = (fh3.variables['islayerdeep'][sltn].filled(np.nan))
            swash = (islayerdeep0 - islayerdeep)/islayerdeep0*100
            fh3.close()
        else:
            swash = None

        xmom = extract_twamomx_terms(geofil,vgeofil,fil,fil2,xs,xe,ys,ye+1,zs,ze,(0,),
                alreadysaved=False,xyasindices=True,calledfrompv=True)[2]
        ymom = extract_twamomy_terms(geofil,vgeofil,fil,fil2,xs,xe,ys,ye,zs,ze,(0,),
                alreadysaved=False,xyasindices=True,calledfrompv=True)[2]

        xmom = xmom[np.newaxis,:,:,:,:]
        ymom = ymom[np.newaxis,:,:,:,:]
        ymom = np.concatenate((ymom,ymom[:,:,:,-1:]),axis=3)

        bxppvflx = np.sum(xmom[:,:,:,:,[0,1,3,4]],axis=4)
        pvhash1,_ = getpv(fhgeo, fh, fh2, xs, xe, ys-1, ye+1)
        sl1 = np.s_[:,zs:ze,ys-1:ye+1,xs:xe]
        vtwa1, h_cv1 = getvtwa(fhgeo,fh,fh2,sl1)
        vtwa1 = 0.5*(vtwa1[:,:,:,:-1] + vtwa1[:,:,:,1:])
        pvhashvtwa = pvhash1*vtwa1
        sl1 = np.s_[:,zs:ze,ys:ye+1,xs:xe]
        h_cu1 = fh.variables['h_Cu'][sl1].filled(np.nan).mean(axis=0,keepdims=True)
        h_cu1 = np.ma.masked_array(h_cu1,mask=(h_cu1<1e-3))
        pvflxx = h_cu1*(pvhashvtwa[:,:,:-1,:]+pvhashvtwa[:,:,1:,:])/2
        
        byppvflx = np.sum(ymom[:,:,:,:,[0,1,3,4]],axis=4)
        pvhash1,_ = getpv(fhgeo, fh, fh2, xs-1, xe, ys, ye)
        sl1 = np.s_[:,zs:ze,ys:ye+1,xs-1:xe]
        utwa1, h_cu1 = getutwa(fhgeo,fh,fh2,sl1)
        utwa1 = 0.5*(utwa1[:,:,:-1,:]+utwa1[:,:,1:,:])
        pvhashutwa = pvhash1*utwa1
        pvhashutwa[:,:,:,-1] = 0.0
        sl1 = np.s_[:,zs:ze,ys:ye,xs:xe]
        h_cv1 = fh.variables['h_Cv'][sl1].mean(axis=0,keepdims=True)
        h_cv1 = np.ma.masked_array(h_cv1,mask=(h_cv1<1e-3))
        pvflxy = h_cv1*(pvhashutwa[:,:,:,:-1]+pvhashutwa[:,:,:,1:])/2
        pvflxy = np.concatenate((pvflxy,pvflxy[:,:,:,-1:]),axis=3)

        bx = bxppvflx - pvflxx
        by = byppvflx + pvflxy

        xmom1 = xmom[:,:,:,:,[2,5,6,7,8,9]]
        xmom1 = np.concatenate((+pvflxx[:,:,:,:,np.newaxis], xmom1, bx[:,:,:,:,np.newaxis]), axis = 4)
        ymom1 = ymom[:,:,:,:,[2,5,6,7,8,9]]
        ymom1 = np.concatenate((-pvflxy[:,:,:,:,np.newaxis], ymom1, by[:,:,:,:,np.newaxis]), axis = 4) 

        #pv = (-np.diff(xmom*dxcu[:,:,np.newaxis],axis=2) + np.diff(ymom*dycv[:,:,np.newaxis],axis=3))/aq[:,:,np.newaxis]
        pv = -np.diff(xmom,axis=2)/dybu[:,:,np.newaxis] + np.diff(ymom,axis=3)/dxbu[:,:,np.newaxis]
        pv1x = -np.diff(xmom1*dxcu[:,:,np.newaxis],axis=2)/aq[:,:,np.newaxis]
        pv1y =  np.diff(ymom1*dycv[:,:,np.newaxis],axis=3)/aq[:,:,np.newaxis]

        slyp = np.s_[:,:,ys:ye+1,xs:xe]
        ah1 = fhgeo.variables['Ah'][slyp[2:]]
        slxmyp = np.s_[:,:,ys:ye+1,xs-1:xe]
        uh = fh2.variables['uh'][slxmyp].filled(0).mean(axis=0,keepdims=True)
        uhx = np.diff(uh,axis=3)/ah1
        uhx = np.concatenate((uhx,uhx[:,:,:,-1:]),axis=3)
        uhx = 0.25*(uhx[:,:,:-1,:-1] + uhx[:,:,:-1,1:] + uhx[:,:,1:,:-1] +
                uhx[:,:,1:,1:])
        pv1y[:,:,:,:,0] += pvhash*uhx

        slymp = np.s_[:,:,ys-1:ye+1,xs:xe]
        vh = fh2.variables['vh'][slymp].mean(axis=0,keepdims=True)
        vhy = np.diff(vh,axis=2)/ah1
        vhy = np.concatenate((vhy,vhy[:,:,:,-1:]),axis=3)
        vhy = 0.25*(vhy[:,:,:-1,:-1] + vhy[:,:,:-1,1:] + vhy[:,:,1:,:-1] +
                vhy[:,:,1:,1:])
        pv1x[:,:,:,:,0] += pvhash*vhy

        wd = fh2.variables['wd'][slyp].mean(axis=0,keepdims=True)
        wdb = np.diff(wd,axis=1)
        wdb = np.concatenate((wdb,wdb[:,:,:,-1:]),axis=3)
        wdb = 0.25*(wdb[:,:,:-1,:-1] + wdb[:,:,:-1,1:] + wdb[:,:,1:,:-1] +
                wdb[:,:,1:,1:])
        pv3 = pvhash*wdb
        pv3 = pv3[:,:,:,:,np.newaxis]
        #hq[hq<1] = np.nan
        pvnew = np.concatenate((pv1y[:,:,:,:,:1],
                                pv1x[:,:,:,:,:1],
                                pv3,
                                pv1x[:,:,:,:,1:-1],
                                pv1y[:,:,:,:,1:-1],
                                pv1x[:,:,:,:,-1:]+pv1y[:,:,:,:,-1:]),axis=4)/hq[:,:,:,:,np.newaxis]

        pv = np.ma.filled(pv.astype(np.float64), np.nan)
        pvnew = np.ma.filled(pvnew.astype(np.float64), np.nan)
        pvhash = np.ma.filled(pvhash.astype(np.float64), np.nan)

        pv = np.nanmean(pv,meanax,keepdims=True)
        pvnew = np.nanmean(pvnew,meanax,keepdims=True)
        pvhash = np.nanmean(pvhash,meanax,keepdims=True)

        X = dimq[keepax[1]]
        Y = dimq[keepax[0]]
        if 1 in keepax:
            dt = fh.variables['average_DT'][:]
            dt = dt[:,np.newaxis,np.newaxis,np.newaxis]
            em = (fh2.variables['e'][0:,zs:ze,ys:ye,xs:xe]*dt).sum(axis=0,keepdims=True)/np.sum(dt)
            em = np.nanmean(em, meanax,keepdims=True)
            z = np.linspace(-3000,0,100)
            Y = z
            P = getvaratzc5(pv.astype(np.float32),
                    z.astype(np.float32),
                    em.astype(np.float32)).squeeze()
            Pnew = getvaratzc5(pvnew.astype(np.float32),
                    z.astype(np.float32),
                    em.astype(np.float32)).squeeze()
            pvhash = getvaratzc(pvhash.astype(np.float32),
                    z.astype(np.float32),
                    em.astype(np.float32)).squeeze()
            if fil3:
                swash = np.nanmean(swash, meanax,keepdims=True)
                swash = getvaratzc(swash.astype(np.float32),
                        z.astype(np.float32),
                        em.astype(np.float32)).squeeze()
        else:
            P = pv.squeeze()
            Pnew = pvnew.squeeze()
            pvhash = pvhash.squeeze()
            swash = swash.squeeze()
        np.savez('twapv_complete_terms', X=X,Y=Y,P=P)
    else:
        npzfile = np.load('twapv_complete_terms.npz')
        X = npzfile['X']
        Y = npzfile['Y']
        P = npzfile['P']
        
    fhgeo.close()
    fh.close()
    fh2.close()
    return (X,Y,P,pvhash,Pnew,swash)

def plot_twapv(geofil,vgeofil,fil,fil2,xstart,xend,ystart,yend,zs,ze,meanax,
        fil3=None,plotterms = [0,1,10,11,12,13], swashperc = 1,
        cmaxpercfactor = 1,cmaxpercfactorpvhash=15,cmaxpercfactorPnew=15, savfil=None,savfilep=None,alreadysaved=False):
    X,Y,P,pvhash,Pnew,swash = extract_twapv_terms(geofil,vgeofil,fil,fil2,
            xstart,xend,ystart,yend,zs,ze,meanax,alreadysaved=alreadysaved,fil3=fil3)
    cmax = np.nanpercentile(P,[cmaxpercfactor,100-cmaxpercfactor])
    cmax = np.max(np.fabs(cmax))
    fig,ax = plt.subplots(np.int8(np.ceil(P.shape[-1]/2)),2,
                          sharex=True,sharey=True,figsize=(12, 9))
    ti = ['(a)','(b)','(c)','(d)','(e)','(f)','(g)','(h)',
            '(i)','(j)','(k)','(l)','(m)','(n)','(o)','(p)','(q)','(r)']
    labx = [ r'$(\hat{u}\hat{u}_{\tilde{x}})_{\tilde{y}}$', 
            r'$(\hat{v}\hat{u}_{\tilde{y}})_{\tilde{y}}$', 
            r'$(\hat{\varpi}\hat{u}_{\tilde{b}})_{\tilde{y}}$', 
            r'$(-f\hat{v})_{\tilde{y}}$', 
            r'$(\overline{m_{\tilde{x}}})_{\tilde{y}}$', 
            r"""$(\frac{1}{\overline{h}}(\overline{h}\widehat{u^{\prime \prime}u^{\prime \prime}}+\frac{1}{2}\overline{\zeta^{\prime 2}})_{\tilde{x}})_{\tilde{y}}$""", 
            r"""$(\frac{1}{\overline{h}}(\overline{h}\widehat{u^{\prime \prime}v^{\prime \prime}})_{\tilde{y}}$""",
            r"""$(\frac{1}{\overline{h}}(\overline{h}\widehat{u^{\prime \prime}\varpi^{\prime \prime}} + \overline{\zeta^{\prime}m_{\tilde{x}}^{\prime}})_{\tilde{b}})_{\tilde{y}}$""",
            r'$(-\widehat{X^H})_{\tilde{y}}$', 
            r'$(-\widehat{X^V})_{\tilde{y}}$']
    laby = [ r'$(-\hat{u}\hat{v}_{\tilde{x}})_{\tilde{x}}$', 
            r'$(-\hat{v}\hat{v}_{\tilde{y}})_{\tilde{x}}$', 
            r'$(-\hat{\varpi}\hat{v}_{\tilde{b}})_{\tilde{x}}$', 
            r'$(-f\hat{u})_{\tilde{x}}$', 
            r'$(-\overline{m_{\tilde{y}}})_{\tilde{x}}$', 
            r"""$(-\frac{1}{\overline{h}}(\overline{h}\widehat{u^{\prime \prime}v^{\prime \prime}})_{\tilde{x}})_{\tilde{x}}$""", 
            r"""$(-\frac{1}{\overline{h}}(\overline{h}\widehat{v^{\prime \prime}v^{\prime \prime}}+\frac{1}{2}\overline{\zeta^{\prime 2}})_{\tilde{y}})_{\tilde{x}}$""",
            r"""$(-\frac{1}{\overline{h}}(\overline{h}\widehat{v^{\prime \prime}\varpi^{\prime \prime}} + \overline{\zeta^{\prime}m_{\tilde{y}}^{\prime}})_{\tilde{b}})_{\tilde{x}}$""",
            r'$(\widehat{Y^H})_{\tilde{x}}$', 
            r'$(\widehat{Y^V})_{\tilde{x}}$']
    for i in range(P.shape[-1]):
        axc = ax.ravel()[i]
        im = m6plot((X,Y,P[:,:,i]),axc,vmax=cmax,vmin=-cmax,ptype='imshow',
                txt=labx[i]+' + '+laby[i], ylim=(-2500,0),
                cmap='RdBu_r', cbar=False)
        if fil3:
            cs = axc.contour(X,Y,swash,np.array([swashperc]), colors='k')
        
        if i % 2 == 0:
            axc.set_ylabel('z (m)')
        if i > np.size(ax)-3:
            xdegtokm(axc,0.5*(ystart+yend))
            
    fig.tight_layout()
    cb = fig.colorbar(im, ax=ax.ravel().tolist())
    cb.formatter.set_powerlimits((0, 0))
    cb.update_ticks() 
    if savfil:
        plt.savefig(savfil+'.eps', dpi=300, facecolor='w', edgecolor='w', 
                    format='eps', transparent=False, bbox_inches='tight')
    else:
        plt.show()

    im = m6plot((X,Y,np.sum(P,axis=2)),vmax=cmax,vmin=-cmax,ptype='imshow',cmap='RdBu_r',ylim=(-2500,0))
    if savfil:
        plt.savefig(savfil+'res.eps', dpi=300, facecolor='w', edgecolor='w', 
                    format='eps', transparent=False, bbox_inches='tight')
    else:
        plt.show()

    fig,ax = plt.subplots(np.int8(np.ceil(len(plotterms)/2)),2,
                          sharex=True,sharey=True,figsize=(12,7))

    cmaxpvhash = np.nanpercentile(pvhash,
            [cmaxpercfactorpvhash,100-cmaxpercfactorpvhash])
    cmaxpvhash = np.max(np.fabs(cmaxpvhash))

    cmax = np.nanpercentile(Pnew,
            [cmaxpercfactorPnew,100-cmaxpercfactorPnew])
    cmax = np.max(np.fabs(cmax))

    lab = [ r"$-\hat{u}\Pi^{\#}_{\tilde{x}}$",
            r"$-\hat{v}\Pi^{\#}_{\tilde{y}}$",
            r"$\Pi^{\#}(\bar{h} \hat{\varpi})_{\tilde{b}}$",
            r"$\frac{(\hat{\varpi}\hat{u}_{\tilde{b}})_{\tilde{y}}}{\bar{h}}$",
            r"""$\frac{1}{\bar{h}}(\frac{1}{\bar{h}}(\bar{h}\widehat{u^{\prime \prime}u^{\prime \prime}}+\frac{1}{2}\overline{\zeta^{\prime 2}})_{\tilde{x}})_{\tilde{y}}$""", 
            r"""$\frac{1}{\bar{h}}(\frac{1}{\bar{h}}(\bar{h}\widehat{u^{\prime \prime}v^{\prime \prime}})_{\tilde{y}})_{\tilde{y}}$""",
            r"""$\frac{1}{\bar{h}}(\frac{1}{\bar{h}}(\bar{h}\widehat{u^{\prime \prime}\varpi^{\prime \prime}} + \overline{\zeta^{\prime}m_{\tilde{x}}^{\prime}})_{\tilde{b}})_{\tilde{y}}$""",
            r'$-\frac{1}{\bar{h}}(\widehat{X^H})_{\tilde{y}}$', 
            r'$-\frac{1}{\bar{h}}(\widehat{X^V})_{\tilde{y}}$',
            r'$-\frac{(\hat{\varpi}\hat{v}_{\tilde{b}})_{\tilde{x}}}{\bar{h}}$', 
            r"""$-\frac{1}{\bar{h}}(\frac{1}{\bar{h}}(\bar{h}\widehat{u^{\prime \prime}v^{\prime \prime}})_{\tilde{x}})_{\tilde{x}}$""", 
            #r"""$-\frac{1}{\bar{h}}(\frac{1}{\bar{h}}(\bar{h}\widehat{v^{\prime \prime}v^{\prime \prime}}+\frac{1}{2}\overline{\zeta^{\prime 2}})_{\tilde{y}})_{\tilde{x}}$""",
            r"""$-\frac{1}{\bar{h}}(\frac{1}{\bar{h}}(\bar{h}\widehat{v^{\prime \prime}v^{\prime \prime}})_{\tilde{y}})_{\tilde{x}}$""",
            #r"""$-\frac{1}{\bar{h}}(\frac{1}{\bar{h}}(\bar{h}\widehat{v^{\prime \prime}\varpi^{\prime \prime}} + \overline{\zeta^{\prime}m_{\tilde{y}}^{\prime}})_{\tilde{b}})_{\tilde{x}}$""",
            r"""$-\frac{1}{\bar{h}}(\frac{1}{\bar{\zeta_{\tilde{b}}}}(\overline{\zeta^{\prime}m_{\tilde{y}}^{\prime}})_{\tilde{b}})_{\tilde{x}}$""",
            r'$\frac{1}{\bar{h}}(\widehat{Y^H})_{\tilde{x}}$', 
            r'$\frac{1}{\bar{h}}(\widehat{Y^V})_{\tilde{x}}$',
            r'$B_{\tilde{x} \tilde{y}} - B_{\tilde{y} \tilde{x}}$']

    matplotlib.rcParams['contour.negative_linestyle'] = 'solid'
    for i,p in enumerate(plotterms):
        axc = ax.ravel()[i]
        im = m6plot((X,Y,Pnew[:,:,p]),axc,vmax=cmax,vmin=-cmax,ptype='imshow',
                ylim=(-1200,0), txt=lab[p], 
                cmap='RdBu_r', cbar=False)
        im2 = axc.contour(X,Y,pvhash,np.logspace(-6,-5.5,5),colors='grey',linewidths=2)
        im2.clabel(inline=True,fmt="%.1e")
        if fil3:
            cs = axc.contour(X,Y,swash,np.array([swashperc]), colors='k')
        if i % 2 == 0:
            axc.set_ylabel('z (m)')
        if i > np.size(ax)-3:
            xdegtokm(axc,0.5*(ystart+yend))
            
    fig.tight_layout()
    cb = fig.colorbar(im, ax=ax.ravel().tolist())
    cb.formatter.set_powerlimits((0, 0))
    cb.update_ticks() 
    if savfil:
        plt.savefig(savfil+'Pnew.eps', dpi=300, facecolor='w', edgecolor='w', 
                    format='eps', transparent=False, bbox_inches='tight')
    else:
        plt.show()
    im = m6plot((X,Y,np.sum(Pnew,axis=2)),ptype='imshow',vmax=cmax,vmin=-cmax,cmap='RdBu_r',ylim=(-2500,0))
    if savfil:
        plt.savefig(savfil+'Pnewres.eps', dpi=300, facecolor='w', edgecolor='w', 
                    format='eps', transparent=False, bbox_inches='tight')
    else:
        plt.show()

    fig,ax = plt.subplots(1,2,sharex=True,sharey=True,figsize=(10, 3))
    im = m6plot((X,Y,np.nansum(Pnew[:,:,:2],axis=2)),
            ax[0],vmax=cmax,vmin=-cmax,ptype='imshow',
            ylim=(-1200,0), txt=lab[0]+lab[1], 
            cmap='RdBu_r', cbar=False)
    im = m6plot((X,Y,np.nansum(Pnew[:,:,12:13],axis=2)),
            ax[1],vmax=cmax,vmin=-cmax,ptype='imshow',
            ylim=(-1200,0), txt=lab[12], 
            cmap='RdBu_r', cbar=False)
    ax[0].set_ylabel('z (m)')
    for axc in ax:
        xdegtokm(axc,0.5*(ystart+yend))
        im2 = axc.contour(X,Y,pvhash,np.logspace(-6,-5.5,5),colors='grey',linewidths=2)
        im2.clabel(inline=True,fmt="%.1e")
        axc.set_ylim(-1200,0)
    fig.tight_layout()
    cb = fig.colorbar(im, ax=ax.ravel().tolist())
    cb.formatter.set_powerlimits((0, 0))
    cb.update_ticks() 
    if savfil:
        plt.savefig(savfil+'Pnewnew.eps', dpi=300, facecolor='w', edgecolor='w', 
                    format='eps', transparent=False, bbox_inches='tight')
    else:
        plt.show()

    cmax = np.nanpercentile(pvhash,
            [cmaxpercfactorpvhash,100-cmaxpercfactorpvhash])
    cmax = np.max(np.fabs(cmax))
    im = m6plot((X,Y,pvhash),ptype='imshow',vmax=cmax,vmin=-cmax,cmap='RdBu_r',ylim=(-2500,0))
    if fil3:
        cs = axc.contour(X,Y,swash,np.array([swashperc]), colors='k')
    if savfil:
        plt.savefig(savfil+'pvhash.eps', dpi=300, facecolor='w', edgecolor='w', 
                    format='eps', transparent=False, bbox_inches='tight')
    else:
        plt.show()
