Unexpected(?) blow up

[ymamoutos]

Greetings,

I am running a 14 year simulation over the Aegean sea but for some reason,
which I am struggling to understand, the model blows up after 4 years. To
be honest I am surprised because the current configuration of the model
appeared to be ok. I have run the same region for 6 years with different
atmospheric forcing (6hr NCEP) and the model run smoothly without any problem.
At the moment I am using 3hr ECMWF era-interim data. The issue that i can't
understand is that the last record of restart file has only NaNs. From
log's file last records it doesn't seem anything unusual.

Code: Select all

...

 175331  1095 19:39:00  2.216486E-03  2.547863E+03  2.547865E+03  5.064988E+13
          (133,142,30)  3.073807E-02  0.000000E+00  2.351532E+01  6.851973E-01
 175332  1095 19:48:00  2.218229E-03  2.547863E+03  2.547865E+03  5.064989E+13
          (133,142,30)  3.073789E-02  0.000000E+00  2.351516E+01  6.855414E-01
 175333  1095 19:57:00  2.219928E-03  2.547863E+03  2.547865E+03  5.064989E+13
          (133,142,30)  3.073771E-02  0.000000E+00  2.351500E+01  6.873933E-01
 175334  1095 20:06:00  2.221593E-03  2.547863E+03  2.547865E+03  5.064990E+13
          (133,142,30)  3.073753E-02  0.000000E+00  2.351484E+01  6.891778E-01
 175335  1095 20:15:00  2.223204E-03  2.547863E+03  2.547865E+03  5.064991E+13
          (133,142,30)  3.073734E-02  0.000000E+00  2.351466E+01  6.909007E-01
 175336  1095 20:24:00  2.224774E-03  2.547863E+03  2.547865E+03  5.064991E+13
          (133,142,30)  3.073714E-02  0.000000E+00  2.351449E+01  6.925654E-01
 175337  1095 20:33:00  2.226283E-03  2.547863E+03  2.547865E+03  5.064992E+13
          (133,142,30)  3.073694E-02  0.000000E+00  2.351431E+01  6.941484E-01
 175338  1095 20:42:00  2.227745E-03  2.547863E+03  2.547865E+03  5.064993E+13
          (133,142,30)  3.073672E-02  0.000000E+00  2.351412E+01  6.956356E-01
 175339  1095 20:51:00  2.229138E-03  2.547863E+03  2.547865E+03  5.064994E+13
          (133,142,30)  3.073651E-02  0.000000E+00  2.351392E+01  6.970196E-01
 175340  1095 21:00:00  2.230479E-03  2.547863E+03  2.547865E+03  5.064995E+13
          (133,142,30)  3.073630E-02  0.000000E+00  2.351371E+01  6.983050E-01
    GET_2DFLD   - surface u-wind component,                  t =  1096 00:00:00
                   (Rec=0002920, Index=1, File: aeg_wind_era_2002.nc)
                   (Tmin=        731.1250 Tmax=       1096.0000)
                   (Min = -2.31089458E+00 Max =  9.14213168E+00)
    GET_2DFLD   - surface v-wind component,                  t =  1096 00:00:00
                   (Rec=0002920, Index=1, File: aeg_wind_era_2002.nc)
                   (Tmin=        731.1250 Tmax=       1096.0000)
                   (Min = -1.24266469E+34 Max =  5.11613192E+33)
    GET_2DFLD   - surface air pressure,                      t =  1096 00:00:00
                   (Rec=0002920, Index=1, File: aeg_Pair_era_2002.nc)
                   (Tmin=        731.1250 Tmax=       1096.0000)
                   (Min =  1.00522447E+03 Max =  1.01613286E+03)
    GET_2DFLD   - cloud fraction,                            t =  1096 00:00:00
                   (Rec=0002920, Index=1, File: aeg_cloud_era_2002.nc)
                   (Tmin=        731.1250 Tmax=       1096.0000)
                   (Min =  1.54924877E-01 Max =  9.99684666E-01)
    GET_2DFLD   - surface air temperature,                   t =  1096 00:00:00
                   (Rec=0002920, Index=1, File: aeg_Tair_era_2002.nc)
                   (Tmin=        731.1250 Tmax=       1096.0000)
                   (Min =  3.44416372E+00 Max =  1.65264554E+01)
    GET_2DFLD   - surface air relative humidity,             t =  1096 00:00:00
                   (Rec=0002920, Index=1, File: aeg_Qair_era_2002.nc)
                   (Tmin=        731.1250 Tmax=       1096.0000)
                   (Min =  7.52880407E-01 Max =  9.95894587E-01)
 175341  1095 21:09:00  2.231745E-03  2.547863E+03  2.547865E+03  5.064995E+13
          (133,142,30)  3.073609E-02  0.000000E+00  2.351352E+01  6.995120E-01
 175342  1095 21:18:00           NaN           NaN           NaN           NaN
          (051,062,01)           NaN           NaN           NaN  0.000000E+00

 Blowing-up: Saving latest model state into  RESTART file

      WRT_RST   - wrote re-start fields (Index=1,1) into time record = 0000003

