Vertical tracer advection

[kate]

I'm having the model blow up after ten months or so of running. The nature of the trouble is isolated peaks in the surface temperature. Some grow then stabilize without making the model blow up, but eventually one grows big enough to trip the density check in diag. I assume that I could reduce the timestep to get through these rough spots, but the temperature peaks are still worrisome to the biologists counting on putting an NPZ model into this simulation. Using the debugger, I can see peaks growing during the vertical advection step. Some cppdefs:

SPLINES Conservative parabolic spline reconstruction.
TS_U3HADVECTION Third-order upstream horizontal advection of tracers.
TS_SVADVECTION Parabolic splines vertical advection of tracers.

You want plots, you say. Here is a small patch of SST(nnew) after the horizontal advection, before any vertical operation, followed by later images of the same patch while looping through j. I stopped it after the vertical advection and before the vertical diffusion, so as to see which was causing the growth:









The vertical advection is mostly a function of W, put into the FC array. Here are mid-depth W and FC corresponding to picture 6 above. The troublesome spot has a near-zero W with a high to one side and a low to the other (218):



This algorithm doesn't depend on a positive definite tracer does it? Does anyone have any suggestions?

[arango]

Why are you using TS_SVADVECTION This option must be only used in shallow, high resolution applications with a lot of vertical levels. It should not be used in coarse, basin scale applications with high stretched vertical grids. I have mentioned this several times in the past. Use TS_C4VADVECTION instead

It has also been reported in this forum that the SPLINES option may affect the downward flux in the vertical diffusion/viscosty terms. Many users get better and the correct results by not activating this option.

Finally, we are aware that TS_U3HADVECTION produces excessive numercal dyapicnal diffusion in highly stretched vertical grids. Jacopo posted a message about this recently. I prefer to use the other advection schemes available and add geopontential or isopycnic diffusion explicitly.

[kate]

Thanks, Hernan, we will try some other options then. I thought we checked our options with you three years ago for the last big run and haven't changed the core ones. Note that the entire patch I show is on the shallow Chukchi shelf with N=60.

[kate]

Just FYI: the cases without SPLINES blow up faster than those with SPLINES. We're keeping our SPLINES for now and trying various advection schemes.

[wilkin]

Kate,

I'm curious what vertical turbulent closure, air-sea flux, and subgridscale lateral mixing options you use in conjunction with these advection options.

John.

[kate]

Hi John,

I'm using a bulk_flux parameterization we got from NCAR:

BULK_FLUXES Surface bulk fluxes parametererization.
CCSM_FLUXES Surface bulk fluxes parametererization.
DIURNAL_SRFLUX Modulate shortwave radiation by the local diurnal cycle.
ICE_THERMO Include ice thermodynamics.
SOLAR_SOURCE Solar Radiation Source Term.


DIFF_GRID Horizontal diffusion coefficient scaled by grid size.
MIX_GEO_TS Mixing of tracers along geopotential surfaces.
MIX_S_UV Mixing of momentum along constant S-surfaces.
RDRG_GRID Spatially varying bottom drag.
SPONGE Enhanced horizontal mixing in the sponge areas.
TS_DIF2 Harmonic mixing of tracers.
UV_VIS2 Harmonic mixing of momentum.
VISC_GRID Horizontal viscosity coefficient scaled by grid size.


LMD_BKPP KPP bottom boundary layer mixing.
LMD_CONVEC LMD convective mixing due to shear instability.
LMD_MIXING Large/McWilliams/Doney interior mixing.
LMD_NONLOCAL LMD convective nonlocal transport.
LMD_RIMIX LMD diffusivity due to shear instability.
LMD_SHAPIRO Shapiro filtering boundary layer depth.
LMD_SKPP KPP surface boundary layer mixing.

The case without SPLINES blew up at mid-depth (N=20ish) in a shallow Aleutian pass.
The SPONGE ramps up the horizontal smoothing around the open boundaries, including the Chukchi shelf. The C4VADVECTION blew up much like the SVADVECTION.

