Blow up with LMD_MIXING and GLS_MIXING but BVF_MIXING works

[aryansafaie]

Dear ROMS Community,

I'm running an idealized simulation of an along-shore uniform, wedge-shaped (simple linear sloped bottom) coastal region, to hopefully observe a thermally-induced exchange flow. I'm aiming for a 4 week simulation time, with a time step of 5 seconds. My primary forcing is an analytical expression for shortwave radiation using ana_srflux.h, and also an expression in ana_stflux.h that removes heat from the surface.

The issue I'm having is that within my CPP options, if I use BVF_MIXING, my simulation runs to completeness and behaves as expected. However, whenever I instead activate LMD_MIXING or GLS_MIXING and ANY of their sub options (all other things being equal), my simulation blows up, usually between 17-22 hours of simulation time.

Before you proceed to read the following details, have any of you ever experienced something similar? I should say that I'm still a novice ROMS user, and it's likely I'm making a rookie mistake.

I've attached the model output for a simulation that blew up after 22 hours of simulation time using GLS mixing (with N2S2_HORAVG and the CANUTO_A stability function also activated) with the gen parameters as described in Warner et al. (2005);
I've also attached the CPP options I used (also copied and pasted below), as well as the .in file.

Please let me know if there's any additional info I can provide to help explain this problem.

Thank you for your time and assistance!
-Aryan

/*
** Application flag: FARROW
** Input script: ocean_farrow.in
*/

/* Error: ocean_frc.nc */

/* Aryan Safaie 06/22/15
Only forcing will be from srflux, and for this we will use an analytical expression, not
bulk_flux.F; therefore we can undefine most everything here
*/

/* First define necessecities for double exponential decay function */
#define NONLINEAR /* Required for nonlinear analytical and active fluxes */
#define SHORTWAVE /* Required for ana_stflux accessing srflux */
#define SOLAR_SOURCE /* Required for lmd_swfrac (double exp decay) */
#define ANA_SRFLUX /* Analytical shortwave radiation flux */
#define ANA_STFLUX /* Analytical kinematic surface temperature flux */


#undef BULK_FLUXES


/* Temperature */
#define ANA_BTFLUX /* Analytical kinematic bottom temperature flux */
#define ANA_SST /* Sea surface temperature SST */

/* Salinity (files not found in Functionals) */
#undef ANA_SSFLUX /* Analytical kinematic surface salinity flux (missing file)*/
#undef ANA_BSFLUX /* Analytical kinematic bottom salinity flux (missing file) */
#define ANA_SSS /* Sea surface salinity SST */

/* Momentum fluxes */
#define ANA_SMFLUX /* Analytical kinematic surface momentum flux */
#undef ANA_BMFLUX /* Analytical kinematic bottom momentum flux (missing file) */

#undef ANA_VMIX /**/


/* Not defining BULK_FLUXES, so the following should not be implemented,
though make sure to double check this
*/
/* Need to be defined in order to avoid calling FRC_FILE
- LONGWAVE *
- PAIR
- TAIR
- HUMIDITY (HUMID) *
- CLOUD *
- RAIN *
- WINDS *
- SRFLUX
*/
#ifdef BULK_FLUXES
#define ANA_WINDS
#define ANA_TAIR
#define ANA_PAIR
#define ANA_HUMIDITY
#define ANA_RAIN
#define LONGWAVE
#define ANA_CLOUD
#endif

/* Undefined analytical files */
#undef ANA_GRID
#define ANA_INITIAL

/* Do these fix the blow up after 66 hours problem? */
#undef NONLIN_EOS /* use if using nonlinear equation of state */
#undef SALINITY

/*
**
*/
#define AVERAGES
#define DIAGNOSTICS_TS
#define DIAGNOSTICS_UV
#define DJ_GRADPS
#define MIX_S_TS
#define MIX_S_UV
#define OUT_DOUBLE
#define SOLVE3D
#define SPLINES
#define TS_U3HADVECTION
#define TS_C4VADVECTION
#define TS_DIF2
#define UV_ADV
#define UV_CORR
#define UV_LDRAG

/* Now choose a Mixing Scheme: GLS, LMD, BVF, or MY25 */

