notebook.community

Edit and run



In [1]:

    
'''
Test to integrate XAD quantities in NEMO from coastline
After discussion with Sarah Wakelin

'''









    Out[1]:





'\nTest to integrate XAD quantities in NEMO from coastline\nAfter discussion with Sarah Wakelin\n\n'



In [2]:

    
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')



In [3]:

    
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.colors as colors
%matplotlib inline
import math



In [4]:

    
%ls c:\Users\tams00\Documents\nerc_ssb\c_fluxes\AMM7-HINDCAST-v0-erosion\1981\01\









    



 Volume in drive C is Windows
 Volume Serial Number is 228F-02A7

 Directory of c:\Users\tams00\Documents\nerc_ssb\c_fluxes\AMM7-HINDCAST-v0-erosion\1981\01

04/06/2017  14:56    <DIR>          .
04/06/2017  14:56    <DIR>          ..
01/02/2017  15:42           548,870 amm7_1m_19810101_19810131_grid_T.nc
04/06/2017  14:54             3,944 amm7_1m_19810101_19810131_grid_T.nc.aux.xml
01/02/2017  15:42        24,558,914 amm7_1m_19810101_19810131_grid_W.nc
04/06/2017  14:56             1,416 amm7_1m_19810101_19810131_grid_W.nc.aux.xml
01/02/2017  15:50       869,269,937 amm7_1m_19810101_19810131_ptrc_T.nc
04/06/2017  14:56             1,382 amm7_1m_19810101_19810131_ptrc_T.nc.aux.xml
               6 File(s)    894,384,463 bytes
               2 Dir(s)  51,350,167,552 bytes free



In [5]:

    
modelpaths='c:/Users/tams00/Documents/nerc_ssb/c_fluxes/AMM7-HINDCAST-v0-erosion/1981/01/amm7_1m_19810101_19810131_grid_T.nc'
modelpaths='c:/Users/tams00/Documents/nerc_ssb/selected_outputs/amm7_1d_20050501_20050531_grid_T_time_counter0.nc'
bathy_path='c:/Users/tams00/Documents/nerc_ssb/selected_outputs/bathy_meter.nc'
xadv_3d='XAD_O3_c_e3t'
yadv_3d='YAD_O3_c_e3t'
bathy = xr.open_dataset(bathy_path)['Bathymetry'].squeeze()
modelout = xr.open_mfdataset([modelpaths])

modelout;

Functions



In [6]:

    
def is_land(bathy,row,col):
    if math.isclose(bathy.values[row,col], 0.0):
        return True
    else:
        return False



In [7]:

    
# integrate varArr from row,col in the given direction
# and save in intArr
#
# XAD[r,c] = adv_x[r,c+1] - adv_x[r,c]       (tested empirically in embayment case)
# adv_x[r,c] = adv_x[r,c+1] - XAD[r,c]       use when going west
# adv_x[r,c+1] = XAD[r,c] + adv_x[r,c]       use when going east
#
# NOTE: imshow shows geographic array upside down
# YAD[r,c] = adv_y[r+1,c] - adv_y[r,c]       (tested empirically in embayment case)
# adv_x[r,c] = adv_x[r+1,c] - YAD[r,c]       use when going south
# adv_x[r+1,c] = YAD[r,c] + adv_x[r,c]       use when going north


# Corrected Arakawa C indexing in NEMO
# https://www.nemo-ocean.eu/wp-content/uploads/NEMO_book.pdf pp.54
#
#    V
#    _
# j  T |U
#    i
#
# XAD[r,c] = adv_x[r,c] - adv_x[r,c-1]       
# adv_x[r,c] = adv_x[r,c-1] + XAD[r,c]       use when going east
# adv_x[r,c] = adv_x[r,c+1] - XAD[r,c+1]     use when going west
#
# YAD[r,c] = adv_y[r,c] - adv_y[r-1,c]       
# adv_y[r,c] = adv_y[r-1,c] + YAD[r,c]       use when going north
# adv_y[r,c] = adv_y[r+1,c] - YAD[r+1,c]     use when going south

# S.W.: as below, which results in opposite direction (postive right and up?)
# YAD[r,c] = adv_y[r-1,c] - adv_y[r,c]
# adv_y[r,c] = adv_y[r-1,c] - YAD[r,c]       use when going north
# adv_y[r,c] = adv_y[r+1,c] + YAD[r+1,c]     use when going south



def integrate_vector(varArr, intArr, bathy, row,col,direction):
    #print('integrate ',row,col,direction)
    [rows,cols]=varArr.shape
    
    if direction=='west':
        for c in range(col,0,-1):
            if is_land(bathy,row,c):
                break
            if not(math.isclose( intArr[row,c],0. )):
                #print("Break in integrate west. Already integrated east!",row,c)
                break
            # Eastern open boundary
            if ( c == cols-2 ):
                intArr[row,c] = +1 * varArr[row,c]
            else:                
                intArr[row,c] = intArr[row,c+1] - varArr[row,c+1] 
            #print('west',row,c,varArr[row,c+1],intArr[row,c])
        #print(row,c,intArr[row,c])
        
    if direction=='east':
        for c in range(col,cols-1):
            if is_land(bathy,row,c):
                break
            # Western open boundary
            if ( c == 1 ):
                intArr[row,c] = +1 * varArr[row,c]
            else:
                intArr[row,c] = varArr[row,c] + intArr[row,c-1]
            #print('east',row,c,varArr[row,c],intArr[row,c])
        #print(row,c,intArr[row,c])
        
    if direction=='south':
        for r in range(row,0,-1):
            if is_land(bathy,r,col):
                break
            if not(math.isclose( intArr[r,col],0. )):
                #print("Break in integrate south. Already integrated north!",r,col)
                break
            # Northern open boundary
            if ( r == rows-2 ):
                intArr[r,col] = +1 * varArr[r,col]
            else:
                intArr[r,col] = intArr[r+1,col] - varArr[r+1,col]
            #print('south adv-YAD=adv',r,col,intArr[r+1,col],varArr[r+1,col],intArr[r,col],is_land(bathy,r,col),bathy[r,col])
        #print(r,col,intArr[r,col])
            
    if direction=='north':
        for r in range(row,rows-1):
            if is_land(bathy,r,col):
                break
            # Southern open boundary
            if ( r == 1 ):
                intArr[r,col] = +1 * varArr[r,col]
            else:
                intArr[r,col] = varArr[r,col] + intArr[r-1,col]
            #print('north',r,col,varArr[r,col],intArr[r,col])
        #print(r,col,intArr[r,col])

    return intArr



In [8]:

    
def integrate_xadv(varArr,bathy):
    # for XAD
    [rows,cols]=varArr.shape
    horArr=np.zeros([rows,cols])
    for row in range(1,rows-1):
        for col in range(1,cols-1):
            #print(row,col,varArr[row,col])
            # W coast
            if ( not(is_land(bathy,row,col)) and is_land(bathy,row,col+1) ):
                horArr=integrate_vector(varArr,horArr,bathy,row,col,'west')
            # E coast
            elif ( is_land(bathy,row,col) and not(is_land(bathy,row,col+1)) ):
                horArr=integrate_vector(varArr,horArr,bathy,row,col+1,'east')

                # Second passage uses boundary flux as reference
    #print('    Second passage using boundary as reference')
    # This is only used where there are no coastal reference as 
    # the error is larger
    for col in [1,cols-2]:
        for row in range(1,rows-1):
            #print('2nd- ',rows,cols,row,col)
            # Open Western boundary
            if (( col == 1 ) and  not(is_land(bathy,row,col))):
                #pass
                horArr=integrate_vector(varArr,horArr,bathy,row,col,'north')
            # Open Eastern boundary
            elif (( col == cols-2 ) and  not(is_land(bathy,row,col))):
                horArr=integrate_vector(varArr,horArr,bathy,row,col,'south')
                                        
    return horArr



