

In [ ]:

    
import numpy as np
import numpy.ma as ma
from scipy.io import loadmat
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap


def curl(lon, lat, taux, tauy):
    """Return midlon, midlat, curl(tau) based on first differences.
    CAVEAT: Lon are columns and Lat are rows.  The differences might change if this is reverse!"""

    re = 6371e3
    latrad = np.deg2rad(lat)
    lonrad = np.deg2rad(lon)
    midlon = 0.5 * (lon[:-1] + lon[1:])
    midlat = 0.5 * (lat[:-1] + lat[1:])

    # See Vallis page 59.
    rfac = (1.0 / re) / np.cos(np.deg2rad(midlat))

    dvdx = np.diff(tauy, axis=-1) / np.diff(lonrad)
    dvdx = 0.5 * (dvdx[1:, ...] + dvdx[:-1, ...])

    du = np.diff(taux * np.cos(latrad)[..., None], axis=0)
    dudy = du / np.diff(latrad)[..., None]
    dudy = 0.5 * (dudy[..., 1:] + dudy[..., :-1])

    curltau = rfac[..., None] * (dvdx - dudy)
    return midlon, midlat, curltau


def make_basemap(llcrnrlat, urcrnrlat, llcrnrlon, urcrnrlon,
                 projection='merc', figsize=(8, 4), resolution='c'):
    fig, ax = plt.subplots(figsize=figsize)
    m = Basemap(projection='merc', llcrnrlat=llcrnrlat, urcrnrlat=urcrnrlat,
                llcrnrlon=llcrnrlon, urcrnrlon=urcrnrlon, resolution=resolution)
    m.ax = ax
    m.drawcoastlines()
    m.fillcontinents(color='0.85')
    dx, dy = 15, 10
    meridians = np.arange(llcrnrlon+dx, urcrnrlon+dx, dx)
    parallels = np.arange(llcrnrlat, urcrnrlat, dy)
    kw = dict(color='k', alpha=0.5)
    m.drawparallels(parallels, labels=[1, 0, 0, 0], **kw)
    m.drawmeridians(meridians, labels=[0, 0, 0, 1], **kw)
    return fig, m

def unit_uv(u, v):
    vec = u + 1j * v
    spd = np.abs(vec)
    ang = np.angle(vec, deg=True)
    ang = np.mod(90 - ang, 360)
    u = np.sin(np.deg2rad(ang))
    v = np.cos(np.deg2rad(ang))
    return u, v

Dados.



In [ ]:

    
AS_vento = loadmat('AS_vento.mat', squeeze_me=True)
lon, lat = AS_vento['lon'], AS_vento['lat']
lat = ma.masked_inside(lat, -4, 4)
u, v = AS_vento['uver'], AS_vento['vver']  # Verão.
#u, v = AS_vento['uinv'], AS_vento['vinv']  # Inverno.
u, v = ma.masked_invalid(u), ma.masked_invalid(v)
lon, lat = np.meshgrid(lon, lat)

rho_ar = 1.226  # kg m^{-3}
rho0 = 1025  # kg m^{-3}
omega = 2 * np.pi / (60*60*24)  # rad s^{-1}
f0 = 2 * omega * np.sin(np.deg2rad(lat))  # rad s^{-1}



In [ ]:

    
cut = slice(0, -1, 15), slice(0, -1, 15)

fig, m = make_basemap(llcrnrlat=lat.min(), urcrnrlat=lat.max(),
                       llcrnrlon=lon.min(), urcrnrlon=lon.max(),
                       projection='merc', figsize=(8, 8), resolution='i')
kw = dict(latlon=True, color='cornflowerblue')
Q = m.quiver(lon[cut], lat[cut], u[cut], v[cut], **kw)
arrow = 10
m.ax.quiverkey(Q, 0.2, 0.65, arrow, r'%s m s$^{-1}$' % arrow, color='k', labelpos='N', coordinates='figure')

Tensão de cisalhamento do vento.