/* Activate GLS mixing (Warner et al./ 2005) */
#define GLS_MIXING /*generic length scale mixing*/
#ifdef GLS_MIXING
# define N2S2_HORAVG /*horiz smoothing of buoyancy/shear*/
# define CANUTO_A /*Canuto-A stability function*/
# undef CRAIG_BANNER /*wave breaking surface flux*/
# undef CHARNOK /*Charnok surface roughness from wind stress*/
# undef KANTHA_CLAYSON /*KC stability function */
#endif

/* Activate Large et al. (1994) KPP Vertical Mixing Scheme */
#undef LMD_MIXING /* Large et al 1994 mixing */
#ifdef LMD_MIXING
# define LMD_RIMIX /* diffusivity due to shear instability */
# define LMD_CONVEC /* convective mixing due to shear instability */
# define LMD_SKPP /* surface boundary layer KPP mixing */
# define LMD_NONLOCAL /*nonlocal transport for active (TS) tracers */
#endif

/* Use Buoyancy Frequency Mixing */
#undef BVF_MIXING

#undef MIX_GEO_TS
#undef MIX_GEO_UV
#undef TS_DIF4
#undef TS_MPDATA

[kate]

I have had ROMS blow up many a time, for all sorts of reasons. I suggest you save a snapshot right as it blows up - perhaps you have one in the restart file already. Investigate the file with say ncview. You can also try a shorter timestep to see if you were right on the edge of stability with the other mixing scheme.

[aryansafaie]

Hi again,

I've included an image of the surface temperature and the mid-depth salinity from the restart file of another similar run. Nearshore (the left side) is full of NaNs.


I've also included a transect of the temperature field through time. Just before the blow up, the nearshore region cools at a noticeably greater rate.
I'm not entirely sure why this blow up occurs, nor why the simulation works if I use the BVF mixing scheme instead of the current GLS (gen) scheme.

The run log is also attached below, as well as the .in file and the CPP defs.

Thanks again for your help!
Best,
-aryan

[kate]

You've got a dense bottom current flowing downslope from the left to the right. Dale Haidvogel has studied a similar problem and many models have various troubles with such a thing. You might find that some GLS options are more stable than others.

I see you have zero depth at the coast. Try deepening that just a bit.

[m.hadfield]

I see you have zero depth at the coast. Try deepening that just a bit.

A note for novices: without wetting and drying enabled, ROMS will blow up if zeta drops below -h at any grid point.

[aryansafaie]

Thank you both for your quick replies.

"#define WET_DRY" ---> wow, silly me, I could have sworn I already activated this. Thank you so much for pointing that out!

I'm investigating GLS parameters that will hopefully bring the model to stability; I'll post when I find a set that works.

With regards to the minimum depth, this grid is set to a minimum depth of 4 m, which may not be deep enough, and in the creation of the .nc grid file, masking was indeed turned on such that rho points shallower than the min depth were masked with 0.
I do have some confusion regarding this however:
When I create a grid (as an .nc file), I can set a minimum depth, and I can activate masking. Concurrently, in the ocean.in file, there is a DCRIT parameter (right now I'm setting it to 4 m) for wetting/drying. As some cells will already be masked with 0 (and will always be dry, with a depth = DCRIT), it seems I need to increase DCRIT to a value even greater than 4 m? However, as I haven't defined WET_DRY yet, this seems like a moot point.

I'll report back. Thank you both again!

Best,
-aryan

[kate]

You can't have a depth of zero, even inside the land mask, unless you have WET_DRY. You shouldn't have to go as deep as 4 m for your DCRIT. I got by with 0.1 m in a Cook Inlet domain.

[aryansafaie]

Quick update:
with the WET_DRY CPP option defined, and a minimum depth of 10 m (which I could probably reduce), the GLS simulation has lasted longer than ever: 8 days 20 hours and counting!
Thank you both for your assistance.

When this is all said and done, I'm still very curious as to why the BVF mixing scheme did not blow up, even with my forgetful omission of "#define WET_DRY". I'm reading about the GLS mixing scheme from Warner et al.,(2005), and soon Umlauf and Burchard (2003), hopefully I'll have an idea soon.

Thanks again,
-aryan

[aryansafaie]

Alright, so the simulation ran to completion after adding #define WET_DRY to the CPP defs.

Thank you both very much for your help.