Carl-
The upwelling code is define to use a constant value of vertical mixing, as you can tell from this code in upwelling.h
#if defined GLS_MIXING || defined MY25_MIXING
# define KANTHA_CLAYSON
# define N2S2_HORAVG
#else
# define ANA_VMIX
#endif
and by noting that neither GLS_MIXING or MY25_MIXING are enabled earlier in the file. Furthermore, on the side of the "upwelling" model case in which there is downwelling, the left side in your plot, light water is being forced under denser water. This leads to a static instability in the water column. With a low and constant viscosity, this leads to convection in the model.
Unfortunately, in a relatively coarsely resolved hydrostatic model, this convection is poorly represented as 2 Delta X upwelling cells, as you see in your plot. To model the convection properly, you would need a horizontal resolution that is small compared to the water depth (!!) and a non-hydrostatic model. This will tend to be computationally unfeasible if you want to use the model to resolve any reasonably large area. Thus you must parameterize convection. MY25_MIXING will do an ok-ish job of doing this by cranking up the eddy viscosity. See Austin & Lentz's paper on the inner shelf for a discription of how the inner shelf will develop in a case like that shown in the "upwelling" model case.
P.s. how well developed is your python plotting code for ROMS output -- is it good enough to share?
Cheers,
Jamie Pringle
University of New Hampshire