In [1]:

import numpy as np

In [2]:

#w=np.array([-.1,-.1,-.1,-.1,-.1, 0])
#w=np.array([-1,-1,-1,-1,-1, 0])
#c=np.array([  1,  1,  1,.01,100, 1e10])*1e-6
e1e2t=214298.32535430687
e3t=np.array([25.70331479, 26.28684983, 26.59728865, 26.75965336, 
              26.84381704, 26.88724213, 26.90959407, 26.92108493,
              26.9269885,  26.93002054, 26.93157752, 26.93237697,])/20
c=np.array([0.31565186, 0.30992895, 0.3301549,  0.423929,
            0.54217166, 0.01, 0.01, 0.01,
            0.01 , 0.0099, 100,   0.])
w=-1*np.ones(np.shape(c))*2.8e-4
w[-1]=0
tmask=np.ones(np.shape(c))
tmask[-1]=0

In [3]:

#1. bottom value zero
zwz=np.zeros(np.shape(c))
#2. upstream advection
zwz[1:]=w[1:]*c[:-1]
zwz

Out[3]:

array([ 0.00000000e+00, -8.83825208e-05, -8.67801060e-05, -9.24433720e-05,
       -1.18700120e-04, -1.51808065e-04, -2.80000000e-06, -2.80000000e-06,
       -2.80000000e-06, -2.80000000e-06, -2.77200000e-06,  0.00000000e+00])

In [4]:

dt=40*2

In [5]:

ztra=np.zeros(np.shape(c))
ztra[:-1]=-1*(zwz[:-1]-zwz[1:])/e1e2t
delc_upwind=ztra.copy() # ** added to tra at this step
zwi=(e3t*c+dt*ztra)/e3t*tmask # set fse3t's to 1
zwi # upwind increases low point from .01 to .0298 in one time step

Out[5]:

array([3.15651834e-01, 3.09928950e-01, 3.30154898e-01, 4.23928993e-01,
       5.42171651e-01, 1.00000414e-02, 1.00000000e-02, 1.00000000e-02,
       1.00000000e-02, 9.90000001e-03, 1.00000000e+02, 0.00000000e+00])

In [6]:

# 3. antidiffusive flux
zwz_sav=zwz.copy()

In [7]:

jnzts=20 # same as SSC

In [8]:

cb=c.copy()

In [9]:

ca=np.zeros(np.shape(c))
zwzts=np.zeros(np.shape(c))

In [10]:

for jl in range(0,jnzts):
    if jl==0:
        jtaken=np.mod(jnzts+1,2)
        print(jtaken)
        zts=dt/jnzts # z_rzts=1/jnzts
        cn=c.copy() # this makes first step euler
    elif jl==1:
        zts=2*dt/jnzts
    zwz[1:]=.5*w[1:]*(cn[1:]+cn[:-1])
    if jtaken==0:
        zwzts[1:]=zwzts[1:]+zwz[1:]*zts
    jtaken=np.mod(jtaken+1,2) # switch on/off
    
    zbtr=1/(e1e2t*e3t)
    ztra[1:-1]=-1*zbtr[1:-1]*(zwz[1:-1]-zwz[2:])
    ca[1:-1]=cb[1:-1]+zts*ztra[1:-1] # step forward
    # swap for next step:
    cb=cn.copy()
    cn=ca.copy()
    ca=np.zeros(np.shape(c))

In [11]:

cn

Out[11]:

array([0.00000000e+00, 3.09928937e-01, 3.30154896e-01, 4.23928992e-01,
       5.42171676e-01, 1.00000207e-02, 1.00000000e-02, 1.00000000e-02,
       1.00000000e-02, 9.89611895e-03, 1.00000004e+02, 0.00000000e+00])

In [12]:

zwzts

Out[12]:

array([ 0.00000000e+00, -3.47120417e-03, -7.16893903e-03, -8.44573961e-03,
       -1.08203274e-02, -6.18432280e-03, -2.24000116e-04, -2.24000000e-04,
       -2.24000000e-04, -2.22858266e-04, -1.12011088e+00,  0.00000000e+00])

In [13]:

# anti-diffusive vertical flux using average flux from sub-timestepping
zwz=zwzts/dt-zwz_sav

In [14]:

zwz

Out[14]:

array([ 0.00000000e+00,  4.49924686e-05, -2.83163186e-06, -1.31283731e-05,
       -1.65539729e-05,  7.45040298e-05, -1.44821411e-12, -1.86347248e-20,
       -2.75210745e-16,  1.42716732e-08, -1.39986140e-02,  0.00000000e+00])

In [15]:

zbig=1e40
zrtrn=1e-15

In [16]:

# 4. nonosc(pbef=ptb,zwx,zwy,pcc=zwz,paft=zwi,p2dt)
# zwi is after value based on upwind
pbef=c.copy()
paft=zwi.copy()
# search local extrema
zbup=np.maximum(pbef*tmask-zbig*(1-tmask),paft*tmask-zbig*(1-tmask))
zbup

Out[16]:

array([ 3.15651860e-01,  3.09928950e-01,  3.30154900e-01,  4.23929000e-01,
        5.42171660e-01,  1.00000414e-02,  1.00000000e-02,  1.00000000e-02,
        1.00000000e-02,  9.90000001e-03,  1.00000000e+02, -1.00000000e+40])

In [17]:

zbdo=np.minimum(pbef*tmask+zbig*(1-tmask),paft*tmask+zbig*(1-tmask))
zbdo

Out[17]:

array([3.15651834e-01, 3.09928950e-01, 3.30154898e-01, 4.23928993e-01,
       5.42171651e-01, 1.00000000e-02, 1.00000000e-02, 1.00000000e-02,
       1.00000000e-02, 9.90000000e-03, 1.00000000e+02, 1.00000000e+40])

In [18]:

zup=np.zeros(np.shape(c))
zdo=zup.copy()
zup[0]=np.maximum(zbup[0],zbup[1])
zup[1:-1]=np.maximum(np.maximum(zbup[:-2],zbup[1:-1]),zbup[2:])
#zup[-1]=np.maximum(zbup[-2],zbup[-1]) this is never calculated
zup

Out[18]:

array([3.15651860e-01, 3.30154900e-01, 4.23929000e-01, 5.42171660e-01,
       5.42171660e-01, 5.42171660e-01, 1.00000414e-02, 1.00000000e-02,
       1.00000000e-02, 1.00000000e+02, 1.00000000e+02, 0.00000000e+00])

In [19]:

zdo[0]=np.minimum(zbdo[0],zbdo[1])
zdo[1:-1]=np.minimum(np.minimum(zbdo[:-2],zbdo[1:-1]),zbdo[2:])
#zdo[-1]=np.minimum(zbdo[-2],zbdo[-1])
zdo

Out[19]:

array([0.30992895, 0.30992895, 0.30992895, 0.3301549 , 0.01      ,
       0.01      , 0.01      , 0.01      , 0.0099    , 0.0099    ,
       0.0099    , 0.        ])

In [20]:

zeros=np.zeros(np.shape(c))
pcc=zwz.copy()
zpos=np.zeros(np.shape(c))
zneg=zpos.copy()
zpos[:-1]=np.maximum(zeros[1:],pcc[1:])-np.minimum(zeros[:-1],pcc[:-1])
zneg[:-1]=np.maximum(zeros[:-1],pcc[:-1])-np.minimum(zeros[1:],pcc[1:])

In [21]:

zpos

Out[21]:

array([4.49924686e-05, 0.00000000e+00, 2.83163186e-06, 1.31283731e-05,
       9.10580028e-05, 0.00000000e+00, 1.44821411e-12, 1.86347248e-20,
       1.42716734e-08, 0.00000000e+00, 1.39986140e-02, 0.00000000e+00])

In [22]:

zneg

Out[22]:

array([0.00000000e+00, 4.78241005e-05, 1.31283731e-05, 1.65539729e-05,
       0.00000000e+00, 7.45040313e-05, 1.86347248e-20, 2.75210745e-16,
       0.00000000e+00, 1.39986283e-02, 0.00000000e+00, 0.00000000e+00])

In [23]:

#! up & down beta terms
zbt = e1e2t*e3t / dt
zbetup = ( zup - paft) / ( zpos + zrtrn ) * zbt
zbetdo = ( paft - zdo ) / ( zneg + zrtrn ) * zbt

In [24]:

zbetup

Out[24]:

array([1.96438478e+00, 7.12108645e+16, 1.17972908e+08, 3.22806733e+07,
       3.63591818e-01, 1.91644677e+18, 1.02905386e+08, 0.00000000e+00,
       0.00000000e+00, 3.60655436e+20, 0.00000000e+00, 0.00000000e+00])

In [25]:

zbetdo

Out[25]:

array([1.97016635e+16, 3.35064296e-02, 5.48825328e+06, 2.03029860e+07,
       1.91335167e+18, 2.00000003e+00, 0.00000000e+00, 0.00000000e+00,
       3.60650534e+14, 2.00019628e-06, 3.60676287e+20, 0.00000000e+00])

