Closed HenryDane closed 1 year ago
Code to reproduce:
import climlab
from climlab import constants as const
import numpy as np
import matplotlib.pyplot as plt
full_state = climlab.surface_state(num_lat=40, num_lon=30, water_depth=2.5)
model = climlab.TimeDependentProcess(state=full_state, timestep=const.seconds_per_hour)
medif = climlab.dynamics.meridional_heat_diffusion.MeridionalHeatDiffusion(state=full_state)
zodif = climlab.dynamics.zonal_heat_diffusion.ZonalHeatDiffusion(state=full_state)
model.add_subprocess('medif', medif)
model.add_subprocess('zodif', zodif)
for i in range(8765):
model.Ts[:,0,:] = 100.0
model.step_forward()
plt.imshow(model.Ts)
plt.imshow(medif.heat_transport)
plt.xlabel('Lon (indices)')
plt.ylabel('Lat (indices)')
plt.title('Meridional heat transport')
plt.colorbar()
plt.imshow(zodif.heat_transport)
plt.xlabel('Lon (indices)')
plt.ylabel('Lat (indices)')
plt.title('Zonal heat transport')
plt.colorbar()
plt.plot(model.Ts[0,:], label='lon index=0')
plt.plot(model.Ts[-1,:], label='lon index=-1')
plt.legend()
plt.gca().ticklabel_format(useOffset=False, style='plain')
plt.xlabel('Longitude index')
plt.ylabel('Temperature')
I think that the bottom left and top right corners of the tridiagonal matrix need to be set to non-zero values, which implies that this wouldn't be tri-diagonal anymore...
Specifically, the top-right should look like $T{l0}$ and the bottom-left should look like $T{u(J-1)}$.
http://www.physics.iisc.ac.in/~prateek/numerical_analysis/hw2.pdf has a relevant example on page 2.
This was fixed in [c602133](https://github.com/HenryDane/climlab/commit/c6021336f444cfca1b998600aacc751dc77235c9)
.
Diffusion/advection does not "wrap around."