In [9]:

    
def integrate_yadv(varArr,bathy):
    # for YAD
    [rows,cols]=varArr.shape
    #print('inegrate_yadv,rows,cols',rows,cols)
    verArr=np.zeros([rows,cols])
    
    # First passage uses coastal zero flux as reference
    # Does't visit last row or col due to the way it checks coastlines
    for row in range(1,rows-1):
        for col in range(1,cols-1):
            #print('1st- ',rows,cols,row,col)
            # N coast
            if ( not(is_land(bathy,row,col)) and is_land(bathy,row+1,col) ):
                verArr=integrate_vector(varArr,verArr,bathy,row,col,'south')
            # S coast
            elif ( is_land(bathy,row,col) and not(is_land(bathy,row+1,col)) ):
                verArr=integrate_vector(varArr,verArr,bathy,row+1,col,'north')
                
    #print('    Second passage using boundary as reference')
    
    # Second passage uses boundary flux as reference
    # This is only used where there are no coastal reference as 
    # the error is larger
    for row in [1,rows-2]:
        for col in range(1,cols-1):
            #print('2nd- ',rows,cols,row,col)
            # Open Southern boundary
            if (( row == 1 ) and  not(is_land(bathy,row,col))):
                #pass
                verArr=integrate_vector(varArr,verArr,bathy,row,col,'north')
            elif (( row == rows-2 ) and  not(is_land(bathy,row,col))):
                verArr=integrate_vector(varArr,verArr,bathy,row,col,'south')
                  
    return verArr

Call functions

Test at Northern boundary - Icelandic Coast



In [11]:

    
varDA=modelout[yadv_3d]
#varDA=modelout[test_2d]

# Test off Portuguese coast
varArr=varDA.values[0,0,367:,52:58].squeeze()
bath=bathy[367:,52:58]


print(varArr.shape)
verArr= integrate_yadv(varArr,bath)

# plot YAD
fig,axs = plt.subplots(1,3,figsize=(12,6),squeeze=False)

plt0 = axs[0,0].imshow(~(bath>0),origin='lower')
axs[0,0].set_title('land')
plt.colorbar(plt0,ax=axs[0,0])

plt1 = axs[0,1].imshow(varArr,origin='lower')
axs[0,1].set_title('YAD')
plt.colorbar(plt1, ax=axs[0,1])

plt2 = axs[0,2].imshow(verArr,origin='lower')
axs[0,2].set_title('adv tracer')
plt.colorbar(plt2,ax=axs[0,2])


#axs[0,0].colorbar()
#print(bath)
#print('YAD')
#print(varArr[:,:])
#print('y adv')
#print(verArr[:,:])









    



(8, 6)






    Out[11]:





<matplotlib.colorbar.Colorbar at 0xc2e8f60>

Test at Southern boundary - Portuguese Coast



In [14]:

    
varDA=modelout[yadv_3d]
#varDA=modelout[test_2d]

# Test off Portuguese coast
varArr=varDA.values[0,0,0:5,98:104].squeeze()
bath=bathy[0:5,98:104]


print(varArr.shape)
verArr= integrate_yadv(varArr,bath)

# plot YAD
fig,axs = plt.subplots(1,3,figsize=(12,6),squeeze=False)

plt0 = axs[0,0].imshow(~(bath>0),origin='lower')
axs[0,0].set_title('land')
plt.colorbar(plt0,ax=axs[0,0])

plt1 = axs[0,1].imshow(varArr,origin='lower')
axs[0,1].set_title('YAD')
plt.colorbar(plt1, ax=axs[0,1])

plt2 = axs[0,2].imshow(verArr,origin='lower')
axs[0,2].set_title('adv tracer')
plt.colorbar(plt2,ax=axs[0,2])


#axs[0,0].colorbar()
#print(bath)
#print('YAD')
#print(varArr[:,:])
#print('y adv')
#print(verArr[:,:])
# there is an error 1000x below value of fluxa









    



(5, 6)






    Out[14]:





<matplotlib.colorbar.Colorbar at 0xd13bb38>

Test for small embayment

Integrate XAD horizontally



In [15]:

    
varDA=modelout[xadv_3d]
#varDA=modelout[test_2d]

# test small embayment
varArr=varDA.values[0,0,208:212,149:154].squeeze()
bath=bathy[208:212,149:154]

# Test off Portuguese coast
#varArr=varDA.values[0,0,0:50,70:120].squeeze()

print(varArr.shape)
horArr = integrate_xadv(varArr,bath)

fig,axs = plt.subplots(1,3,figsize=(12,6),squeeze=False)

plt0 = axs[0,0].imshow(~(bath>0),origin='lower')
axs[0,0].set_title('land')
plt.colorbar(plt0,ax=axs[0,0])

