s: streamfunction
def rhs(nx,ny,dx,dy,re,w,s):
aa = 1.0/(re*dx*dx)
bb = 1.0/(re*dy*dy)
gg = 1.0/(4.0*dx*dy)
hh = 1.0/3.0
f = np.empty((nx+3,ny+3))
for i in range(1,nx+2):
for j in range(1,ny+2):
#Arakawa
j1 = gg*( (w[i+1,j]-w[i-1,j])*(s[i,j+1]-s[i,j-1]) \
j2 = gg*( w[i+1,j]*(s[i+1,j+1]-s[i+1,j-1]) \
+ w[i,j-1]*(s[i+1,j-1]-s[i-1,j-1]))
j3 = gg*( w[i+1,j+1]*(s[i,j+1]-s[i+1,j]) \
+ w[i+1,j-1]*(s[i+1,j]-s[i,j-1]) )
jac = (j1+j2+j3)*hh
#Central difference for Laplacian
f[i,j] = -jac + aa*(w[i+1,j]-2.0*w[i,j]+w[i-1,j]) \
+ bb*(w[i,j+1]-2.0*w[i,j]+w[i,j-1])
```
✝ ✆ **Listing 31.** Computation of right hand side function in Python with vectorization. ✞ ☎

```