potential mass conservation problem on dry cells

[schen]

Hi all,

There may be a potential problem in the ROMS fill value (spval = 1.0E+37) for dry cells when wetting/drying is activated.

In Nonlinear/wetdry.F, when the local water depth (zeta + h) < Dcrit, the cell is considered dry, and rmask_wet and rmask_full of that cell is set to zero. Then, say for writing out zeta (in Utility/wrt_his.F), rmask_full is passed to nf_fwrite2d.F so that zeta at the dry cell is set equal to the 1.0E+37.

But, there is still small water mass at the dry cell, and the water mass is not included when computing global mass conservation. This can potentially introduce large error when the modeled domain has a large intertidal area. In my simple beach test case, I did see that total mass decreases a little bit during ebbs. But the error is relatively small due to the small intertidal area.

Is my statement above correct? If so, then we may want to write out the computed zeta at the dry cells instead of filling in 1.0E+37. Suggestion is appreciated!!

---- Shih-Nan

[arango]

Yes, good catch. We need to avoid the _FillValue attribute and masking when WET_DRY is activated. This is avoided easily by adding the optional SetFillVal = .FALSE. argument to some of the I/O processing routines.

Please update. I hope that this solves your mass conservation problems.

[schen]

Thank John and Hernan for implementing the corrections. I have tested it out, and zeta + h at dry cells is < Dcrit now. Plus, total mass conserves very well.

[dhson]

I've tried to use WET_DRY for a small island that its top is at about 2 m above MSL. To reserve its original topography, I've got to use hmin = -2 m (in SEAGRID). Start computing with zeta =0 and I 've expected that most of cells will be "WET" in flood tide. But I couldn't compute with any value of hc.

Probably, the main problem is how to set correct values of hmin, Hc, and Tcline.

Would any body get me out? Thanks

[jcwarner]

so what is the problem? Does it blow up on first time step? does it initialize correctly? There should not be a limit on himn.Try tcline = hc =0.

[dhson]

Sorry, with hc=.1m, hmin=-5m, I've got (SVN Revision : 342M):

NETCDF_CHECK_VAR - inconsistent value of variable: hc 1.000000E-01 -5.000000E+00
in file: ROMS_FILES/roms_ini.nc

with hc=0, I've got the same results.

If I use hmin=.1m in SEAGRID, I can run ROMS but that's not what I want since SEAGRID will flatten all grids above hmin to hmin.

[jcwarner]

viewtopic.php?f=13&t=1233

[dhson]

Thanks for the information, but not much help.
Using
Vtransform == 1 ! transformation equation
Vstretching == 1 ! stretching function
with the combination hmin=-5m (as explained above),
and (hc=-5m, Tcline=0m) or (hc=-5m, Tcline=-5m) I can run but:

STEP Day HH:MM:SS KINETIC_ENRG POTEN_ENRG TOTAL_ENRG NET_VOLUME

0 0 00:00:00 0.000000E+00 1.261459E+02 1.261459E+02 8.199222E+10
DEF_HIS - creating history file: slr_his_0001.nc
WRT_HIS - wrote history fields (Index=1,1) into time record = 0000001
DEF_STATION - creating stations file: slr_sta.nc
1 0 00:01:00 NaN NaN NaN NaN

Blowing-up: Saving latest model state into RESTART file

Thanks in advance,

[jcwarner]

well, it does seem to be help since you are now past that problem and onto the next problem. So i ask that you look ito the history file and see if all looks ok.

"Blowing-up: Saving latest model state into RESTART file"
So now look into that RESTART file and see where did it blow up? At the boundary? At some interior point? You need to use your deductive skills here, be a detective.

[kate]

There are some brand new check in ROMS on reading in Tcline and friends. It may be that they are overly restrictive in your case.

[dhson]

I think probably compare Hmin to Tcline and hc is not a good idea since Hmin is elevation and the later are thickness.
I've commented out both Tcline and hc in mod_netcdf.F:

! exit_flag=5
! EXIT

but no improvement in ROMS.

[jcwarner]

"Blowing-up: Saving latest model state into RESTART file"
So now look into that RESTART file and see where did it blow up? At the boundary? At some interior point? You need to use your deductive skills here, be a detective.