[wilkin]

Since it sounds like you are trying many different combinations of options, I'd suggest stepping away from KPP and trying something else. Long ago I stopped using KPP in shallow water though I don't remember the precise reasons, but it wasn't instability so much as unreasonable surface temperatures stemming from too much vertical stratification (I think).

What's your rationale for a grid size dependent (if that's what it is) bottom drag: RDRG_GRID

[kate]

Since it sounds like you are trying many different combinations of options, I'd suggest stepping away from KPP and trying something else. Long ago I stopped using KPP in shallow water though I don't remember the precise reasons, but it wasn't instability so much as unreasonable surface temperatures stemming from too much vertical stratification (I think).

OK, we can try that.

What's your rationale for a grid size dependent (if that's what it is) bottom drag: RDRG_GRID

It's actually depth-dependent, for improving the tidal simulation. Not having Sasha's code for keeping quadratic bottom drag in bounds, I'm using linear bottom drag.

[arango]

I agree that you are using a very complex set-up and it is hard to diagnose where the problem is. I am not familiar with your CCSM_FLUXES algorithm but it can be a candidate. Bulk fluxes parameterizations are based on similarity theory and look at specified stability tables to determine the turbulent air-sea fluxes. In very extreme forcing conditions, these fluxes maybe unstable and violate model time-steping.

Your problems seems to be associated with surface forcing. Why do you need two different bulk fluxes parameterizations?

[kate]

I knew you wouldn't like that CCSM_FLUXES thing. I'm actually not using two bulk fluxes, just using the old tag to get ROMS to read all the atmospheric fields. Bill Large doesn't like the Fairall scheme because it was designed for the tropics and doesn't taper off properly in the face of gale-force winds. I did look at the fluxes from both side by side so I don't think there are any glaring errors there, but I can try the Fairall.

I'm not having much luck with GLS_MIXING. It blows up in one of three places at the bottom within just a few timesteps, even with half the timestep I was using.

It's beginning to look like a miracle that NEP4 was as good as it was. We know we can make this run by deepening hmin and/or by turning off the tides, but well, we're being greedy.

[ggerbi]

Kate--

What settings are you using for GLS? I have found k-epsilon to be less stable than the other versions. Did you try k-omega or gen?

[kate]

Err, it was nothing about the advection after all. I checked the averages files for fields not in the restart file and found the culprit - shortwave radiation. Somehow, our input radiation has some negative values which became transformed into large positive values up near the 24-hour night boundary. Clamping it to be positive before the DIURNAL_SRFLUX give a smooth field.

By the way, ana_srflux should have "cff=cff1" instead of "cff=cff1*2.0_r8*pi" for the all day case.

[wilkin]

By the way, ana_srflux should have "cff=cff1" instead of "cff=cff1*2.0_r8*pi" for the all day case.

Are you sure about this Kate? It's a long time ago that I wrote this, and I don't have my notes with me to double-check, but I think the normalization was conceived as the area under the curve:

max ( 0, cff1 + cff2 * cos(angle) )

from angle = 0 to 2*pi.

In the case that cff1 > cff2 the integral doesn't get tangled with the nasty max function because the integral of the cosine part is zero, and the integral of the constant is just cff1*2*pi.

So the area under the curve, and the factor, should be cff = cff1*2*pi.

I thought I checked that these integrals were continuous as one passed through the full range of cff1 and cff2, but I still could have made a mistake putting it in ROMS.

If the factor is wrong because the integration limits are not [0,2*pi], then I think the other part of the formulation would be wrong too.

[kate]

I'm positive that I get smooth fields without the factor of 2 pi. Here are three-hourly averages for one summer day (hope these aren't too big):

[wilkin]

If you extract a time series of the srflx as modulated by the DIURNAL_SRFLUX option, does it have the same integral over one day as the input netcdf data? If you're right and that code is wrong, then it's likely there are other errors in that code.

[wilkin]

Kate, you're right. In the more typical case |cff2| > |cff1| where we get an actual night-time, i.e. the gnarly formula:

cff=(cff1*ACOS(-cff1/cff2)+SQRT(cff2*cff2-cff1*cff1))/pi

the limit as cff2 approaches cff1 gives a normalization factor of cff1 because acos(-1) = pi. As |cff2| decreases further it disappears from the integral, and cff should stay as cff1.

I logged the trac ticket.

This thread is now in the wrong place in the forum, oh well.

[OcGaBy]

Hi everyone

I thought I had a problem similar to Kate's... however she already fix it and I can't!.


My model is a realistic simulation of the east coast of Africa.

I test the grid with analytical profiles uniformly in the horizontal, and surface forcing and it worked well. The problem arose when I added realistic initial conditions, and boundary forcing. I have looked at all the input files fields and they look smooth...

The model blows up because one single cell, or a group of a few neighbor cells get very different SSH values (i.e. -6 and +5m).

I changed form TS_SVADVECTION to TS_C4ADVECTION as Hernan suggested, and that helped me to extend my run from less than 1 day to 138 days... but then it blows up similarly.

I am using monthly climatology surface forcing, I checked and all my short wave radiation values are positive. I haven't change the 2*pi factor.

I will appreciate any suggestion.
Many thanks in advance.

GaBy

These are the CPPdefs I am using and some print screen figures.

ANA_BSFLUX Analytical kinematic bottom salinity flux.
ANA_BTFLUX Analytical kinematic bottom temperature flux.
ASSUMED_SHAPE Using assumed-shape arrays.
AVERAGES Writing out time-averaged fields.
AVERAGES_FLUXES Writing out time-averaged surface fluxes.
DIURNAL_SRFLUX Modulate shortwave radiation by the local diurnal cycle.
DJ_GRADPS Parabolic Splines density Jacobian (Shchepetkin, 2002).
DOUBLE_PRECISION Double precision arithmetic.
EAST_FSGRADIENT Eastern edge, free-surface, gradient condition.
EAST_M2NUDGING Eastern edge, 2D momentum, passive/active outflow/inflow.
EAST_M2RADIATION Eastern edge, 2D momentum, radiation condition.
EAST_M3NUDGING Eastern edge, 3D momentum, passive/active outflow/inflow.
EAST_M3RADIATION Eastern edge, 3D momentum, radiation condition.
EAST_TNUDGING Eastern edge, tracers, passive/active outflow/inflow.
EAST_TRADIATION Eastern edge, tracers, radiation condition.
EAST_VOLCONS Eastern edge, enforce mass conservation.
LMD_CONVEC LMD convective mixing due to shear instability.
LMD_MIXING Large/McWilliams/Doney interior mixing.
LMD_NONLOCAL LMD convective nonlocal transport.
LMD_RIMIX LMD diffusivity due to shear instability.
LMD_SKPP KPP surface boundary layer mixing.
MASKING Land/Sea masking.
MIX_S_TS Mixing of tracers along constant S-surfaces.
MIX_S_UV Mixing of momentum along constant S-surfaces.
MPI MPI distributed-memory configuration.
NONLINEAR Nonlinear Model.
NONLIN_EOS Nonlinear Equation of State for seawater.
NORTH_VOLCONS Northern edge, enforce mass conservation.
NORTH_FSGRADIENT Northern edge, free-surface, gradient condition.
NORTH_M2NUDGING Northern edge, 2D momentum, passive/active outflow/inflow.
NORTH_M2RADIATION Northern edge, 2D momentum, radiation condition.
NORTH_M3NUDGING Northern edge, 3D momentum, passive/active outflow/inflow.
NORTH_M3RADIATION Northern edge, 3D momentum, radiation condition.
NORTH_TNUDGING Northern edge, tracers, passive/active outflow/inflow.
NORTH_TRADIATION Northern edge, tracers, radiation condition.
OUT_DOUBLE Double precision output fields in NetCDF files.
POWER_LAW Power-law shape time-averaging barotropic filter.
PROFILE Time profiling activated .
!RST_SINGLE Double precision fields in restart NetCDF file.
SALINITY Using salinity.
SOLAR_SOURCE Solar Radiation Source Term.
SOLVE3D Solving 3D Primitive Equations.
SOUTH_FSGRADIENT Southern edge, free-surface, gradient condition.
SOUTH_M2NUDGING Southern edge, 2D momentum, passive/active outflow/inflow.
SOUTH_M2RADIATION Southern edge, 2D momentum, radiation condition.
SOUTH_M3NUDGING Southern edge, 3D momentum, passive/active outflow/inflow.
SOUTH_M3RADIATION Southern edge, 3D momentum, radiation condition.
SOUTH_TNUDGING Southern edge, tracers, passive/active outflow/inflow.
SOUTH_TRADIATION Southern edge, tracers, radiation condition.
SOUTH_VOLCONS Southern edge, enforce mass conservation.
SPLINES Conservative parabolic spline reconstruction.
SPONGE Enhanced horizontal mixing in the sponge areas.
TS_U3HADVECTION Third-order upstream horizontal advection of tracers.
TS_C4VADVECTION Fourth-order centered vertical advection of tracers.
TS_DIF2 Harmonic mixing of tracers.
UV_ADV Advection of momentum.
UV_COR Coriolis term.
UV_U3HADVECTION Third-order upstream horizontal advection of 3D momentum.
UV_SADVECTION Parabolic splines vertical advection of momentum.
UV_QDRAG Quadratic bottom stress.
UV_VIS2 Harmonic mixing of momentum.
VAR_RHO_2D Variable density barotropic mode.
WESTERN_WALL Wall boundary at Western edge.

