First of all, sorry for the lengthy post. I tried to make it as short as possible without disregarding important info.
I know that a lot of posts are about ROMS blowing up, and I also know that the solution is very user-specific. I do not expect someone to solve my problem, but maybe give me clues on how to identify the reason why my model is blowing up at the second timestep.
My setup is a channel (2 closed boundary, one clamped boundary upstream, one radiation boundary downstream), with a clamped 2D flow at both the upstream and downstream boundary. Here is what ROMS prints (I striped off what I think is non relevant, just to make the post less overwhelming)
Code: Select all
Resolution, Grid 01: 0050x0200x030, Parallel Threads: 1, Tiling: 001x001
Physical Parameters, Grid: 01
=============================
86400 ntimes Number of timesteps for 3-D equations.
231.000 dt Timestep size (s) for 3-D equations.
36 ndtfast Number of timesteps for 2-D equations between
each 3D timestep.
1 ERstr Starting ensemble/perturbation run number.
1 ERend Ending ensemble/perturbation run number.
0 nrrec Number of restart records to read from disk.
T LcycleRST Switch to recycle time-records in restart file.
200 nRST Number of timesteps between the writing of data
into restart fields.
1 ninfo Number of timesteps between print of information
to standard output.
T ldefout Switch to create a new output NetCDF file(s).
200 nHIS Number of timesteps between the writing fields
into history file.
1 ntsAVG Starting timestep for the accumulation of output
time-averaged data.
200 nAVG Number of timesteps between the writing of
time-averaged data into averages file.
1 ntsDIA Starting timestep for the accumulation of output
time-averaged diagnostics data.
200 nDIA Number of timesteps between the writing of
time-averaged data into diagnostics file.
0.0000E+00 nl_tnu2(01) NLM Horizontal, harmonic mixing coefficient
(m2/s) for tracer 01: temp
0.0000E+00 nl_tnu2(02) NLM Horizontal, harmonic mixing coefficient
(m2/s) for tracer 02: salt
0.0000E+00 nl_visc2 NLM Horizontal, harmonic mixing coefficient
(m2/s) for momentum.
5.0000E-06 Akt_bak(01) Background vertical mixing coefficient (m2/s)
for tracer 01: temp
5.0000E-06 Akt_bak(02) Background vertical mixing coefficient (m2/s)
for tracer 02: salt
5.0000E-06 Akv_bak Background vertical mixing coefficient (m2/s)
for momentum.
5.0000E-06 Akk_bak Background vertical mixing coefficient (m2/s)
for turbulent energy.
5.0000E-06 Akp_bak Background vertical mixing coefficient (m2/s)
for turbulent generic statistical field.
3.0000E-04 rdrg Linear bottom drag coefficient (m/s).
3.0000E-03 rdrg2 Quadratic bottom drag coefficient.
2.0000E-02 Zob Bottom roughness (m).
2 Vtransform S-coordinate transformation equation.
4 Vstretching S-coordinate stretching function.
7.0000E+00 theta_s S-coordinate surface control parameter.
2.5000E+00 theta_b S-coordinate bottom control parameter.
30.000 Tcline S-coordinate surface/bottom layer width (m) used
in vertical coordinate stretching.
1025.000 rho0 Mean density (kg/m3) for Boussinesq approximation.
0.000 dstart Time-stamp assigned to model initialization (days).
0.00 time_ref Reference time for units attribute (yyyymmdd.dd)
0.0000E+00 Tnudg(01) Nudging/relaxation time scale (days)
for tracer 01: temp
0.0000E+00 Tnudg(02) Nudging/relaxation time scale (days)
for tracer 02: salt
0.0000E+00 Znudg Nudging/relaxation time scale (days)
for free-surface.
0.0000E+00 M2nudg Nudging/relaxation time scale (days)
for 2D momentum.
0.0000E+00 M3nudg Nudging/relaxation time scale (days)
for 3D momentum.
0.0000E+00 obcfac Factor between passive and active
open boundary conditions.
F VolCons(1) NLM western edge boundary volume conservation.
F VolCons(2) NLM southern edge boundary volume conservation.
F VolCons(3) NLM eastern edge boundary volume conservation.
F VolCons(4) NLM northern edge boundary volume conservation.
9.000 T0 Background potential temperature (C) constant.
32.000 S0 Background salinity (PSU) constant.
1025.000 R0 Background density (kg/m3) used in linear Equation
of State.
1.7000E-04 Tcoef Thermal expansion coefficient (1/Celsius).
7.6000E-04 Scoef Saline contraction coefficient (1/PSU).
-1.000 gamma2 Slipperiness variable: free-slip (1.0) or
no-slip (-1.0).
...
Output/Input Files:
Output Restart File: ocean_rst.nc
Output History File: ocean_his.nc
Output Averages File: ocean_avg.nc
Output Diagnostics File: ocean_dia.nc
Tile partition information for Grid 01: 0050x0200x0030 tiling: 001x001
tile Istr Iend Jstr Jend Npts
0 1 50 1 200 300000
Tile minimum and maximum fractional grid coordinates:
(interior points only)
tile Xmin Xmax Ymin Ymax grid
0 0.50 51.50 0.50 201.50 RHO-points
0 0.00 51.00 0.50 201.50 U-points
0 0.50 51.50 0.00 201.00 V-points
Lateral Boundary Conditions: NLM
============================
Variable Grid West Edge South Edge East Edge North Edge
--------- ---- ---------- ---------- ---------- ----------
zeta 1 Closed Radiation Closed Clamped
ubar 1 Closed Radiation Closed Clamped
vbar 1 Closed Radiation Closed Clamped
u 1 Closed Radiation Closed Clamped
v 1 Closed Radiation Closed Clamped
temp 1 Closed Radiation Closed Clamped
salt 1 Closed Radiation Closed Clamped
tke 1 Closed Radiation Closed Clamped
Activated C-preprocessing Options:
COASTAL_CURRENT Coastal Current project
ANA_BSFLUX Analytical kinematic bottom salinity flux.
ANA_BTFLUX Analytical kinematic bottom temperature flux.
ANA_FSOBC Analytical free-surface boundary conditions.
ANA_GRID Analytical grid set-up.
ANA_INITIAL Analytical initial conditions.
ANA_M2OBC Analytical 2D momentum boundary conditions.
ANA_M3OBC Analytical 3D momentum boundary conditions.
ANA_SMFLUX Analytical kinematic surface momentum flux.
ANA_SSFLUX Analytical kinematic surface salinity flux.
ANA_STFLUX Analytical kinematic surface temperature flux.
ANA_TOBC Analytical tracers boundary conditions.
ASSUMED_SHAPE Using assumed-shape arrays.
AVERAGES Writing out time-averaged nonlinear model fields.
DIAGNOSTICS_TS Computing and writing tracer diagnostic terms.
DIAGNOSTICS_UV Computing and writing momentum diagnostic terms.
DJ_GRADPS Parabolic Splines density Jacobian (Shchepetkin, 2002).
DOUBLE_PRECISION Double precision arithmetic.
KANTHA_CLAYSON Kantha and Clayson stability function formulation.
MIX_S_TS Mixing of tracers along constant S-surfaces.
MIX_S_UV Mixing of momentum along constant S-surfaces.
MY25_MIXING Mellor/Yamada Level-2.5 mixing closure.
NONLINEAR Nonlinear Model.
!NONLIN_EOS Linear Equation of State for seawater.
POWER_LAW Power-law shape time-averaging barotropic filter.
PROFILE Time profiling activated .
K_GSCHEME Third-order upstream advection of TKE fields.
RADIATION_2D Use tangential phase speed in radiation conditions.
!RST_SINGLE Double precision fields in restart NetCDF file.
SALINITY Using salinity.
SOLVE3D Solving 3D Primitive Equations.
SPLINES Conservative parabolic spline reconstruction.
!STOCH_OPT_WHITE Stochastic Optimals, red noise processes.
TS_A4HADVECTION Fourth-order Akima horizontal advection of tracers.
TS_A4VADVECTION Fourth-order Akima 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_C4VADVECTION Fourth-order centered vertical advection of momentum.
UV_LDRAG Linear bottom stress.
UV_VIS2 Harmonic mixing of momentum.
VAR_RHO_2D Variable density barotropic mode.
Process Information:
Thread # 0 (pid= 7384) is active.
INITIAL: Configuring and initializing forward nonlinear model ...
Vertical S-coordinate System:
...
ndtfast, nfast = 36 50 nfast/ndtfast = 1.38889
Centers of gravity and integrals (values must be 1, 1, approx 1/2, 1, 1):
1.000000000000 1.042155637315 0.521077818658 1.000000000000 1.000000000000
Power filter parameters, Fgamma, gamma = 0.28400 0.20511
Minimum X-grid spacing, DXmin = 3.00000000E+00 km
Maximum X-grid spacing, DXmax = 3.00000000E+00 km
Minimum Y-grid spacing, DYmin = 3.01500000E-03 km
Maximum Y-grid spacing, DYmax = 3.01500000E-03 km
Minimum Z-grid spacing, DZmin = 6.29252209E-01 m
Maximum Z-grid spacing, DZmax = 2.16989607E+01 m
Minimum barotropic Courant Number = 4.71348238E+01
Maximum barotropic Courant Number = 1.05396670E+02
Maximum Coriolis Courant Number = 2.31000000E-02
Maximum grid stiffness ratios: rx0 = 8.000000E-02 (Beckmann and Haidvogel)
rx1 = 3.029772E+00 (Haney)
Initial basin volumes: TotVolume = 1.8286630435E+10 m3
MinVolume = 5.6915862271E+03 m3
MaxVolume = 1.9626709935E+05 m3
Max/Min = 3.4483725892E+01
NL ROMS/TOMS: started time-stepping: (Grid: 01 TimeSteps: 00000001 - 00086400)
STEP Day HH:MM:SS KINETIC_ENRG POTEN_ENRG TOTAL_ENRG NET_VOLUME
C => (i,j,k) Cu Cv Cw Max Speed
0 0 00:00:00 0.000000E+00 1.095119E+03 1.095119E+03 1.828663E+10
(00,000,00) 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00
DEF_HIS - creating history file: ocean_his.nc
WRT_HIS - wrote history fields (Index=1,1) into time record = 0000001
DEF_AVG - creating average file: ocean_avg.nc
DEF_DIAGS - creating diagnostics file: ocean_dia.nc
1 0 00:03:51 NaN NaN NaN NaN
(00,000,00) 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00
Blowing-up: Saving latest model state into RESTART file
WRT_RST - wrote re-start fields (Index=1,2) into time record = 0000001
Elapsed CPU time (seconds):
Thread # 0 CPU: 1.292
Total: 1.292
Nonlinear model elapsed time profile:
Allocation and array initialization .............. 0.302 (23.3794 %)
Ocean state initialization ....................... 0.033 ( 2.5643 %)
Reading of input data ............................ 0.000 ( 0.0005 %)
Processing of input data ......................... 0.000 ( 0.0122 %)
Processing of output time averaged data .......... 0.001 ( 0.0463 %)
Computation of vertical boundary conditions ...... 0.000 ( 0.0087 %)
Computation of global information integrals ...... 0.013 ( 1.0436 %)
Writing of output data ........................... 0.084 ( 6.5188 %)
Model 2D kernel .................................. 0.225 (17.3987 %)
2D/3D coupling, vertical metrics ................. 0.014 ( 1.0921 %)
Omega vertical velocity .......................... 0.007 ( 0.5286 %)
Equation of state for seawater ................... 0.012 ( 0.9484 %)
My2.5 vertical mixing parameterization ........... 0.058 ( 4.4564 %)
3D equations right-side terms .................... 0.084 ( 6.4930 %)
3D equations predictor step ...................... 0.151 (11.6745 %)
Pressure gradient ................................ 0.026 ( 1.9761 %)
Harmonic mixing of tracers, S-surfaces ........... 0.041 ( 3.1514 %)
Harmonic stress tensor, S-surfaces ............... 0.046 ( 3.5648 %)
Corrector time-step for 3D momentum .............. 0.060 ( 4.6148 %)
Corrector time-step for tracers .................. 0.084 ( 6.5136 %)
Total: 1.240 95.9859
All percentages are with respect to total time = 1.292
ROMS/TOMS - Output NetCDF summary for Grid 01:
number of time records written in HISTORY file = 00000001
number of time records written in RESTART file = 00000001
Analytical header files used:
ROMS/Functionals/ana_btflux.h
ROMS/Functionals/ana_fsobc.h
ROMS/Functionals/ana_grid.h
ROMS/Functionals/ana_initial.h
ROMS/Functionals/ana_m2obc.h
ROMS/Functionals/ana_m3obc.h
ROMS/Functionals/ana_nudgcoef.h
ROMS/Functionals/ana_smflux.h
ROMS/Functionals/ana_stflux.h
ROMS/Functionals/ana_tobc.hCode: Select all
STEP Day HH:MM:SS KINETIC_ENRG POTEN_ENRG TOTAL_ENRG NET_VOLUME
C => (i,j,k) Cu Cv Cw Max Speed
0 0 00:00:00 0.000000E+00 1.095119E+03 1.095119E+03 1.828663E+10
(00,000,00) 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00
DEF_HIS - creating history file: ocean_his.nc
WRT_HIS - wrote history fields (Index=1,1) into time record = 0000001
DEF_AVG - creating average file: ocean_avg.nc
DEF_DIAGS - creating diagnostics file: ocean_dia.nc
1 0 00:03:51 NaN NaN NaN NaN
(00,000,00) 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00- Change the timestep (make it shorter, my original values are dt=432 and ndtfast =72)
- Change my bathymetry to a flat bottom
- Change the different cpp definition
- Check the first record of the history file to check my initial conditions (bathy, temp, salt, u, v, zeta are OK). ubar and vbar are zero, which is not right because I have imposed an incoming flow (vbar) of -1.0 m/s at the upstream (north) boundary and -1.0m/s at the downstream (south) boundary in ANA_M2OBC. Could it come from the fact that I clamped a zero zeta but a non-zero 2D flow at the boundary?
Code: Select all
#elif defined COASTAL_CURRENT
IF (LBC(ieast,isUbar,ng)%acquire.and. &
& LBC(ieast,isVbar,ng)%acquire.and. &
& DOMAIN(ng)%Eastern_Edge(tile)) THEN
DO j=JstrR,JendR
BOUNDARY(ng)%ubar_east(j)=0.0_r8
END DO
DO j=Jstr,JendR
BOUNDARY(ng)%vbar_east(j)=0.0_r8
END DO
END IF
IF (LBC(iwest,isUbar,ng)%acquire.and. &
& LBC(iwest,isVbar,ng)%acquire.and. &
& DOMAIN(ng)%Western_Edge(tile)) THEN
DO j=JstrR,JendR
BOUNDARY(ng)%ubar_west(j)=0.0_r8
END DO
DO j=Jstr,JendR
BOUNDARY(ng)%vbar_west(j)=0.0_r8
END DO
END IF
IF (LBC(isouth,isUbar,ng)%acquire.and. &
& LBC(isouth,isVbar,ng)%acquire.and. &
& DOMAIN(ng)%Southern_Edge(tile)) THEN
DO i=Istr,IendR
BOUNDARY(ng)%ubar_south(i)=0.0_r8
END DO
DO i=IstrR,IendR
BOUNDARY(ng)%vbar_south(i)=-1.0_r8
END DO
END IF
IF (LBC(inorth,isUbar,ng)%acquire.and. &
& LBC(inorth,isVbar,ng)%acquire.and. &
& DOMAIN(ng)%Northern_Edge(tile)) THEN
DO i=Istr,IendR
BOUNDARY(ng)%ubar_north(i)=0.0_r8
END DO
DO i=IstrR,IendR
BOUNDARY(ng)%vbar_north(i)=-1.0_r8
END DO
END IF
The other thing I have picked up in what ROMS prints is that my minimum and maximum "barotropic Courant Number" seems to be too big, from what I could understand from other posts. But I have no idea how to correct for this, any suggestions?
Code: Select all
Minimum barotropic Courant Number = 4.71348238E+01
Maximum barotropic Courant Number = 1.05396670E+02
Maximum Coriolis Courant Number = 2.31000000E-02Mat.