ROMS is blowing up after just one time step

[mdessert]

Hi modellers,
I'm trying to model the ocean dynamics in a “small” domain along the coast of West South Africa in the Benguela Current. I chose it a bit randomly as a first test domain as I want to study other domains afterwards. At the fist time step (t=1), ROMS blows up while indicating that :

Code: Select all

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

 Blowing-up: Saving latest model state into  RESTART file

Just above, it indicates that :

Code: Select all

 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

         0 2018-01-01 12:00:00.00  3.072290E-03  2.027491E+04  2.027491E+04  2.728826E+15
                       (36,07,70)  2.846040E-02  1.058142E-02  0.000000E+00  5.856957E-01
      DEF_HIS     - creating  history      file, Grid 01: ocean_his.nc
      WRT_HIS     - wrote history     fields (Index=1,1) in record = 0000001
         1 2018-01-01 12:15:00.00           NaN           NaN           NaN           NaN
                       (00,00,00)  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00

I first tried to diminish the time step from 5 minutes to 60 secondes. But that was not efficient.
Then, I thought about the Courant Number.

Code: Select all

 Metrics information for Grid 01:
 ===============================

 Minimum X-grid spacing, DXmin =  1.56642103E+01 km
 Maximum X-grid spacing, DXmax =  1.68372871E+01 km
 Minimum Y-grid spacing, DYmin =  2.42245568E+01 km
 Maximum Y-grid spacing, DYmax =  2.60174165E+01 km
 Minimum Z-grid spacing, DZmin =  1.36165981E-01 m
 Maximum Z-grid spacing, DZmax =  1.14109522E+02 m

 Minimum barotropic Courant Number =  2.10206234E-02
 Maximum barotropic Courant Number =  5.07703716E-01
 Maximum Coriolis   Courant Number =  6.13830395E-02

But, I've read somewhere on the forum that having a Courant number <1 was enough.

As my maximum Z-grid spacing was around several hundreds of meters and that I've only seen examples where this value was lower, I've tested to run ROMS with a more bottomward accurated grid (modifying thetab). But that did not solve my problem neither.

Finally, ROMS wrote that the roughness of my grid was higher than those I've seen in examples.

Code: Select all

Maximum grid stiffness ratios:  rx0 =   9.925375E-01 (Beckmann and Haidvogel)
                                 rx1 =   1.301136E+02 (Haney)

So, I examined my grid more precisely. I first note that the depth increases very rapidly near the coast. (I was wondering if that could be caused by my Hmin which is very low -10meters-).



Secondly, while the dmde -not showed- (which stands for "XI derivative of inverse metric factor pn") changes smoothly, the dndx (which stands for 'ETA derivative of inverse metric factor pm') shows some abrupt modifications.



I am now wondering if my problem could come from my bathymetry and eventually if you add some clues for me to seek better (or eventually, I could choose a more simple domain -a little southward- for a first test) or if I am looking hard in a totally bad direction…

I've read a lot of posts about similar problems, but that did not solved mine. I hope you could give me some clues...
Thanks a lot!

[kate]

Well, it could still be many things. You really should smooth your bathymetry, but that might not be the only problem.

[mdessert]

Thank you Kate for your answer.
I tried to smooth (even exageratedly) but, as you expected, that does not solve anything.
I tried to compute an other domain with a smoother bathymetry as joining the coast but that does not change anything neither.
With a colleague, we have noticed that my ocean_his.nc (that ROMS just computes before it blows up) showed Akt and Akv =0 at the surface but Nan in the lower levels (easily seen with a ncdump ocean_his.nc | less or with ferret).

I'm looking for something as scrutizing the code. And I have seen that there are two different variables named Akt and AKt. Especially, in lmd_skpp.F at line 897, I've found :

Code: Select all

              Akv(i,j,k)=depth*wm(i,j)*(1.0_r8+sigma*Gm)
              Akt(i,j,k,itemp)=depth*ws(i,j)*(1.0_r8+sigma*Gt)
#  ifdef SALINITY
              AKt(i,j,k,isalt)=depth*ws(i,j)*(1.0_r8+sigma*Gs)
#  endif

Is there a problem in my code? Or does the both names really stand for different variables in this file?

[kate]

Your code is fine - Akt is a four-dimensional array and you are setting different parts of it. Having zero at the surface is to be expected, NaN values below that is not OK.

[mdessert]

That does not matter the fact that the K from (Akt) is in capital letter in one line and not in the other?
(Anyway, I tested with Akt at the place of AKt, and that does not solve my problem... )
I keep on searching...

[jivica]

Have you checked your initial condition's file values?
This line looks really suspicious:
1 2018-01-01 12:15:00.00 NaN NaN NaN NaN

Do you have missing values in the ini.nc? What about bry.nc? atmo.nc? you have to check all those files to be sure you are not giving NaN to ROMS model...