[kate]

So it always blows up along an open boundary? Your boundary condition for the free surface is essentially a wall. Do you have values you can provide as boundary conditions to go along with the rest of your boundary fields?

With or without tides, my boundary conditions on each open edge are:

WEST_FSCHAPMAN Western edge, free-surface, Chapman condition.
WEST_M2FLATHER Western edge, 2D momentum, Flather condition.
WEST_M3NUDGING Western edge, 3D momentum, passive/active outflow/inflow.
WEST_M3RADIATION Western edge, 3D momentum, radiation condition.
WEST_TNUDGING Western edge, tracers, passive/active outflow/inflow.
WEST_TRADIATION Western edge, tracers, radiation condition.

[OcGaBy]

Hi Kate

Thanks for your answer.

Since I change the advection scheme to TS_C4ADVECTION the blow up tend to be close to the coast.

My boundary conditions file includes monthly values for zeta. Is my understanding that ROMS interpolate between the values on the input files for each time steps it needs right...?

Anyway I am going to try the open boundary conditions you recommend.

Greetings

GaBy

[bibi951]

Hello,Kate
I do tide simulation,and I think bottom friction really matters a lot.I use RDRG2 but it seems not good for my case,since my area is too large(the whole China Sea).
I notice that you use spatial distributed bottom drag.

What's your rationale for a grid size dependent (if that's what it is) bottom drag: RDRG_GRID
It's actually depth-dependent, for improving the tidal simulation. Not having Sasha's code for keeping quadratic bottom drag in bounds, I'm using linear bottom drag.

and you say ,"My current implementation uses RDRG_GRID and expects to read a 2-D field from the grid file - a more complete version would also have an ANA_RDRG option and ana_rdrg to back it up." https://www.myroms.org/blog/
Is it available now in Roms? Or can you release a copy or something if it is open-source.

Thanks.

[kate]

It's all on my svn branch, available by emailing robertson@marine.rutgers.edu. User beware - it's my working code, not trying to be all that tidy. I ended up asking Sasha about the limiter and I now have DRAG_LIMITER as well. I'm still using linear bottom drag. Does the China Sea have large tides? I had to go to these lengths to tame some unphysical tides in Bristol Bay in the Bering Sea.