[dhson]

Thank you for the response,

RESTART file did not follow any rule with differences of dt and all blow-up in the first step:
- with dt of several second: patchy pattern of NaN and real value (*.E-6) without specific relation with topography;
- with smaller dt (0.001 seconds): no NaN observed but one very big value (*.E240) somewhere in the middle of domain also without specific relation with topography.

ROMS was compiled with MPICH2, gfortran 4.3, netcdf4 on x86_64 GNU/Linux 2.6.28-11.

I've used this version of ROMS for other configuration (positive Hmin, no WET_DRY), It worked well.

[jcwarner]

"with smaller dt (0.001 seconds): no NaN observed but one very big value (*.E240) somewhere in the middle of domain also without specific relation with topography."

Lets keep pushing on this. Where? What does bathy look like there ? Did it create a history file and write the initial conditions for time step 0 into the his file before it blew up? So what did that area look like in the his file? Can you make some plots?
Was zeta always > h ?
keep digging.

[dhson]

Hi! It's me again,
Since my domain was a bit big (200x600) so I was afraid of some uncontrolled grid cells. It took me sometime to prepared more ideal bathymetry as attached with the same size (so in thi run only netcdf grid file was different).
ROMS died at the 2nd timestep (a step better than the original)
zeta, ubar, vbar in RESTART and my impression also included there.
My header file as:
#undef WRF_COUPLING



#undef SWAN_COUPLING

#ifdef SWAN_COUPLING

# define MCT_LIB

#endif

#define UV_ADV

#define UV_COR

#define UV_VIS2

#define MIX_S_UV



#define TS_U3HADVECTION

#define TS_C4VADVECTION

#define TS_DIF2

#define MIX_S_TS

#undef TS_MPDATA

#define DJ_GRADPS



#undef UV_PSOURCE

#undef TS_PSOURCE



#define SOLVE3D

#define SALINITY

#define NONLIN_EOS

#define CURVGRID

#define SPLINES

#define MASKING

#define STATIONS

#define WET_DRY



#define WESTERN_WALL



#undef SPONGE

#define RADIATION_2D



#define RAMP_TIDES

#define SSH_TIDES

#ifdef SSH_TIDES

# define ADD_FSOBC

# define EAST_FSCHAPMAN

# define NORTH_FSCHAPMAN

# define SOUTH_FSCHAPMAN

#else

# define EAST_FSGRADIENT

# define NORTH_FSGRADIENT

# define SOUTH_FSGRADIENT
#endif

#define UV_TIDES

#ifdef UV_TIDES

# define ADD_M2OBC

# define EAST_M2FLATHER

# define NORTH_M2FLATHER

# define SOUTH_M2FLATHER
#else

# define EAST_M2RADIATION

# define NORTH_M2RADIATION

# define SOUTH_M2RADIATION
#endif

#if defined SSH_TIDES || defined UV_TIDES

# undef EAST_VOLCONS

# undef NORTH_VOLCONS

# undef SOUTH_VOLCONS

#else

# define EAST_VOLCONS

# define NORTH_VOLCONS

# define SOUTH_VOLCONS

#endif



#define EAST_M3RADIATION

#define EAST_TRADIATION

#define EAST_TNUDGING



#define NORTH_M3RADIATION

#define NORTH_TRADIATION

#define NORTH_TNUDGING



#define SOUTH_M3RADIATION

#define SOUTH_TRADIATION

#define SOUTH_TNUDGING


/* define only one of the following 5 */

#undef UV_LOGDRAG

#define UV_QDRAG

#undef MB_BBL

#undef SG_BBL



#ifdef SOLVE3D

# undef LMD_MIXING

# ifdef LMD_MIXING

# define LMD_RIMIX

# define LMD_CONVEC

# define LMD_SKPP

# define LMD_NONLOCAL

# endif

# undef MY25_MIXING

# ifdef MY25_MIXING

# define N2S2_HORAVG

# define KANTHA_CLAYSON

# undef K_C2ADVECTION

# undef K_C4ADVECTION

# endif

# define GLS_MIXING

# ifdef GLS_MIXING

# define KANTHA_CLAYSON