Cheers
I.

[mdessert]

Thank you for your advice!
None of these files seems to be weird (No Nan, and values are coherent).
But I noticed something concerning my wind input. On the netCDF in the initial file, all values seem OK. But as I make my ROMS model write the values of wind it reads in netCDF (adding WRITE(*,*) in the bulk_flux.F), it shows me that the wind is zero everywhere...
I suspect a problem as reading the wind forcing...

Actually, I thought this line suspicious too. So, I added WRITE(*,*) in several lines (one after the other) in my routines to understand from where these values could become Nan. I went back up from the kinetic energy (the variable which is showed -if i understood well- in my error message from ROMS) to the reading of the null wind... And There I am now!

I tried to add these lines in my bulk_flux.F routine :

Code: Select all

      WRITE(*,*) ' Test on wind ',FORCES(ng) % Uwind
      WRITE(*,*) ' Test on temperature ',FORCES(ng) % Tair

And Roms returns me null fields (for both temperature and Uwind)...
I realized too that ROMS tried to read inputs for a weird date : 2434 ?!?

GET_2DFLD - surface u-wind component, 2434-01-23 00:00:00.00
[...]
GET_2DFLD - surface v-wind component, 2434-01-23 00:00:00.00
[...]
GET_2DFLD - surface air pressure, 2434-01-23 00:00:00.00
[...]
GET_2DFLD - solar shortwave radiation flux, 2434-01-23 00:00:00.00
[...]
GET_2DFLD - downwelling longwave radiation flux, 2434-01-23 00:00:00.00
[...]
GET_2DFLD - surface air temperature, 2434-01-23 00:00:00.00
[...]
GET_2DFLD - surface air relative humidity, 2434-01-23 00:00:00.00
[...]
GET_2DFLD - rain fall rate, 2434-01-23 00:00:00.00

.. whereas my netcdf file shows time units that seems correct...
Can it be caused by the units I used? "hours since 01-01-2000" ?

[kate]

ROMS cannot parse your time units string. It looks for "seconds" or "days", not hours. One could change that, of course.

[mdessert]

Hello modellers,
I've made some modifications to ROMS take into account the "hour" units as time axis. Now, here is what I've added.

In mod_scalars.F :
I added line 705

real(r8), parameter :: day2sec = 86400.0_r8
real(r8), parameter :: hour2day = 1.0 / 24.0_r8 ! LINE 705
real(r8), parameter :: sec2day = 1.0_r8 / 86400.0_r8

In check_multifile.F :
I added line 655

Tunits=TRIM(var_Achar(i))
IF (INDEX(TRIM(var_Achar(i)),'day').ne.0) THEN
Tscale=day2sec
ELSE IF (INDEX(TRIM(var_Achar(i)),'hour').ne.0) THEN ! LINE 655
Tscale=hour2sec
ELSE IF (INDEX(TRIM(var_Achar(i)),'second').ne.0) THEN
Tscale=1.0_r8
END IF

In get_state.F :
I added line 311

IF (INDEX(TRIM(var_Achar(i)),'day').ne.0) THEN
time_scale=day2sec
ELSE IF (INDEX(TRIM(var_Achar(i)),'hour').ne.0) THEN ! LINE 311
time_scale=hour2sec
ELSE IF (INDEX(TRIM(var_Achar(i)),'second').ne.0) THEN
time_scale=1.0_r8
ELSE
IF (Master) WRITE(stdout,11) TRIM(var_Achar(i)),TRIM(tvarnam)
exit_flag=4 ! Here I was not sure about the value of exit_flag to use
END IF

[…]
and at the end of the file

11 FORMAT (/,' GET_STATE - unable to get units for attribute: ',a, &
& /,13x,'in variable: ',a)

Now, I could not have yet to test totally these added code because my simulation does not yet work. I still have a problem as reading my atmospheric inputs.

Actually, I've isolated my problem in the reading of longitude_rho and latitude_rho by the nf_nfread2d.F in regrid.F (I use ONLINE option to interpolate my atmosphere input by ROMS. In my atmosphere input files, my latitudes were written decreasing (from -10 to -44)...

[zhaoqiang]

When I face such problem of 1st-step-blow-up, the first thing I do is to check all the forcing files to see if there is any NaN in them.

[mdessert]

Thank you for your advices.

Actually, I had some problems with my atmospheric input files because of interpolation online. In fact, my input files were written with the latitude decreasing instead of increasing. So when ROMS was scanning (in subroutine hindices in interpolate.F) my latitude values to find relations between atmospheric grid and destination grid, it does not find and the interpolation computed atmospheric fields to 0.
I modified my initial netcdf file (as atmospheric files) by reversing the latitudes axis and its does not blow up at all now!

