This code solves time-dependent diffusion equations for zonally averaged QGPV in a two-level atmosphere.
The resulting surface stress curl is used to compute heat transport in an underlying ocean layer, assuming that transport is dominated by wind-driven gyres.

Reference:

Rose and Marshall, 2009: Ocean heat transport, sea ice, and multiple climate states: insights from energy balance models, J. Atmos. Sci. 66, 2828--2843

For more details, see also Rose, B.E.J., 2010: Oceanic Control of the Sea Ice Edge and Multiple Equilibria in the Climate System, PhD thesis, Massachusetts Institute of Technology.

Code is provided as-is. Proceed with caution.
I am not actively developing or maintaining the code at this point, but I hope it will be useful to others.
Please feel free to raise issues here on github.

Three versions of the code are provided:

• time-stepping codes for Fortran 95
• time-stepping code for MATLAB,
• MATLAB code for the steady-state boundary value problem

Fortran code was used to produce the calculations in Rose and Marshall (2009). MATLAB code is slower but should do the same thing. However this code is older and there are no guarantees that it is free of bugs.

• emombm.f95: the main model code. Notation generally follows what's laid out in Rose and Marshall (2009), so should be understandable. Note that these are just subroutines that need to be invoked by a main program.
• run_emombm.f95: an example program that invokes the model and integrates it for a specified number of timesteps.
• iceline_emombm.f95: a program that iterates over an array of parameters, in this case solar constant, timestepping the model out to equilibrium for each parameter value.
• tridiag.f95: a useful subroutine for solving a tridiagonal linear system, which arises from the implicit timestepping scheme on the diffusion operator.
• SIZE.h: a header file that specifies the grid size.
• paramfile: contains parameters that are read in at run-time. Not everything in this file is used in these examples.
• extract_results.f95: a standalone program that I used to read particular results from a large output file, such as generated by iceline_emombm.f95

I was able to compile and run this code with gfortran on several different machines.

• computeku.m
• current_diagnostics.m
• EMomBM.m
• fastSolveU.m
• final_diagnostics.m
• fluxes.m
• heating.m
• iceline_qg.m
• Kvectors.m
• marshallY.m
• noflow.m
• pvflux.m
• stress.m
• temperature.m

Try just running iceline_qg at the MATLAB prompt.

• run_EMomBM.m
• YetAnotherEMomBM.m

Type run_EMomBM at the MATLAB prompt.

This is offered as an example of how to set up the boundary value problem for surface wind stress, using the built-in MATLAB BVP solvers. There is no ocean component in this version, but one could be added by adding another ODE and appropriate boundary conditions. For details on the steady-state problem please see Chapter 3 of my thesis cited above.