In [26]:

#! monotonic flux in the k direction, i.e. pcc
ones=np.ones(np.shape(c))
za=zeros.copy()
zb=zeros.copy()
zc=zeros.copy()
za[:-1] = np.minimum(np.minimum( ones[1:], zbetdo[1:]), zbetup[:-1] )
zb[:-1] = np.minimum(np.minimum( ones[1:], zbetup[1:]), zbetdo[:-1] )
zc[:-1]=np.where(pcc[1:]<0,0,1)
#zc[:-1]=.5+np.abs(pcc[1:])/pcc[1:]*.5
#zc =       ( 0.5  + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
# SIGN(A,B) returns value of A with sign of B
pcc[1:]=pcc[1:]*(zc[:-1]*za[:-1]+(1-zc[:-1])*zb[:-1])
#pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb )

In [27]:

pcc

Out[27]:

array([ 0.00000000e+00,  1.50753698e-06, -9.48778735e-08, -1.31283731e-05,
       -6.01888912e-06,  2.70890557e-05, -1.44821411e-12, -0.00000000e+00,
       -0.00000000e+00,  0.00000000e+00, -0.00000000e+00,  0.00000000e+00])

In [28]:

# before nonosc, zwz contains difference between centered substepped leapfrog and upwind:
zwz

Out[28]:

array([ 0.00000000e+00,  4.49924686e-05, -2.83163186e-06, -1.31283731e-05,
       -1.65539729e-05,  7.45040298e-05, -1.44821411e-12, -1.86347248e-20,
       -2.75210745e-16,  1.42716732e-08, -1.39986140e-02,  0.00000000e+00])

In [29]:

# after nonosc, zwz (pcc) is adjusted so that only one element changes from upwind 

In [30]:

# 5. final trend with corrected fluxes
ztra=0*ztra
ztra2=ztra.copy()
ztra[:-1]=-1*zbtr[:-1]*(pcc[:-1]-pcc[1:])
ztra2[:-1]=-1*zbtr[:-1]*(zwz[:-1]-zwz[1:])

In [31]:

# total change is ztra+delc, ztra is after nonosc adjustment
ztra

Out[31]:

array([ 5.47381429e-12, -5.68915331e-12, -4.57335326e-11,  2.47952685e-11,
        1.15106320e-10, -9.40283574e-11,  5.02269547e-18, -0.00000000e+00,
        0.00000000e+00, -0.00000000e+00,  0.00000000e+00,  0.00000000e+00])

In [32]:

ztra2 # this is before nonosc version

Out[32]:

array([ 1.63366087e-10, -1.69792884e-10, -3.61304733e-11, -1.19472339e-11,
        3.16581162e-10, -2.58609654e-10,  5.02269541e-18, -9.54013762e-22,
        4.94650364e-14, -4.85132104e-08,  4.85103562e-08,  0.00000000e+00])

In [33]:

delc_upwind

Out[33]:

array([-4.12427492e-10,  7.47749567e-12, -2.64270194e-11, -1.22524280e-10,
       -1.54494650e-10,  6.95330048e-10, -0.00000000e+00, -0.00000000e+00,
       -0.00000000e+00,  1.30658977e-13,  1.29352387e-11,  0.00000000e+00])

In [34]:

delc_upwind+ztra

Out[34]:

array([-4.06953677e-10,  1.78834236e-12, -7.21605520e-11, -9.77290114e-11,
       -3.93883298e-11,  6.01301690e-10,  5.02269547e-18, -0.00000000e+00,
        0.00000000e+00,  1.30658977e-13,  1.29352387e-11,  0.00000000e+00])

In [35]:

delc_upwind+ztra2

Out[35]:

array([-2.49061405e-10, -1.62315388e-10, -6.25574927e-11, -1.34471514e-10,
        1.62086512e-10,  4.36720394e-10,  5.02269541e-18, -9.54013762e-22,
        4.94650364e-14, -4.85130797e-08,  4.85232915e-08,  0.00000000e+00])

In [36]:

Out[36]:

array([3.1565186e-01, 3.0992895e-01, 3.3015490e-01, 4.2392900e-01,
       5.4217166e-01, 1.0000000e-02, 1.0000000e-02, 1.0000000e-02,
       1.0000000e-02, 9.9000000e-03, 1.0000000e+02, 0.0000000e+00])

In [37]:

# nonosc and diatoms???? maybe code shouldn't have to be monotonic