My DT=540, ndtfat=50, Vtransform=2, Vstretching=4, theta_s=7.0, theta_b=2.0, Tcline=200 and Hmin=2
For bathymetry smoothing i used Mathieu Dutour Sikiric LP_bathymetry package

Other values concerning Barotropic Courant number etc are

Code: Select all

 Centers of gravity and integrals (values must be 1, 1, approx 1/2, 1, 1):

    1.000000000000 1.033396903681 0.516698451840 1.000000000000 1.000000000000

 Power filter parameters, Fgamma, gamma =  0.28400   0.22720

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

 Minimum X-grid spacing, DXmin =  2.38738876E+00 km
 Maximum X-grid spacing, DXmax =  2.54891064E+00 km
 Minimum Y-grid spacing, DYmin =  2.68266913E+00 km
 Maximum Y-grid spacing, DYmax =  2.86416860E+00 km
 Minimum Z-grid spacing, DZmin =  6.60088890E-02 m
 Maximum Z-grid spacing, DZmax =  1.47386036E+02 m

 Minimum barotropic Courant Number =  5.94298262E-02
 Maximum barotropic Courant Number =  6.98174811E-01
 Maximum Coriolis   Courant Number =  5.16545468E-02

 ...

  Basin information for Grid 01:

 Maximum grid stiffness ratios:  rx0 =   6.000344E-02 (Beckmann and Haidvogel)
                                 rx1 =   3.279927E+00 (Haney)

 Initial basin volumes: TotVolume =  4.7390984047E+13 m3
                        MinVolume =  2.3233030571E+06 m3
                        MaxVolume =  8.6057460458E+08 m3
                          Max/Min =  3.7040996523E+02

My model configuration and OBC options are the above

Code: Select all

 Lateral Boundary Conditions: NLM
 ============================

 Variable               Grid  West Edge    South Edge   East Edge    North Edge
 ---------              ----  ----------   ----------   ----------   ----------

 zeta                     1   Closed       Chapman Exp  Closed       Closed

 ubar                     1   Closed       Shchepetkin  Closed       Closed

 vbar                     1   Closed       Shchepetkin  Closed       Closed

 u                        1   Closed       Clamped      Closed       Closed

 v                        1   Closed       Clamped      Closed       Closed

 temp                     1   Closed       Clamped      Closed       Closed

 salt                     1   Closed       Clamped      Closed       Closed

 tke                      1   Closed       Gradient     Closed       Closed

 Activated C-preprocessing Options:

 AEGEAN              Aegean Sea 14 years experiment
 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 nonlinear model fields.
 BULK_FLUXES         Surface bulk fluxes parameterization.
 CURVGRID            Orthogonal curvilinear grid.
 DIFF_3DCOEF         Horizontal, time-dependent 3D diffusion coefficient.
 DJ_GRADPS           Parabolic Splines density Jacobian (Shchepetkin, 2002).
 DOUBLE_PRECISION    Double precision arithmetic.
 EMINUSP             Compute Salt Flux using E-P.
 FLOATS              Simulated Lagrangian drifters.
 KANTHA_CLAYSON      Kantha and Clayson stability function formulation.
 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.
 MY25_MIXING         Mellor/Yamada Level-2.5 mixing closure.
 NONLINEAR           Nonlinear Model.
 NONLIN_EOS          Nonlinear Equation of State for seawater.
 N2S2_HORAVG         Horizontal smoothing of buoyancy and shear.
 _OPENMP             OpenMP parallel shared-memory directives.
 POWER_LAW           Power-law shape time-averaging barotropic filter.
 PROFILE             Time profiling activated .
 QCORRECTION         Surface net heat flux correction.
 K_C4ADVECTION       Fourth-order centered differences 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.
 STATIONS            Writing out station data.
 TS_U3HADVECTION     Third-order upstream horizontal advection of tracers.
 TS_C4VADVECTION     Fourth-order centered vertical advection of tracers.
 TS_DIF2             Harmonic mixing of tracers.
 TS_SMAGORINSKY      Smagorinksy-like time-dependent diffusion coefficients.
 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_QDRAG            Quadratic bottom stress.
 UV_VIS2             Harmonic mixing of momentum.
 UV_SMAGORINSKY      Smagorinksy-like time-dependent viscosity coefficients.
 VAR_RHO_2D          Variable density barotropic mode.
 VISC_3DCOEF         Horizontal, time-dependent 3D viscosity coefficient.

A question i have is based on this post viewtopic.php?f=17&t=2457&hilit=restart+file+NaN concerning Minimum Z-grid spacing and the second number
of Centers of gravity and integrals. Is it possible that these parameters create the problem. For the 6 year run these values where the same and the model run flawlessly. I am also attaching the output from the records
of restart file. Any help will be highly appreciated.

Giannis

[kate]

It's weird that it went from fine to NaN in just one step. Coincidence that it happened just after reading the forcing? I don't think so:

GET_2DFLD - surface v-wind component, t = 1096 00:00:00
(Rec=0002920, Index=1, File: aeg_wind_era_2002.nc)
(Tmin= 731.1250 Tmax= 1096.0000)
(Min = -1.24266469E+34 Max = 5.11613192E+33)

[ymamoutos]

Kate thank you so much once more.

I have check everything except from that...
The weird is that when i open the aeg_wind_era_2002.nc file
on that record (2920) doesn't have any extreme value.How
is that possible? To be honest i feel bit silly but I've never
seen anything like that. Is it a netcdf library problem? I am
using version 4.1.3

Giannis

[lanerolle]

My personal experience with NaN blow-ups is that your H+zeta at some point in your grid has become zero or < 0 for some reason or other. Do you cut-off your minimum bathymetry at some depth - say minimum bathy value is 2m? I noticed that your CPP options indicate that you do not use wetting-drying and so we can eliminate that. I would suggest that you (1) monitor/plot H+zeta and see what it does at the grid points in your domain and/or (2) play with various minimum cut-off bathy values (e.g. 2m, 3m, 5m, etc.). Minimum bathy depths close to zero should be avoided if you are not using wetting-drying.

[kate]

I doubt that it's a Netcdf library problem. You checked v-wind in your file? Can you view it with ncview? Or maybe query the min/max values for that record using a scripting language. I'd use Python and ask about record 2919 because of how Python counts.

[ymamoutos]

Greetings,

Lanerolle thanks for your advise, I'll check H+zeta and I'll make a new grid
with different hmin. Kate i don't have any problem viewing my wind file with
ncview (i'll attach a figure). The minimum value for all records is -20.12764 m/s
the maximum 18.314148 m/s and for the last record min is 0.22686242 m/s and max
11.466743 m/s. Thanks both of you again.

Giannis

[arango]

My experience with blows-up of this type after several years running is due to atmospheric forcing, specially wind stress due to strong storms. By the way, this happen to us in the Mediterranean part of our grids too. In the Aegean, there are strong winds sometimes. What we usually do is restart with a record before the model blows-up and reduce the time step and run for just a couple of days. Then, restart again at the regular time-step. Alternatively, we smooth the wind stress forcing due to strong storms that may be present during blow-up period. This implies that you are running the model at the stability limit for moderate forcing.

In your case, Kate is completely correct. The V-component of the wind is huge or all masked with a huge special value. Otherwise, ROMS input routine is reporting huge Min and Max at day 1096. By the way, that the last time available for your dataset (Tmax = 1096):

Code: Select all

    GET_2DFLD   - surface v-wind component,                  t =  1096 00:00:00
                   (Rec=0002920, Index=1, File: aeg_wind_era_2002.nc)
                   (Tmin=        731.1250 Tmax=       1096.0000)
                   (Min = -1.24266469E+34 Max =  5.11613192E+33)

You need further time records to continue running. Perhaps, the problem is in your script to generate the forcing fields NetCDF file. It wrote the U-component correctly but the V-component was set to the Fill_Value that you have.

[ymamoutos]

Hernan,

thank you very much for making clear to me all this.
The script I've used for the forcing files is d_ecmwf2roms.m.
Concerning Tmax, the total number of days for this run this
is 5114, 1096 is the last day of the aeg_wind_era_2002.nc
file.

Code: Select all

...
Input Forcing File 07:  /home/yiannis/ROMS/aegean/Input/12yrs/forcing/aeg_wind_era_2000.nc
                        /home/yiannis/ROMS/aegean/Input/12yrs/forcing/aeg_wind_era_2001.nc
                        /home/yiannis/ROMS/aegean/Input/12yrs/forcing/aeg_wind_era_2002.nc
                        /home/yiannis/ROMS/aegean/Input/12yrs/forcing/aeg_wind_era_2003.nc
                        /home/yiannis/ROMS/aegean/Input/12yrs/forcing/aeg_wind_era_2004.nc
                        /home/yiannis/ROMS/aegean/Input/12yrs/forcing/aeg_wind_era_2005.nc
                        /home/yiannis/ROMS/aegean/Input/12yrs/forcing/aeg_wind_era_2006.nc
                        /home/yiannis/ROMS/aegean/Input/12yrs/forcing/aeg_wind_era_2007.nc
                        /home/yiannis/ROMS/aegean/Input/12yrs/forcing/aeg_wind_era_2008.nc
                        /home/yiannis/ROMS/aegean/Input/12yrs/forcing/aeg_wind_era_2009.nc
                        /home/yiannis/ROMS/aegean/Input/12yrs/forcing/aeg_wind_era_2010.nc
                        /home/yiannis/ROMS/aegean/Input/12yrs/forcing/aeg_wind_era_2011.nc
                        /home/yiannis/ROMS/aegean/Input/12yrs/forcing/aeg_wind_era_2012.nc
                        /home/yiannis/ROMS/aegean/Input/12yrs/forcing/aeg_wind_era_2013.nc
...

I'll try to stop and restart the model as you proposed
and if this solution doesn't work I'll try smoothing the wind
force fields to see what is going to happen.

Thanks again and greetings from Greece.

Giannis