[mdessert]

Hi everyone,
I keep working on my configuration and I'm facing for several days a new blowing-up event after four time step (one time step is 180 s).
I look deeper into my output files (netcdf and log error) and I had some clues. But I'm afraid that I have not the keys to understand well.

Actually, I wondered whether that could be caused by my grid (especially my vertical definition).

My grid was built with Vtransform=2 (and Vstretching=4) and I choose ThetaS=3 and Theta_B = 2 to increase my vertical resolution on surface. I have 70 vertical levels (but I tried also with 50 or 60). My vertical decomposition is something like that :

(bottom is on the left and surface on the right)

ROMS model wrote me :

Metrics information for Grid 01:
===============================

Minimum X-grid spacing, DXmin = 2.39340275E+00 km
Maximum X-grid spacing, DXmax = 2.78118434E+00 km
Minimum Y-grid spacing, DYmin = 2.28795991E+00 km
Maximum Y-grid spacing, DYmax = 2.65497320E+00 km
Minimum Z-grid spacing, DZmin = 6.97149813E-02 m
Maximum Z-grid spacing, DZmax = 1.20089084E+02 m

Minimum barotropic Courant Number = 1.17094996E-02
Maximum barotropic Courant Number = 4.07452685E-01
Maximum Coriolis Courant Number = 1.68657029E-02


Minimum horizontal diffusion coefficient = 1.06786179E+08 m4/s
Maximum horizontal diffusion coefficient = 1.67206682E+08 m4/s

Minimum horizontal viscosity coefficient = 1.06786179E+08 m4/s
Maximum horizontal viscosity coefficient = 1.00000000E+20 m4/s

As studying my output netcdf file, I realized that my u-velocity increase anormally (more than 40m/S) at the ground of my ocean (but the v-velocity stay normal). The profile at some points I identify show some instability from the bottom which provoke "oscillation" (vertical level to vertical level, u-velocity changes from -40m/S to 40 m/S and then -38m/s)... These points take place near modification of bathymetry (seamounts or near plateau breaking).

(bottom is on the left and surface is on the right)

I have activated the LMD mixing scheme with the definition of surface and bottom bondary layers :

LMD_BKPP KPP bottom boundary layer mixing.
LMD_CONVEC LMD convective mixing due to shear instability.
LMD_MIXING Large/McWilliams/Doney interior mixing.
LMD_RIMIX LMD diffusivity due to shear instability.
LMD_SKPP KPP surface boundary layer mixing.

(I attached the whole ocean.h)

I suspect many things but I think I need some helps to point the more adapted solution...
Is my Courant number too high? I read on the forum that it must be lower than 0.1...
Are my vertical layer too narrowed in surface and not enough in depths?
Are my coefficient of diffusivity (or viscosity) something like AKt (or AKv)? In this case they look too high, don't they? (In my netcdf output file, I saw AKV values which reach 0.175)

[kate]

Using quadratic bottom drag without a drag limiter can lead to such troubles. Try #define LIMIT_BSTRESS.

[mdessert]

