39.9.11 Discrete density-meridional streamfunction

The streamfunction computed with potential density as the vertical coordinate requires a bit more work. The most important part of the process is how to partition the potential density layers. In general, how the layers are placed will largely determine the quality of the diagnosed overturning streamfunction. Currently, the density range is divided into km layers, since this is the maximum number of horizontally uniform density layers which can be resolved in a level model. Using fewer density layers will generally produce a smoother streamfunction, but the value of the streamfunction is typically reduced in magnitude. More layers typically produces a noisy solution, and will make the streamfunction have larger values for its extrema. More sophisticated schemes are possible, in which various forms of interpolation are applied. These have not been tried. If someone has a better scheme, please feel free to suggest its use. Currently, the simplest algorithm was chosen, in which zero interpolation is used. Furthermore, since each model solution likely will require its own particular partitioning, it is currently left to the researcher to determine the density layers. The placement of the density layers occurs inside the routine diag.F. Look for the USER INPUT section within the ifdef option meridional_overturning.

Given a partitioning of the density coordinate $\sigma_{m}$, for m=1,km, the algorithm for computing the streamfunction is the following:

$$sigmsf_{jrow,m} = -\csuj \; \sum_{i=2}^{imt-1} \sum_{k=1}^{kmt(i,... ...)} \dxui \; dzt_{m} \; {\cal H}(\sigma_{i,k,j}-\sigma_{m}) \, u_{i,m,j,2,\tau},$$ (39.157)

where ${\cal H}(x)$ is a Heaviside step function (equal to zero for x<0 and unity for x>0.

As coded in MOM, there are two potential densities which are used for defining the streamfunction: potential density referenced to the surface, and potential density referenced to 2000m. The framework for using these potential densities can easily be extended.