# define N2S2_HORAVG

# endif

# undef SEDIMENT

# ifdef SEDIMENT

# define SUSPLOAD

# undef BEDLOAD_SOULSBY

# undef BEDLOAD_MPM

# define SED_MORPH

# endif

# if defined SEDIMENT || defined SG_BBL || defined MB_BBL || defined SSW_BBL

# define ANA_SEDIMENT

# endif


# undef BULK_FLUXES
# ifdef BULK_FLUXES
# define LONGWAVE
# define ANA_CLOUD
# define ANA_RAIN
# endif


# define ANA_BSFLUX

# define ANA_SSFLUX

# define ANA_BTFLUX

#endif



#define QCORRECTION

#define SOLAR_SOURCE



#undef BIO_FASHAM

#ifdef BIO_FASHAM

# define CARBON

# define DENITRIFICATION

# define BIO_SEDIMENT

# undef DIAGNOSTICS_BIO

#endif

#undef ECOSIM

#if defined BIO_FASHAM || defined ECOSIM

# define ANA_BIOLOGY

# define ANA_SPFLUX

# define ANA_BPFLUX

#endif

[jcwarner]

OK. The way wet_dry works is that the free surface is draped over the topography. So you had a comment "strange zeta: the left haft of the domain shouldn't be covered with water"
Well, what is happening is that the water level is being set to be greater than the depth, by an amount of Dcrit. For further description please seee here-->
https://www.myroms.org/wiki/index.php/WET_DRY
(our wet_dry paper is almost complete. it should be out soon).
It looks to me that the problem is along the northern (left ) boundary. There could be issues with certain boundary conditions. With all the possible combinations of types of BCs, I can not check them all. So maybe you have a configuration that was not tested. I recently (yesterday) posted a fix for set_data - what was happening is that the user was imposing a BC of zeta that was less than the dpeth. So i had to code in a catch to prevent that. It sounds like a similar issue to what you have. This wet_dry option is relatively new. So maybe you are finding something that others did not cathch yet. Sorry. But can you at least update the ROMS/Nonlinear/set_data.F routine and see if that works. IF not, I may need to grab your input files and take a quick look.
-j

[dhson]

wow! It works!!!
But I believed that there're still some bugs related to BC (when the water level is physically below bottom of the boundary cells).
The attached are output of HISTORY after 1 hour (in fact ROMS now can lives much more than that).

Thanks

[jcwarner]

You are making progress - good. But we are not done yet.
That northern boundary is still not looking correct. In the figure, vbar along the northern edge, left side, is negative (in a triangle region). So is the water level along that side much greater than the land? If so, then water will come in. Like a tidal wave. I can not prevent water from coming in. You must do that. Is it realistic for the water level to raise above the land there on the left side?
If it is not realistic, then you must not let it happen. This can be dealt with in several ways. 1) use land mask = 0 so that those cells never get wet. 2) limit the boundary zeta_north to not go above the land (just on the left side). The latest 'catch' i coded in should then set the water level to be Dcrit higher than the land, and all should be ok.

[dhson]

Well, I didn't make any specific condition along all boundaries. They all start from initial condition (zeta=0) then driven by tidal constants.
If BC at wet_dry cells was coded correctly then, Probably, I've got to check tidal constants that were interpolated form TPXO7.1. I'm not sure how the inland cells (dry) were dealt with.

So thank you and if you have any suggestion pls let me know.

[dhson]

Hi!
I didn't make ROMS run longer 1day and 15 hours.
Our tidal amplitude is <1m. That means (if you said correctly "the water level to raise above the land" ) the water can come only on boundary cells with h>-1m. But my figure indicated that the water come from all boundary cells.
Even though, as you suggested, I've set all north and south boundary cells with h<0m to land to prevent water come from.
ROMS still was blow up from not very special cell (maybe the shallowest among h-positive cells but not at the boudaries).
I've also tried to change dt, viscosity but unsuccessful.
Any comment or suggestion?

[jcwarner]

i found an 'issue' with v2dbc_im.F. There is a RETURN call that prevents execution of the wet_dry boundary calls. I posted a fix for this. Should be out soon, and then please give it a try - let me know if it works.