In [ ]:

    
def drag_coeff(Ustar):
    """Mellor (2004)."""
    return 7.5e-4 + 6.7e-5 * Ustar



In [ ]:

    
mag = np.sqrt(u**2 + v**2)  # m s^{-1}
Cd = drag_coeff(mag)
taux = rho_ar * Cd * mag * u  # N m^{-2}
tauy = rho_ar * Cd * mag * v  # N m^{-2}



In [ ]:

    
cut = slice(0, -1, 15), slice(0, -1, 15)

fig, m = make_basemap(llcrnrlat=lat.min(), urcrnrlat=lat.max(),
                      llcrnrlon=lon.min(), urcrnrlon=lon.max(),
                      projection='merc', figsize=(8, 8), resolution='i')
kw = dict(latlon=True, color='cornflowerblue')
Q = m.quiver(lon[cut], lat[cut], taux[cut], tauy[cut], **kw)
arrow = 0.1
m.ax.quiverkey(Q, 0.2, 0.65, arrow, r'%s N m$^{-2}$' % arrow, color='k', labelpos='N', coordinates='figure')

Transporte de Ekman.

\begin{align*} 0 &= -fu + \dfrac{1}{\rho}\dfrac{\partial{\tau_y}}{\partial z} \\ 0 &= +fv + \dfrac{1}{\rho}\dfrac{\partial{\tau_x}}{\partial z} \\ \end{align*}



In [ ]:

    
Ue = +tauy / (rho0 * f0)  # m^2 s^{-1}
Ve = -taux / (rho0 * f0)  # m^2 s^{-1}



In [ ]:

    
cut = slice(0, -1, 15), slice(0, -1, 15)

fig, m = make_basemap(llcrnrlat=lat.min(), urcrnrlat=lat.max(),
                       llcrnrlon=lon.min(), urcrnrlon=lon.max(),
                       projection='merc', figsize=(8, 8), resolution='i')
kw = dict(latlon=True, color='cornflowerblue')
Q = m.quiver(lon[cut], lat[cut], Ue[cut], Ve[cut], **kw)
arrow = 2
m.ax.quiverkey(Q, 0.2, 0.65, arrow, r'%s m$^2$ s$^{-1}$' % arrow, color='k', labelpos='N', coordinates='figure')

Componente vertical do rotacional do vetor tensão de cisalhamento do vento.

Bombeamento de Ekman.



In [ ]:

    
midlon, midlat, rot = curl(lon[0, :], lat[:, 0], taux, tauy)
midlon, midlat = np.meshgrid(midlon, midlat)



In [ ]:

    
f0c = 2 * omega * np.sin(np.deg2rad(midlat))
We = rot / (rho0 * f0c)  # m s^{-1}



In [ ]:

    
cut = slice(0, -1, 15), slice(0, -1, 15)

fig, m = make_basemap(llcrnrlat=lat.min(), urcrnrlat=lat.max(),
                      llcrnrlon=lon.min(), urcrnrlon=lon.max(),
                      projection='merc', figsize=(8, 8), resolution='i')
scale = 1e5
cs = m.pcolormesh(midlon, midlat, We*scale, cmap=plt.cm.RdBu, vmin=-0.8, vmax=0.8, latlon=True)
cbar = fig.colorbar(cs, extend='both', orientation='vertical', shrink=0.65, pad=0.05)
cbar.ax.set_ylabel(r'[m s$^{-1}$] We $\times 10^-5$', rotation=-90, labelpad=20)

tx, ty = unit_uv(taux, tauy)
U, V = unit_uv(Ue, Ve)
m.quiver(lon[cut], lat[cut], tx[cut], ty[cut], latlon=True, color='green', scale=30)
m.quiver(lon[cut], lat[cut], U[cut], V[cut], latlon=True, color='black', scale=30)



In [ ]:

    
fig, ax = plt.subplots()
ax.plot(We.mean(axis=1)*scale, midlat[:, 0])
ax.vlines(0, midlat.min(), midlat.max())



In [ ]: