Dear all,

I am working on a case study of simulation the water current of estuary without wind forcing. I have created my own model grid(by using GridBuilder), boundary and initial condition(by using ROMS matlab tool). The model blew up after 1st step. The standard output showing this:

**Code:**

Basin information for Grid 01:

Maximum grid stiffness ratios: rx0 = 2.000000E-01 (Beckmann and Haidvogel)

rx1 = 5.855392E+00 (Haney)

Initial domain volumes: TotVolume = 9.9844783620E+07 m3

MinCellVol = 2.0501004235E+02 m3

MaxCellVol = 4.0155433950E+03 m3

Max/Min = 1.9587057048E+01

NL ROMS/TOMS: started time-stepping: (Grid: 01 TimeSteps: 000000000002 - 000000008641)

TIME-STEP YYYY-MM-DD hh:mm:ss.ss KINETIC_ENRG POTEN_ENRG TOTAL_ENRG NET_VOLUME

C => (i,j,k) Cu Cv Cw Max Speed

1 0001-01-01 00:00:01.00 5.573662E-02 1.564449E+02 1.565007E+02 4.488730E+07

(013,017,01) 1.373405E-02 9.012729E-03 0.000000E+00 7.432341E-01

DEF_HIS - creating history file, Grid 01: ocean_his.nc

2 0001-01-01 00:00:02.00 NaN NaN NaN NaN

(049,043,01) 6.110394E-03 7.331021E-03 3.991203E-02 7.416617E-01

Found Error: 01 Line: 298 Source: ROMS/Nonlinear/main3d.F

Found Error: 01 Line: 298 Source: ROMS/Drivers/nl_ocean.h

Below is my CPP options:

**Code:**

ESTUARYNW Estuary current simulation without wind force

ANA_BSFLUX Analytical kinematic bottom salinity flux

ANA_BTFLUX Analytical kinematic bottom temperature flux

ANA_SMFLUX Analytical kinematic surface momentum flux

ANA_SSFLUX Analytical kinematic surface salinity flux

ANA_STFLUX Analytical kinematic surface temperature flux

ASSUMED_SHAPE Using assumed-shape arrays

AVERAGES Writing out time-averaged nonlinear model fields

!BOUNDARY_ALLGATHER Using mpi_allreduce in mp_boundary routine

CANUTO_A Canuto A-stability function formulation

!COLLECT_ALL... Using mpi_isend/mpi_recv in mp_collect routine

CURVGRID Orthogonal curvilinear grid

DEBUGGING Internal debugging switch activated

DIAGNOSTICS_UV Computing and writing momentum diagnostic terms

DOUBLE_PRECISION Double precision arithmetic numerical kernel.

GLS_MIXING Generic Length-Scale turbulence closure

HDF5 Creating NetCDF-4/HDF5 format files

INLINE_2DIO Processing 3D IO level by level to reduce memory needs

MASKING Land/Sea masking

MPI MPI distributed-memory configuration

NONLINEAR Nonlinear Model

!NONLIN_EOS Linear Equation of State for seawater

N2S2_HORAVG Horizontal smoothing of buoyancy and shear

OUT_DOUBLE Double precision output fields in NetCDF files

PARALLEL_IO Parallel I/O processing

PERFECT_RESTART Processing perfect restart variables

POWER_LAW Power-law shape time-averaging barotropic filter

PRSGRD31 Standard density Jacobian formulation (Song, 1998)

PROFILE Time profiling activated

K_C2ADVECTION Second-order centered differences advection of TKE fields

RADIATION_2D Use tangential phase speed in radiation conditions

READ_WATER Reading data at water points only

REDUCE_ALLGATHER Using mpi_allgather in mp_reduce routine

RI_SPLINES Parabolic Spline Reconstruction for Richardson Number

RHO_SURF Include difference between rho0 and surface density

!RST_SINGLE Double precision fields in restart NetCDF file

SOLVE3D Solving 3D Primitive Equations

TS_C4HADVECTION Fourth-order centered horizontal advection of tracers

TS_C4VADVECTION Fourth-order centered vertical advection of tracers

UV_ADV Advection of momentum

UV_COR Coriolis term

UV_C2ADVECTION Second-order centered differences advection of momentum

UV_QDRAG Quadratic bottom stress

VAR_RHO_2D Variable density barotropic mode

VISC_3DCOEF Horizontal, time-dependent 3D viscosity coefficient

After reading some posts that having the same problem, I tried:

1)shorten my DT from 10 seconds to 1 second. The standard output is shown below. The maximum barotropic Courant Number is smaller than 1.

**Code:**

Minimum barotropic Courant Number = 1.27269811E-02

Maximum barotropic Courant Number = 4.75991568E-02

Maximum Coriolis Courant Number = 9.50472208E-05

2)check if bathymetry requiring smoothing:

The standard output is shown below. It seems the rx0 and rx1 satisfying 0 < rx0 < 0.4 and 3 < rx1 < 7

**Code:**

Minimum X-grid spacing, DXmin = 4.01181852E-02 km

Maximum X-grid spacing, DXmax = 4.01577087E-02 km

Minimum Y-grid spacing, DYmin = 4.02366393E-02 km

Maximum Y-grid spacing, DYmax = 4.02661984E-02 km

Minimum Z-grid spacing, DZmin = -2.17046894E-01 m

Maximum Z-grid spacing, DZmax = 2.48532885E+00 m

...

...

Horizontal mixing scaled by grid size, GRDMAX = 4.02119170E-02 km

Minimum horizontal viscosity coefficient = 0.00000000E+00 ▒▒▒

Maximum horizontal viscosity coefficient = 0.00000000E+00 ▒▒▒

...

...

...

Basin information for Grid 01:

[b] Maximum grid stiffness ratios: rx0 = 2.000000E-01 (Beckmann and Haidvogel)

rx1 = 5.855392E+00 (Haney)[/b]

3)check if boundary/initial/grid nc file has NAN values:

I didn't find nan value.

I tried to plot the restart file and come up with some questions..

1) The figures for temperature and salinity looks ok, same as my initial conditions. However, figures for u and v looks weird, it seems the mask doesn't work. I attached the screenshot of it.

The temperature and salinity in restart files is NaN in land area, while u and v is 0 on land area. Since I actived "READ_WATER" option, which said 'use if only reading water points data', all the value on land should be NaN in restart file (am I right?). I am a little confused about it.

2) I mentioned the dimension of some variables (temp, salt, u and v) in ocean_rst.nc file, is 'Dimensions: xi_v,eta_v,s_rho,two,ocean_time'. What does the 'two' mean? Is this because that I set NAT=2?

3) I saw in one post, people suggested to check the forcing file if ROMS blow up at 1st step. I don't have forcing file right now, only set five analytical options in header file. They are:

**Code:**

#define ANA_SMFLUX

#define ANA_STFLUX

#define ANA_SSFLUX

#define ANA_BTFLUX

#define ANA_BSFLUX

Could this be the reason of blowing up? Do I have to create a ocean_frc.nc file for ROMS, and set the wind field to be zero?

Thanks for your attention and help!

Best regards

Yifan