Hi Kate, thank you for your answer.
I've tried to activate LIMIT_BSTRESS as you suggested. Unfortunately, that does not work either. (It blows up with the same abnormal velocity values).
I then tried # define LIMIT_VVISC
#define RI_SPLINES
I tried to change my drag scheme (with #define UV_LOGDRAG and then UV_LDRAG).

But test after test, that fails...

Moreover, I add some advice/warning that ROMS write back to me :

CHECKADJ - use caution when activating: TS_U3ADV_SPLIT
REASON: stability problems, WARNING.

CHECKADJ - use caution when activating: UV_U3ADV_SPLIT
REASON: stability problems, WARNING.

As I try to simulate coastal regions, I choose this advection scheme because of a steep bathymetry. I understood it was the best solution. But, i now have a doubt...

Here are my cpp keys : (extracted from my log file)

Activated C-preprocessing Options:

CARETTA Caretta
ANA_BSFLUX Analytical kinematic bottom salinity flux.
ANA_BTFLUX Analytical kinematic bottom temperature flux.
ANA_SPONGE Analytical enhanced viscosity/diffusion sponge.
ASSUMED_SHAPE Using assumed-shape arrays.
AVERAGES Writing out time-averaged nonlinear model fields.
!BOUNDARY_ALLGATHER Using mpi_allreduce in mp_boundary routine.
BULK_FLUXES Surface bulk fluxes parameterization.
!COLLECT_ALL... Using mpi_isend/mpi_recv in mp_collect routine.
CURVGRID Orthogonal curvilinear grid.
DIAGNOSTICS_TS Computing and writing tracer diagnostic terms.
DIAGNOSTICS_UV Computing and writing momentum diagnostic terms.
DIFF_3DCOEF Horizontal, time-dependent 3D diffusion coefficient.
DOUBLE_PRECISION Double precision arithmetic.
EMINUSP Compute Salt Flux using E-P.
LIMIT_BSTRESS Limit bottom stress to maintain bottom velocity direction.
LIMIT_VVISC Impose an upper limit on vertical viscosity coefficient.
LMD_BKPP KPP bottom boundary layer mixing.
LMD_CONVEC LMD convective mixing due to shear instability.
LMD_MIXING Large/McWilliams/Doney interior mixing.
LMD_RIMIX LMD diffusivity due to shear instability.
LMD_SKPP KPP surface boundary layer mixing.
LONGWAVE_OUT Compute outgoing longwave radiation internally.
MASKING Land/Sea masking.
MIX_ISO_TS Mixing of tracers along isopycnal surfaces.
MIX_GEO_UV Mixing of momentum along geopotential surfaces.
MPI MPI distributed-memory configuration.
NONLINEAR Nonlinear Model.
NONLIN_EOS Nonlinear Equation of State for seawater.
POWER_LAW Power-law shape time-averaging barotropic filter.
PRSGRD31 Standard density Jacobian formulation (Song, 1998).
PROFILE Time profiling activated .
RADIATION_2D Use tangential phase speed in radiation conditions.
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.
SALINITY Using salinity.
SOLVE3D Solving 3D Primitive Equations.
SPLINES_VDIFF Parabolic Spline Reconstruction for Vertical Diffusion.
SPHERICAL Spherical grid configuration.
TS_C4HADVECTION Fourth-order centered horizontal advection of tracers.
TS_C4VADVECTION Fourth-order centered vertical advection of tracers.
TS_U3ADV_SPLIT Split third-order upstream advection of tracers.
TS_DIF4 Biharmonic mixing of tracers.
UV_ADV Advection of momentum.
UV_COR Coriolis term.
UV_C4ADVECTION Fourth-order centered differences advection of momentum.
UV_U3ADV_SPLIT Split third-order upstream advection of momentum.
UV_LOGDRAG Logarithmic bottom stress.
UV_VIS4 Biharmonic mixing of momentum.
VAR_RHO_2D Variable density barotropic mode.
VISC_3DCOEF Horizontal, time-dependent 3D viscosity coefficient.

[ymamoutos]

Greetings

Do you have any specific reason to use
the 3rd-order upstream split advection
scheme for momentum and tracers? I am
asking because I recall that whenever
it tried to use these options in horizontal
resolutions similar or higher than yours
- 2km to 1km - I had stability issues.

Undefine them and use the default option
(or AKIMA for tracres?) for the advection
plus the harmonic horizontal mixing option
(TS_DIF2 and UV_VIS2). Recompile and try
to run the model.

If you must use the 3rd-order upstream
split advection scheme change MIX_ISO_TS
to MIX_GEO_TS as is recommended to ccpdefs.h
file.

Giannis

[mdessert]

Thank you for your advice.

Actually, I understood I had to use the advection scheme if I wanted the more efficient solution for configuration with steep bathymetry, following some advices from K. Hedstrom. 2016. Technical Manual for a Coupled Sea-Ice/Ocean Circulation Model (Version
4). U.S. Dept. of the Interior, Bureau of Ocean Energy Management, Alaska OCS Region. OCS
Study BOEM 2016-037. 176 pp.
and in Wikiroms from roms.org
(But I may have misinterpreted...)

The 3rd-order upstream split advection (TS_U3ADV_SPLIT) can be used to cor-
rect for the spurious diapycnal diffusion of the advection operator in terrain-following
coordinates. If this is chosen, the advection operator is split in advective and diffusive
components and several internal flags are activated in globaldefs.h. Notice that hori-
zontal and vertical advection of tracer is 4th-order centered plus biharmonic diffusion to
correct for spurious diapycnal mixing. The total time-dependent horizontal mixing coef-
ficient are computed in hmixing.F. It is also recommended to use the rotated mixing
tensor along geopotentials (MIX_GEO_TS) for the biharmonic operator.

and

The 3rd-order upstream split advection (UV_U3ADV_SPLIT) can be used to correct
for the spurious mixing of the advection operator in terrain-following coordinates. If this
is chosen, the advection operator is split into advective and viscosity components and
several internal flags are activated in globaldefs.h. Notice that horizontal and vertical
advection of momentum is 4th-order centered plus biharmonic viscosity to correct for
spurious mixing.

I realized that I was not enough careful as I miss the #define MIX_GEO_TS

I add it but that blowed up the same way...

As you explained, I changed my advection scheme tot the default one and that's running (my ocean is computed for 14 time steps! ! ! )

In this case, I have a question : what are actually the TS_U3ADV_SPLIT and UV_U3ADV_SPLIT options made for?

Morgane