plt1 = axs[0,1].imshow(varArr,origin='lower')
axs[0,1].set_title('XAD')
plt.colorbar(plt1, ax=axs[0,1])

plt2 = axs[0,2].imshow(horArr,origin='lower')
axs[0,2].set_title('adv tracer')
plt.colorbar(plt2,ax=axs[0,2])


#axs[0,0].colorbar()
#print(varArr[:,:])
#print(horArr[:,:])
# there is an error 1000x below value of flux









    



(4, 5)






    Out[15]:





<matplotlib.colorbar.Colorbar at 0xd2cd358>

Integrate YAD N-S



In [18]:

    
varDA=modelout[yadv_3d]
#varDA=modelout[test_2d]

# test small embayment
varArr=varDA.values[0,0,208:212,149:154].squeeze()
bath=bathy[208:212,149:154]

# Test off Portuguese coast
#varArr=varDA.values[0,0,0:50,70:120].squeeze()

print(varArr.shape)
verArr = integrate_yadv(varArr,bath)

# plot YAD
fig,axs = plt.subplots(1,3,figsize=(12,6),squeeze=False)

plt0 = axs[0,0].imshow(~(bath>0))
axs[0,0].set_title('land')
plt.colorbar(plt0,ax=axs[0,0])

plt1 = axs[0,1].imshow(varArr)
axs[0,1].set_title('YAD')
plt.colorbar(plt1, ax=axs[0,1])

plt2 = axs[0,2].imshow(verArr)
axs[0,2].set_title('adv tracer')
plt.colorbar(plt2,ax=axs[0,2])

#axs[0,0].colorbar()
#print(varArr[:,:])
#print(verArr[:,:])
# there is an error 1000x below value of flux









    



(4, 5)






    Out[18]:





<matplotlib.colorbar.Colorbar at 0xcd9acf8>

Calculate for whole image



In [19]:

    
varDA=modelout[xadv_3d]
#varDA=modelout[test_2d]
varArr=varDA.values[0,0,:,:].squeeze()
horArr = integrate_xadv(varArr,bathy)


varDA=modelout[yadv_3d]
varArr=varDA.values[0,0,:,:].squeeze()
verArr = integrate_yadv(varArr,bathy)

print(varArr.shape)

Plot maps



In [20]:

    
fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(12,12))

lim1=max(abs(np.array([verArr.max(),verArr.min(),horArr.max(),horArr.min()])))/5
lim2=max(abs(np.array([modelout[xadv_3d][0,0,:,:].max(),modelout[xadv_3d][0,0,:,:].min()])))/2

#axs[0,0].imshow(np.log10(bathy))
axs[0,0].imshow(np.log10(bathy),origin='lower')
axs[0,0].set_title('bathymetry from bathy_meter.nc')

im = axs[0,1].imshow(horArr,origin='lower',vmax=lim1,vmin=-1*lim1,cmap='seismic')
axs[0,1].set_title('abs. hor. advection')
fig.colorbar(im, ax=axs[0,1])

im = axs[1,1].imshow(modelout[xadv_3d][0,0,:,:],origin='lower',vmax=lim2,vmin=-1*lim2,cmap='seismic')
axs[1,1].set_title(xadv_3d)
fig.colorbar(im, ax=axs[1,1])

im = axs[0,2].imshow(verArr,origin='lower',vmax=lim1,vmin=-1*lim1,cmap='seismic')
axs[0,2].set_title('abs. ver. advection')
fig.colorbar(im, ax=axs[0,2])

im = axs[1,2].imshow(modelout[xadv_3d][0,0,:,:],origin='lower',vmax=lim2,vmin=-1*lim2,cmap='seismic')
axs[1,2].set_title(yadv_3d)
fig.colorbar(im, ax=axs[1,2])









    



C:\Users\tams00\AppData\Local\Continuum\anaconda3\lib\site-packages\ipykernel_launcher.py:7: RuntimeWarning: divide by zero encountered in log10
  import sys






    Out[20]:





<matplotlib.colorbar.Colorbar at 0xe5c5ac8>

Table of Contents

Functions

Call functions

Test at Northern boundary - Icelandic Coast

Test at Southern boundary - Portuguese Coast

Test for small embayment

Integrate XAD horizontally

Integrate YAD N-S

Calculate for whole image

Plot maps