WET_DRY and Gridbuilder Problem

[zhangtianhao]

Hi there,

I met a problem with WET_DRY option.

When I use WET_DRY to run a simple coastal ocean case, ROMS blows up quickly (19 steps, 1.9s). I set DCRIT=0.1m.

Notice: I use Gridbuilder and set h of land area as negative value(though people usually set all h>0), for example,-15m. And rx0 and rx1 is very fine.

Methods I have tried to fix it:

(1)When I set all h > 0 and turn off WET_DRY, the model runs successfully. But I still want to use WET_DRY.

(2)Changing dt and Ndfast only delays the time of blow-up, for example,1.9s to 45s.

(3)Decreasing rx0 and rx1 does not work.

(4)I checked the value of u,v and w in the whole area before blow-up, it seems good, and no large value appeared.

I consider some other possible reasons:

(1) I import myself bathymetry into Gridbuilder but still use Gridbuilder's GSHHG coastline, maybe the value of h does not match mask (for exmample h<0 but mask=1, but dry place cannot make calculation) ? --------So if I import myself bathymetry, should I also input corresponding coastline (I don't have coastline data, I can only use matlab to generate it according to my bath data)
_
(2) The zeta of land area changed with time, which is weird.

Here is result.out, thanks!

Model Input Parameters: ROMS/TOMS version 3.7
Saturday - January 11, 2020 - 2:59:20 PM
-----------------------------------------------------------------------------

simpletest

Operating system : Linux
CPU/hardware : x86_64
Compiler system : gfortran
Compiler command : /usr/bin/mpif90
Compiler flags : -frepack-arrays -O3 -ffast-math -ftree-vectorize -ftree-loop-linear -funroll-loops -w -ffree-form -ffree-line-length-none -frecord-marker=4 -fconvert=big-endian

Input Script :

SVN Root URL : https:://myroms.org/svn/src
SVN Revision : 65M

Local Root : /home/ubuntu/app/COAWST_V3.2
Header Dir : /home/ubuntu/app/COAWST_V3.2/Projects/simpletest
Header file : simpletest.h
Analytical Dir: /home/ubuntu/app/COAWST_V3.2/Projects/simpletest

Resolution, Grid 01: 0176x0279x025, Parallel Nodes: 1, Tiling: 001x001


Physical Parameters, Grid: 01
=============================

29000 ntimes Number of timesteps for 3-D equations.
0.100 dt Timestep size (s) for 3-D equations.
20 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.
120 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).
1 nHIS Number of timesteps between the writing fields
into history file.
1 ntsAVG Starting timestep for the accumulation of output
time-averaged data.
360 nAVG Number of timesteps between the writing of
time-averaged data into averages file.
2.0000E-01 nl_tnu2(01) NLM Horizontal, harmonic mixing coefficient
(m2/s) for tracer 01: temp
2.0000E-01 nl_tnu2(02) NLM Horizontal, harmonic mixing coefficient
(m2/s) for tracer 02: salt
1.0000E-01 nl_visc2 NLM Horizontal, harmonic mixing coefficient
(m2/s) for momentum.
F LuvSponge Turning OFF sponge on horizontal momentum.
F LtracerSponge(01) Turning OFF sponge on tracer 01: temp
F LtracerSponge(02) Turning OFF sponge on tracer 02: salt
1.0000E-06 Akt_bak(01) Background vertical mixing coefficient (m2/s)
for tracer 01: temp
1.0000E-06 Akt_bak(02) Background vertical mixing coefficient (m2/s)
for tracer 02: salt
1.0000E-05 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.000 gls_p GLS stability exponent.
1.500 gls_m GLS turbulent kinetic energy exponent.
-1.000 gls_n GLS turbulent length scale exponent.
7.6000E-06 gls_Kmin GLS minimum value of turbulent kinetic energy.
1.0000E-12 gls_Pmin GLS minimum value of dissipation.
5.4770E-01 gls_cmu0 GLS stability coefficient.
1.4400E+00 gls_c1 GLS shear production coefficient.
1.9200E+00 gls_c2 GLS dissipation coefficient.
-4.0000E-01 gls_c3m GLS stable buoyancy production coefficient.
1.0000E+00 gls_c3p GLS unstable buoyancy production coefficient.
1.0000E+00 gls_sigk GLS constant Schmidt number for TKE.
1.3000E+00 gls_sigp GLS constant Schmidt number for PSI.
1400.000 charnok_alpha Charnok factor for Zos calculation.
0.500 zos_hsig_alpha Factor for Zos calculation using Hsig(Awave).
0.250 sz_alpha Factor for Wave dissipation surface tke flux .
100.000 crgban_cw Factor for Craig/Banner surface tke flux.
0.000 wec_alpha WEC factor for roller/breaking energy distribution.
3.0000E-04 rdrg Linear bottom drag coefficient (m/s).
2.5000E-02 rdrg2 Quadratic bottom drag coefficient.
2.0000E-02 Zob Bottom roughness (m).
2.0000E-02 Zos Surface roughness (m).
1.0000E-01 Dcrit Minimum depth for wetting and drying (m).
1 levsfrc Deepest level to apply surface stress as a
bodyforce.
1 levbfrc Shallowest level to apply bottom stress as a
bodyforce.
2 Vtransform S-coordinate transformation equation.
4 Vstretching S-coordinate stretching function.
6.0000E+00 theta_s S-coordinate surface control parameter.
4.0000E-01 theta_b S-coordinate bottom control parameter.
10.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 Tnudg_SSS Nudging/relaxation time scale (days)
for sea surface salinity.
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.
10.000 T0 Background potential temperature (C) constant.
35.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).
F LuvSrc Turning OFF momentum point Sources/Sinks.
F LwSrc Turning OFF volume influx point Sources/Sinks.
F LtracerSrc(01) Turning OFF point Sources/Sinks on tracer 01: temp
F LtracerSrc(02) Turning OFF point Sources/Sinks on tracer 02: salt
F LsshCLM Turning OFF processing of SSH climatology.
F Lm2CLM Turning OFF processing of 2D momentum climatology.
F Lm3CLM Turning OFF processing of 3D momentum climatology.
F LtracerCLM(01) Turning OFF processing of climatology tracer 01: temp
F LtracerCLM(02) Turning OFF processing of climatology tracer 02: salt
F LnudgeM2CLM Turning OFF nudging of 2D momentum climatology.
F LnudgeM3CLM Turning OFF nudging of 3D momentum climatology.
F LnudgeTCLM(01) Turning OFF nudging of climatology tracer 01: temp
F LnudgeTCLM(02) Turning OFF nudging of climatology tracer 02: salt
T Hout(idFsur) Write out free-surface.
T Hout(idUvel) Write out 3D U-momentum component.
T Hout(idVvel) Write out 3D V-momentum component.
T Hout(idWvel) Write out W-momentum component.


Output/Input Files:

Output Restart File: simpletest_rst.nc
Output History File: simpletest_his.nc
Output Averages File: simpletest_avg.nc
Input Grid File: /home/ubuntu/app/COAWST_V3.2/Projects/simpletest/simpletest_grid.nc
Input Nonlinear Initial File: /home/ubuntu/app/COAWST_V3.2/Projects/simpletest/simpletest_init.nc
Input Forcing File 01: /home/ubuntu/app/COAWST_V3.2/Projects/simpletest/simpletest_stress.nc

Tile partition information for Grid 01: 0176x0279x0025 tiling: 001x001

tile Istr Iend Jstr Jend Npts

Number of tracers: 2
0 1 176 1 279 1227600

Tile minimum and maximum fractional coordinates for Grid 01:
(interior points only)

tile Xmin Xmax Ymin Ymax grid

0 0.50 176.50 0.50 279.50 RHO-points

0 1.00 176.00 0.50 279.50 U-points

0 0.50 176.50 1.00 279.00 V-points

Maximum halo size in XI and ETA directions:

HaloSizeI(1) = 374
HaloSizeJ(1) = 578
TileSide(1) = 283
TileSize(1) = 51223


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

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

zeta 1 Chapman Imp Chapman Imp Chapman Imp Chapman Imp

ubar 1 Closed Flather Closed Flather

vbar 1 Closed Flather Closed Flather

u 1 Closed Gradient Closed Gradient

v 1 Closed Gradient Closed Gradient

temp 1 Closed Closed Closed Closed

salt 1 Closed Closed Closed Closed

tke 1 Closed Gradient Closed Gradient

Activated C-preprocessing Options:

PALU palu
ANA_BSFLUX Analytical kinematic bottom salinity flux.
ANA_BTFLUX Analytical kinematic bottom temperature flux.
ANA_FSOBC Analytical free-surface boundary conditions.
ANA_M2OBC Analytical 2D momentum boundary conditions.
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.
BODYFORCE Momentum stresses as body-forces.
CURVGRID Orthogonal curvilinear grid.
DJ_GRADPS Parabolic Splines density Jacobian (Shchepetkin, 2002).
DOUBLE_PRECISION Double precision arithmetic.
GLS_MIXING Generic Length-Scale turbulence closure.
KANTHA_CLAYSON Kantha and Clayson stability function formulation.
MASKING Land/Sea masking.
MIX_GEO_TS Mixing of tracers along geopotential surfaces.
MIX_S_UV Mixing of momentum along constant S-surfaces.
MPI MPI distributed-memory configuration.
NONLINEAR Nonlinear Model.
!NONLIN_EOS Linear Equation of State for seawater.
N2S2_HORAVG Horizontal smoothing of buoyancy and shear.
POWER_LAW Power-law shape time-averaging barotropic filter.
PROFILE Time profiling activated .
K_GSCHEME Third-order upstream advection of TKE fields.
RI_SPLINES Parabolic Spline Reconstruction for Richardson Number.
!RST_SINGLE Double precision fields in restart NetCDF file.
SALINITY Using salinity.
SOLVE3D Solving 3D Primitive Equations.
TS_U3HADVECTION Third-order upstream horizontal advection of tracers.
TS_C4VADVECTION Fourth-order centered vertical advection of tracers.
TS_DIF2 Harmonic mixing of tracers.
UV_ADV Advection of momentum.
UV_U3HADVECTION Third-order upstream horizontal advection of 3D momentum.
UV_C4VADVECTION Fourth-order centered vertical advection of momentum.
UV_LOGDRAG Logarithmic bottom stress.
UV_VIS2 Harmonic mixing of momentum.
VAR_RHO_2D Variable density barotropic mode.
WET_DRY Wetting and drying activated.

Process Information:

Node # 0 (pid= 13025) is active.

INITIAL: Configuring and initializing forward nonlinear model ...
*******

Vertical S-coordinate System, Grid 01:

level S-coord Cs-curve Z at hmin at hc half way at hmax

25 0.0000000 0.0000000 0.000 0.000 0.000 0.000
24 -0.0400000 -0.0001749 -1.192 -0.201 -0.459 -0.538
23 -0.0800000 -0.0007098 -2.368 -0.404 -1.061 -1.368
22 -0.1200000 -0.0016353 -3.526 -0.608 -1.816 -2.516
21 -0.1600000 -0.0030046 -4.665 -0.815 -2.747 -4.025
20 -0.2000000 -0.0048963 -5.780 -1.024 -3.883 -5.960
19 -0.2400000 -0.0074189 -6.866 -1.237 -5.269 -8.407
18 -0.2800000 -0.0107165 -7.918 -1.454 -6.960 -11.485
17 -0.3200000 -0.0149770 -8.926 -1.675 -9.031 -15.347
16 -0.3600000 -0.0204424 -9.880 -1.902 -11.577 -20.189
15 -0.4000000 -0.0274214 -10.766 -2.137 -14.721 -26.262
14 -0.4400000 -0.0363057 -11.566 -2.382 -18.616 -33.886
13 -0.4800000 -0.0475899 -12.258 -2.638 -23.458 -43.462
12 -0.5200000 -0.0618955 -12.815 -2.909 -29.492 -55.497
11 -0.5600000 -0.0799996 -13.200 -3.200 -37.024 -70.622
10 -0.6000000 -0.1028685 -13.371 -3.514 -46.436 -89.623
9 -0.6400000 -0.1316956 -13.274 -3.858 -58.198 -113.473
8 -0.6800000 -0.1679410 -12.843 -4.240 -72.887 -143.358
7 -0.7200000 -0.2133714 -11.998 -4.667 -91.199 -180.716
6 -0.7600000 -0.2700910 -10.646 -5.150 -113.964 -227.259
5 -0.8000000 -0.3405555 -8.675 -5.703 -142.152 -284.985
4 -0.8400000 -0.4275481 -5.960 -6.338 -176.860 -356.158
3 -0.8800000 -0.5340940 -2.366 -7.070 -219.282 -443.241
2 -0.9200000 -0.6632759 2.247 -7.916 -270.633 -548.740
1 -0.9600000 -0.8179050 8.006 -8.890 -332.022 -674.944
0 -1.0000000 -1.0000000 15.000 -10.000 -404.247 -823.494

Time Splitting Weights for Grid 01: ndtfast = 20 nfast = 29
==================================

Primary Secondary Accumulated to Current Step

1-0.0009651193358779 0.0500000000000000-0.0009651193358779 0.0500000000000000
2-0.0013488780126037 0.0500482559667939-0.0023139973484816 0.1000482559667939
3-0.0011514592651644 0.0501156998674241-0.0034654566136460 0.1501639558342180
4-0.0003735756740661 0.0501732728306823-0.0038390322877122 0.2003372286649003
5 0.0009829200513762 0.0501919516143856-0.0028561122363360 0.2505291802792859
6 0.0029141799764308 0.0501428056118168 0.0000580677400949 0.3006719858911027
7 0.0054132615310267 0.0499970966129953 0.0054713292711215 0.3506690825040980
8 0.0084687837865133 0.0497264335364439 0.0139401130576348 0.4003955160405419
9 0.0120633394191050 0.0493029943471183 0.0260034524767398 0.4496985103876601
10 0.0161716623600090 0.0486998273761630 0.0421751148367488 0.4983983377638231
11 0.0207585511322367 0.0478912442581626 0.0629336659689855 0.5462895820219857
12 0.0257765478740990 0.0468533167015507 0.0887102138430846 0.5931428987235364
13 0.0311633730493853 0.0455644893078458 0.1198735868924699 0.6387073880313822
14 0.0368391158442262 0.0440063206553765 0.1567127027366961 0.6827137086867586
15 0.0427031802506397 0.0421643648631652 0.1994158829873358 0.7248780735499238
16 0.0486309868367616 0.0400292058506332 0.2480468698240974 0.7649072794005570
17 0.0544704302037592 0.0375976565087951 0.3025173000278565 0.8025049359093521
18 0.0600380921294286 0.0348741349986072 0.3625553921572851 0.8373790709079593
19 0.0651152103984763 0.0318722303921357 0.4276706025557614 0.8692513013000950
20 0.0694434033194839 0.0286164698722119 0.4971140058752453 0.8978677711723069
21 0.0727201499285570 0.0251442997062377 0.5698341558038023 0.9230120708785446
22 0.0745940258796570 0.0215082922098099 0.6444281816834592 0.9445203630883545
23 0.0746596950216179 0.0177785909158270 0.7190878767050771 0.9622989540041815
24 0.0724526566618460 0.0140456061647461 0.7915405333669231 0.9763445601689277
25 0.0674437485167025 0.0104229733316538 0.8589842818836255 0.9867675335005816
26 0.0590334053485720 0.0070507859058187 0.9180176872321975 0.9938183194064003
27 0.0465456732896125 0.0040991156383901 0.9645633605218099 0.9979174350447905
28 0.0292219798521905 0.0017718319739095 0.9937853403740003 0.9996892670187000
29 0.0062146596259993 0.0003107329813000 0.9999999999999997 1.0000000000000000

ndtfast, nfast = 20 29 nfast/ndtfast = 1.45000

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

1.000000000000 1.060707743385 0.530353871693 1.000000000000 1.000000000000

Power filter parameters, Fgamma, gamma = 0.28400 0.14200

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

Minimum X-grid spacing, DXmin = 1.00256452E-01 km
Maximum X-grid spacing, DXmax = 1.00262565E-01 km
Minimum Y-grid spacing, DYmin = 9.95304436E-02 km
Maximum Y-grid spacing, DYmax = 9.95306247E-02 km
Minimum Z-grid spacing, DZmin = -1.63106062E+04 m
Maximum Z-grid spacing, DZmax = 5.81940781E+04 m

Minimum barotropic Courant Number = 4.04695013E-05
Maximum barotropic Courant Number = 6.36221117E-03
Maximum Coriolis Courant Number = 2.24399702E-07


NLM: GET_STATE - Read state initial conditions, t = 0 00:00:00
(Grid 01, File: palu_init.nc, Rec=0001, Index=1)
- free-surface
(Min = 0.00000000E+00 Max = 0.00000000E+00)
- vertically integrated u-momentum component
(Min = 0.00000000E+00 Max = 0.00000000E+00)
- vertically integrated v-momentum component
(Min = 0.00000000E+00 Max = 0.00000000E+00)
- u-momentum component
(Min = 0.00000000E+00 Max = 0.00000000E+00)
- v-momentum component
(Min = 0.00000000E+00 Max = 0.00000000E+00)
- potential temperature
(Min = 1.00000000E+01 Max = 1.00000000E+01)
- salinity
(Min = 3.50000000E+01 Max = 3.50000000E+01)
GET_2DFLD - surface u-momentum stress, t = 0 00:00:00
(Rec=0000001, Index=1, File: palu_stress.nc)
(Tmin= 0.0000 Tmax= 3001.0000)
(Min = 0.00000000E+00 Max = 0.00000000E+00)
GET_2DFLD - surface v-momentum stress, t = 0 00:00:00
(Rec=0000001, Index=1, File: palu_stress.nc)
(Tmin= 0.0000 Tmax= 3001.0000)
(Min = 0.00000000E+00 Max = 0.00000000E+00)

Basin information for Grid 01:

Maximum grid stiffness ratios: rx0 = 1.000000E-01 (Beckmann and Haidvogel)
rx1 = 4.833424E+00 (Haney)

Initial basin volumes: TotVolume = 3.1491821686E+10 m3
MinVolume = 1.3255664177E+01 m3
MaxVolume = 1.4784641003E+06 m3
Max/Min = 1.1153451691E+05

NL ROMS/TOMS: started time-stepping: (Grid: 01 TimeSteps: 00000001 - 00029000)

GET_2DFLD - surface u-momentum stress, t = 0 02:24:00
(Rec=0000002, Index=2, File: simpletest_stress.nc)
(Tmin= 0.0000 Tmax= 3001.0000)
(Min = 0.00000000E+00 Max = 0.00000000E+00)
GET_2DFLD - surface v-momentum stress, t = 0 02:24:00
(Rec=0000002, Index=2, File: simpletest_stress.nc)
(Tmin= 0.0000 Tmax= 3001.0000)
(Min = 0.00000000E+00 Max = 0.00000000E+00)

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 2.168775E+03 2.168775E+03 3.152385E+10
(000,000,00) 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00
DEF_HIS - creating history file, Grid 01: palu_his.nc
WRT_HIS - wrote history fields (Index=1,1) into time record = 0000001
DEF_AVG - creating average file, Grid 01: palu_avg.nc
1 0 00:00:00 6.758682E-35 2.168775E+03 2.168775E+03 3.152385E+10
(055,279,01) 2.956034E-20 2.295978E-20 0.000000E+00 6.972517E-17
WRT_HIS - wrote history fields (Index=2,2) into time record = 0000002
2 0 00:00:00 2.663834E-34 2.168775E+03 2.168775E+03 3.152385E+10
(053,278,25) 1.491527E-22 2.527561E-22 1.101529E-16 1.607954E-15
WRT_HIS - wrote history fields (Index=1,1) into time record = 0000003
3 0 00:00:00 5.868335E-34 2.168775E+03 2.168775E+03 3.152385E+10
(053,278,25) 5.375947E-22 9.461364E-22 2.172604E-16 4.785947E-15
WRT_HIS - wrote history fields (Index=2,2) into time record = 0000004
4 0 00:00:00 1.009185E-33 2.168775E+03 2.168775E+03 3.152385E+10
(055,189,25) 8.942073E-19 1.499638E-18 3.354848E-16 1.092499E-14
WRT_HIS - wrote history fields (Index=1,1) into time record = 0000005
5 0 00:00:00 1.500774E-33 2.168775E+03 2.168775E+03 3.152385E+10
(055,189,25) 9.227794E-19 1.559021E-18 5.960988E-16 1.806866E-14
WRT_HIS - wrote history fields (Index=2,2) into time record = 0000006
6 0 00:00:00 2.036808E-33 2.168775E+03 2.168775E+03 3.152385E+10
(055,189,25) 9.549458E-19 1.633358E-18 8.576736E-16 2.514534E-14
WRT_HIS - wrote history fields (Index=1,1) into time record = 0000007
7 0 00:00:00 2.613172E-33 2.168775E+03 2.168775E+03 3.152385E+10
(055,189,25) 9.910394E-19 1.725476E-18 1.119126E-15 3.213089E-14
WRT_HIS - wrote history fields (Index=2,2) into time record = 0000008
8 0 00:00:00 3.300044E-33 2.168775E+03 2.168775E+03 3.152385E+10
(055,189,25) 1.031782E-18 1.836636E-18 1.380427E-15 3.901778E-14
WRT_HIS - wrote history fields (Index=1,1) into time record = 0000009
9 0 00:00:00 4.200563E-33 2.168775E+03 2.168775E+03 3.152385E+10
(055,189,25) 1.077752E-18 1.967280E-18 1.641572E-15 4.580510E-14
WRT_HIS - wrote history fields (Index=2,2) into time record = 0000010
10 0 00:00:00 5.564856E-33 2.168775E+03 2.168775E+03 3.152385E+10
(055,189,25) 1.129313E-18 2.117548E-18 1.902560E-15 5.249325E-14
WRT_HIS - wrote history fields (Index=1,1) into time record = 0000011
11 0 00:00:01 7.152720E-33 2.168775E+03 2.168775E+03 3.152385E+10
(055,189,25) 1.186620E-18 2.287415E-18 2.163388E-15 5.908300E-14
WRT_HIS - wrote history fields (Index=2,2) into time record = 0000012
12 0 00:00:01 9.194294E-33 2.168775E+03 2.168775E+03 3.152385E+10
(055,189,25) 1.249703E-18 2.476755E-18 2.424058E-15 6.557535E-14
WRT_HIS - wrote history fields (Index=1,1) into time record = 0000013
13 0 00:00:01 1.197251E-32 2.168775E+03 2.168775E+03 3.152385E+10
(055,189,25) 1.318613E-18 2.685478E-18 2.684569E-15 7.197098E-14
WRT_HIS - wrote history fields (Index=2,2) into time record = 0000014
14 0 00:00:01 1.512995E-32 2.168775E+03 2.168775E+03 3.152385E+10
(055,189,25) 1.393348E-18 2.913450E-18 2.944922E-15 7.827065E-14
WRT_HIS - wrote history fields (Index=1,1) into time record = 0000015
15 0 00:00:01 1.864678E-32 2.168775E+03 2.168775E+03 3.152385E+10
(055,189,25) 1.473913E-18 3.160544E-18 3.205115E-15 8.447506E-14
WRT_HIS - wrote history fields (Index=2,2) into time record = 0000016
16 0 00:00:01 2.289543E-32 2.168775E+03 2.168775E+03 3.152385E+10
(055,189,25) 1.560307E-18 3.426630E-18 3.465147E-15 9.058493E-14
WRT_HIS - wrote history fields (Index=1,1) into time record = 0000017
17 0 00:00:01 2.740619E-32 2.168775E+03 2.168775E+03 3.152385E+10
(055,189,25) 1.652521E-18 3.711577E-18 3.725015E-15 9.660098E-14
WRT_HIS - wrote history fields (Index=2,2) into time record = 0000018
18 0 00:00:01 NaN NaN NaN 3.152385E+10
(055,189,25) 1.750576E-18 4.015282E-18 3.984722E-15 1.025239E-13

Blowing-up: Saving latest model state into RESTART file

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

Elapsed CPU time (seconds):

Node # 0 CPU: 19.912
Total: 19.912

Nonlinear ocean model elapsed time profile, Grid: 01

Allocation and array initialization .............. 0.124 ( 0.6211 %)
Ocean state initialization ....................... 0.015 ( 0.0760 %)
Reading of input data ............................ 0.008 ( 0.0379 %)
Processing of input data ......................... 0.005 ( 0.0255 %)
Processing of output time averaged data .......... 0.010 ( 0.0480 %)
Computation of vertical boundary conditions ...... 0.029 ( 0.1456 %)
Computation of global information integrals ...... 0.130 ( 0.6551 %)
Writing of output data ........................... 2.399 (12.0497 %)
Model 2D kernel .................................. 3.545 (17.8024 %)
2D/3D coupling, vertical metrics ................. 0.141 ( 0.7105 %)
Omega vertical velocity .......................... 0.104 ( 0.5225 %)
Equation of state for seawater ................... 0.119 ( 0.5991 %)
GLS vertical mixing parameterization ............. 8.878 (44.5861 %)
3D equations right-side terms .................... 0.645 ( 3.2384 %)
3D equations predictor step ...................... 0.666 ( 3.3463 %)
Pressure gradient ................................ 0.381 ( 1.9151 %)
Harmonic mixing of tracers, geopotentials ........ 0.819 ( 4.1122 %)
Harmonic stress tensor, S-surfaces ............... 0.315 ( 1.5815 %)
Corrector time-step for 3D momentum .............. 0.519 ( 2.6069 %)
Corrector time-step for tracers .................. 0.570 ( 2.8615 %)
Total: 19.423 97.5414

Nonlinear model message Passage profile, Grid: 01

Message Passage: 2D halo exchanges ............... 0.010 ( 0.0516 %)
Message Passage: 3D halo exchanges ............... 0.001 ( 0.0033 %)
Message Passage: 4D halo exchanges ............... 0.000 ( 0.0007 %)
Message Passage: data broadcast .................. 0.001 ( 0.0060 %)
Message Passage: data reduction .................. 0.000 ( 0.0008 %)
Message Passage: data gathering .................. 0.289 ( 1.4521 %)
Message Passage: data scattering.................. 0.022 ( 0.1105 %)
Message Passage: boundary data gathering ......... 0.000 ( 0.0005 %)
Total: 0.324 1.6255

All percentages are with respect to total time = 19.912


ROMS/TOMS - Output NetCDF summary for Grid 01:
number of time records written in HISTORY file = 00000018
number of time records written in RESTART file = 00000001

Analytical header files used:

ROMS/Functionals/ana_btflux.h
/home/ubuntu/app/COAWST_V3.2/Projects/simpletest/ana_fsobc.h
/home/ubuntu/app/COAWST_V3.2/Projects/simpletest/ana_m2obc.h
ROMS/Functionals/ana_stflux.h

ROMS/TOMS: DONE... Saturday - January 11, 2020 - 2:59:43 PM

[zhangtianhao]

What I do these 2 days:

(1)My COAWST version does not need to define limit_bstress.

(2)I found that many people's model blowed up because of large T,S and velocities. So I checked T and S, they seemed also fine before blowing up, the same as uvw.

(3)My model blowed up with kinetic_enrg, poten_enrg and total_enrg suddenly becoming NaN. I don't understand why velocities are small but energy is large.

(4)I activated WET_DRY and my hmin was ~-100. I turn hmin to -15, result did not change. But when I turn it to -5, it can run without blowing up.

So I kept to turn it to 0, it can also run.(I use gridbuilder to generate grid)

Noticed that rx0 and rx1 are always small(I use gridbuilder to check them).

I can't explain this phenomenen.


Hoty

2020-01-14

[jcwarner]

-changing hmin to zero is basically removing all the 'land' parts of your grid, so that removed your wet dry problem.
- an hmin of -100 means that you have a small hill about 100 m above sea level. that might be troubles for roms. so maybe you could set landmask=0 for cells that are above, lets say, -10 m elevation.
so landmask (h<-10) = 0;
- for wet_dry, zeta is draped over the bathy.
- what value of Dcrit did you use? this grid seems rather course for wetdry, but it should work.
- use ncview or similar to step through the netcdf output file. where is the model blowing up?
-j

[zhangtianhao]

-changing hmin to zero is basically removing all the 'land' parts of your grid, so that removed your wet dry problem.
- an hmin of -100 means that you have a small hill about 100 m above sea level. that might be troubles for roms. so maybe you could set landmask=0 for cells that are above, lets say, -10 m elevation.
so landmask (h<-10) = 0;
- for wet_dry, zeta is draped over the bathy.
- what value of Dcrit did you use? this grid seems rather course for wetdry, but it should work.
- use ncview or similar to step through the netcdf output file. where is the model blowing up?
-j

Thanks for your reply,
(1) I changed hmin to zero, it removed my wet dry problem. But when I changed hmin to -5, it can also removed my wet dry problem.
But, when I changed hmin to -10m, model blowed up. I wonder the prople range of hmin.
(2) I activated MASK and WET_DRY and all cells which h<0 have landmask value of 0. I wonder the meaning of set landmask (h<-10) = 0?
(3)Thanks, so in landmask cell, zeta is draped over the bathy; but in oceanmask cell, zeta seemed to be calculated at the base of 0 elevation.
(4)I use Dcrit=0.1m. You mean I should refine my grid? Now horiziontal resolution is 100m.
(5)Model blowed up at oceanmask place, I checked uvts, all well.

Thanks, you help me learn a lot.

Hoty

[kate]

You mean it died inside the land mask? Do you know where in the code this happened?

Dcrit is the depth of the water in the "dry" places. You can never remove that last bit of water from a cell. That's what keeps things stable (if all is working right).

[johnluick]

The author of Gridbuilder is on leave this week. Maybe he will have a suggestion, though it does not appear to be a Gridbuilder problem as such.

[jcwarner]

(4)I use Dcrit=0.1m. You mean I should refine my grid? Now horiziontal resolution is 100m.
this is clear that you are not fully aware of how to use wet dry. it is simple. can you get our paper in Comp and Geosci?

https://darchive.mblwhoilibrary.org/bit ... sequence=1

[zhangtianhao]

Thanks for suggestions from Kate, johnluick and jcwarner, that's really useful.

Firstly, I read this papaer https://darchive.mblwhoilibrary.org/bit ... sequence=1. I found I misunderstood some parameters.

I confused land/sea mask with wer/dry mask.

Now I have a question. In previous case, I set landmask(h<0)=0 using Gridbuilder. I consider it's a mistake, according to jcwarner's suggestion( set landmask(h>-10)=1 )and I shoud set intertidal place which -10<h<0 as landmask=1.

But I found when using Gridbuilder, I can't set sea mask on land, when I change the mask, the depth changes to adapt it.

For example, I set sea mask on land(-10<h<0), but depth of these cells changed from negative to positive values automatically.

The consequence is h<0~~~~ landmask=0 or h>0~~~~ landmask=1 .

So I generate a normal grid by gridbuilder firstly, and then use matlab to modify the mask ( mask_rho(h>-10)=0) , but after modifying mask, I found rx1 became very large. And model blowed up...

Thanks again and hope anyone can help me out

Hoty
2020-01-16

[CharlesJames]

I think an old version of GridBuilder had an issue with resetting depths with land masks to 0m when exporting to a grid nc file regardless of users intent but that was fixed a while ago; you can set the min depth to anything on export. I am working with version 1.2 of GridBuilder and I am able to create a wet cell on land with depth = -5m and export it and it everything seems to be right in the netcdf file. If you modify the grid shape or reset the mask or depths, then the mask and/or depths are recomputed and any changes are lost (you can always undo though) .

The algorithm that determines depth from the Set Min and Set Max is fairly straight forward, in the routine setGridDepths.m you will find the following lines:

depths(depths>maxDepth)=maxDepth;
depths(depths<minDepth)=minDepth;

Since the modifications are always applied to the current GridBuilder depths if you previously 'Set Min' to -5m then setting it to something lower -15m it will not change any values as there are now no depths <-15 to select, the same logic would apply for changes to maximum depth. The option to reset the depths is always there but you will lose any filtering. I can see how a user might expect changing 'Set Min' to force the shallowest depth to be that value so perhaps modifying the code to add the following would be more useful.

indMax=(depths==max(depths(:)))|(depths>maxDepth);
indMin=(depths==min(depths(:)))|(depths<minDepth);
depths(indMax)=maxDepth;
depths(indMin)=minDepth;

I'm testing it out now and it looks like it doesn't create any adverse affects as long as you keep in mind that each time you edit the limits any information about topography outside the limits is lost until you hit Reset and start over, but that is true for both approaches. If this method doesn't cause any other unwanted side effects I may retain it for future releases.

[zhangtianhao]

I think an old version of GridBuilder had an issue with resetting depths with land masks to 0m when exporting to a grid nc file regardless of users intent but that was fixed a while ago; you can set the min depth to anything on export. I am working with version 1.2 of GridBuilder and I am able to create a wet cell on land with depth = -5m and export it and it everything seems to be right in the netcdf file. If you modify the grid shape or reset the mask or depths, then the mask and/or depths are recomputed and any changes are lost (you can always undo though) .

The algorithm that determines depth from the Set Min and Set Max is fairly straight forward, in the routine setGridDepths.m you will find the following lines:

depths(depths>maxDepth)=maxDepth;
depths(depths<minDepth)=minDepth;

Since the modifications are always applied to the current GridBuilder depths if you previously 'Set Min' to -5m then setting it to something lower -15m it will not change any values as there are now no depths <-15 to select, the same logic would apply for changes to maximum depth. The option to reset the depths is always there but you will lose any filtering. I can see how a user might expect changing 'Set Min' to force the shallowest depth to be that value so perhaps modifying the code to add the following would be more useful.

indMax=(depths==max(depths(:)))|(depths>maxDepth);
indMin=(depths==min(depths(:)))|(depths<minDepth);
depths(indMax)=maxDepth;
depths(indMin)=minDepth;

I'm testing it out now and it looks like it doesn't create any adverse affects as long as you keep in mind that each time you edit the limits any information about topography outside the limits is lost until you hit Reset and start over, but that is true for both approaches. If this method doesn't cause any other unwanted side effects I may retain it for future releases.

Dear CharlesJames,

Thanks for your reply. I am using Gridbuilder 1.2 now by following the guide https://austides.com/wp-content/uploads ... -v0.99.pdf. But when I import bathymetry data, the mask will change according to the bathymetry.

For example, I import bathymetry data of a coastal area and set hmin as -100m, the h data seems right, but mask(h>0) always equals 1 and mask(h<0) always equals 0.

I wonder that how to achieve this：for a coastal area, set hmin as -100m, mask(h>-10)=1 and mask(h<-10)=0.

Thanks,

Hoty

[CharlesJames]

Are you using GSHHG coastlines to generate the mask? Under the Mask Settings menu choose "Use GSHHG Coastlines (More Accurate)" - default is "Use Topography (Faster)" which is very useful for quickly placing a grid and getting the resolution correct. Once you are happy with the basic grid I suggest that you use GSHHG (or your own imported coastline) to get a much better mask and better editing. Once the mask is correct then work on the topography.

I have noticed one inconsistent behaviour that I will fix up, when the Mask is computed by topography it goes back to the raw bathymetry data (eTopo or imported) to estimate the coastline rather than any modifications that the user might have made (i.e. new minimum depth). Since the grid doesn't change we should just use the current bathymetry data without recomputing it. The bathymetry computed for the fast-masking routines isn't stored over the existing bathymetry so it shouldn't change the modifications being exported but does lead to some odd results during the mask editing. I'll fix this up in the next revision.

[zhangtianhao]

Are you using GSHHG coastlines to generate the mask? Under the Mask Settings menu choose "Use GSHHG Coastlines (More Accurate)" - default is "Use Topography (Faster)" which is very useful for quickly placing a grid and getting the resolution correct. Once you are happy with the basic grid I suggest that you use GSHHG (or your own imported coastline) to get a much better mask and better editing. Once the mask is correct then work on the topography.

I have noticed one inconsistent behaviour that I will fix up, when the Mask is computed by topography it goes back to the raw bathymetry data (eTopo or imported) to estimate the coastline rather than any modifications that the user might have made (i.e. new minimum depth). Since the grid doesn't change we should just use the current bathymetry data without recomputing it. The bathymetry computed for the fast-masking routines isn't stored over the existing bathymetry so it shouldn't change the modifications being exported but does lead to some odd results during the mask editing. I'll fix this up in the next revision.

Thanks a lot! Yes, I am using my own coastline.

I found the issue is that, just as what you said above, the bathymetry is always recomputed to adapt the coastline. And basic bathymetry in the coastline is set to 0 automatically. But I want it set to -10 but not 0.

Here is my step to generate a grid:

1. Firstly I choose the grid area and then modify my bathymetry(by filter) and make rx0 and rx1 small enough.

2. Secondly, I import my own coastline, and generate a better mask. (My coastline is -10m contour line) But after I import coastline, I found the rx0 and rx1 become so large again. So I have to recompute it by positive filter Once I recompute it, depth of all wet cell becomes positive automatically(mask(h>0)=1), all dry cell becomes negative. That means the bathymetry recomputed with h of coastline area equals 0.

I wonder how to solve this?

Thanks agian!

Hoty

2020-02-13

[CharlesJames]

I think the problem is the bathymetric filtering. The positive and negative adjustment filters only work on wet cells and rx0 is only calculated on wet cells - if you do filtering first then alter the mask it may tag previously dry cells as wet - the topography of these previously dry cells has not been filtered so may create a gradient jump with adjacent wet cells causing the rx0 value to jump up - the topography shouldn't have changed but the new wet cells have unmasked topography which are poorly conditioned. along the coastline. Applying the positive (and negative) adjustment algorithms shouldn't change cells that are already ok, but they will probably make previously dry (h<0) cells wet (h>0) in order to satisfy the rx0 limit. One solution, if you are happy using your coastline file to define the mask, is to create an x,y,z bathymetry file then raise sea-level so more of the bathymetry is wet to start with, import this bathymetry file and your coastline file, create the mask with the coastline file and then filter the bathymetry file. You'd have to post-process the grid file to remove the sea-level offset but it might give you something closer to what you want.

This has all given me an idea which I will try and implement in the next release, I can add a sea-level offset that can be applied as part of the bathymetric editing which would set the relative 0m level - this would also affect the fast topography algorithm so you could play with various high or low sea-level scenarios. It would also make your problem a bit easier to solve!

[zhangtianhao]

I think the problem is the bathymetric filtering. The positive and negative adjustment filters only work on wet cells and rx0 is only calculated on wet cells - if you do filtering first then alter the mask it may tag previously dry cells as wet - the topography of these previously dry cells has not been filtered so may create a gradient jump with adjacent wet cells causing the rx0 value to jump up - the topography shouldn't have changed but the new wet cells have unmasked topography which are poorly conditioned. along the coastline. Applying the positive (and negative) adjustment algorithms shouldn't change cells that are already ok, but they will probably make previously dry (h<0) cells wet (h>0) in order to satisfy the rx0 limit. One solution, if you are happy using your coastline file to define the mask, is to create an x,y,z bathymetry file then raise sea-level so more of the bathymetry is wet to start with, import this bathymetry file and your coastline file, create the mask with the coastline file and then filter the bathymetry file. You'd have to post-process the grid file to remove the sea-level offset but it might give you something closer to what you want.

This has all given me an idea which I will try and implement in the next release, I can add a sea-level offset that can be applied as part of the bathymetric editing which would set the relative 0m level - this would also affect the fast topography algorithm so you could play with various high or low sea-level scenarios. It would also make your problem a bit easier to solve!

Thanks for your helpful suggestion,

I think add a sea-level offset is a good way to solve this problem and I will take a try.

Hoty

[zhangtianhao]

I think the problem is the bathymetric filtering. The positive and negative adjustment filters only work on wet cells and rx0 is only calculated on wet cells - if you do filtering first then alter the mask it may tag previously dry cells as wet - the topography of these previously dry cells has not been filtered so may create a gradient jump with adjacent wet cells causing the rx0 value to jump up - the topography shouldn't have changed but the new wet cells have unmasked topography which are poorly conditioned. along the coastline. Applying the positive (and negative) adjustment algorithms shouldn't change cells that are already ok, but they will probably make previously dry (h<0) cells wet (h>0) in order to satisfy the rx0 limit. One solution, if you are happy using your coastline file to define the mask, is to create an x,y,z bathymetry file then raise sea-level so more of the bathymetry is wet to start with, import this bathymetry file and your coastline file, create the mask with the coastline file and then filter the bathymetry file. You'd have to post-process the grid file to remove the sea-level offset but it might give you something closer to what you want.

This has all given me an idea which I will try and implement in the next release, I can add a sea-level offset that can be applied as part of the bathymetric editing which would set the relative 0m level - this would also affect the fast topography algorithm so you could play with various high or low sea-level scenarios. It would also make your problem a bit easier to solve!

Sorry to ask that, but

I have tried this way, but I failed.

Here is my step:

1.Raised sea-level: add 15 meters to previous depth in bathymetry file and generate corresponding coastline file.

2.Import bathymetry file coastline file into Gridbuilder, create the mask with the coastline file and then filter the bathymetry.

3.Export grid.nc file which rx0 and rx1 are fine from Gridbuilder

4.Remove the sea-level offset from the grid.nc file: remove 15 meters to previous h(1:lon_rho,1:lat_rho) and change hmin and hmax, I don't change all of other variables;

But when I used this grid.nc file to run my case, ROMS told me rx0 and rx1 is too large.

Here is my ROMS result and here is the matlab scipt which use to remove the sea-level offset from the grid.nc

[zhangtianhao]

clear;

modis3='C:\Users\xzz\Desktop\ROMS.nc';
ncdisp(modis3);
ncid3=netcdf.open(modis3);

%1) Enter name of netcdf grid file to be created.
grid_file='ROMS_new.nc';

%2) Obtain grid information.
LP = 178;
MP = 281;
L = LP-1;
M = MP-1;
xi_psi = L;
xi_rho = LP;
xi_u = L;
xi_v = LP;
eta_psi = M;
eta_rho = MP;
eta_u = MP;
eta_v = M;
spherical = 'T';

%3)
%
h(1:xi_rho,1:eta_rho) = ncread(modis3,'h')-15;
depthmin = min(min(h));
depthmax = max(max(h));
f(1:xi_rho,1:eta_rho) = ncread(modis3,'f');
mask_rho(1:xi_rho,1:eta_rho) = ncread(modis3,'mask_rho');
mask_psi(1:xi_psi,1:eta_psi) = ncread(modis3,'mask_psi');
mask_u(1:xi_u,1:eta_u) = ncread(modis3,'mask_u');
mask_v(1:xi_v,1:eta_v) = ncread(modis3,'mask_v');
lon_psi(1:xi_psi,1:eta_psi) = ncread(modis3,'lon_psi');
lat_psi(1:xi_psi,1:eta_psi) = ncread(modis3,'lat_psi');
lon_rho(1:xi_rho,1:eta_rho) = ncread(modis3,'lon_rho');
lat_rho(1:xi_rho,1:eta_rho) = ncread(modis3,'lat_rho');
lon_u(1:xi_u,1:eta_u) = ncread(modis3,'lon_u');
lat_u(1:xi_u,1:eta_u) = ncread(modis3,'lat_u');
lon_v(1:xi_v,1:eta_v) = ncread(modis3,'lon_v');
lat_v(1:xi_v,1:eta_v) = ncread(modis3,'lat_v');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% END of USER INPUT %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%create grid file
create_roms_netcdf_grid_file(grid_file,LP,MP)

%now write the data from the arrays to the netcdf file
disp(' ## Filling Variables in netcdf file with data...')


ncwrite(grid_file,'spherical',spherical);
ncwrite(grid_file,'depthmin',depthmin);
ncwrite(grid_file,'depthmax',depthmax);
ncwrite(grid_file,'f',f);
ncwrite(grid_file,'h',h);

ncwrite(grid_file,'mask_rho',mask_rho);
ncwrite(grid_file,'mask_psi',mask_psi);
ncwrite(grid_file,'mask_u',mask_u);
ncwrite(grid_file,'mask_v',mask_v);
ncwrite(grid_file,'lon_psi',lon_psi);
ncwrite(grid_file,'lat_psi',lat_psi);
ncwrite(grid_file,'lon_rho',lon_rho);
ncwrite(grid_file,'lat_rho',lat_rho);
ncwrite(grid_file,'lon_u',lon_u);
ncwrite(grid_file,'lat_u',lat_u);
ncwrite(grid_file,'lon_v',lon_v);
ncwrite(grid_file,'lat_v',lat_v);


%close file
disp(['created ', grid_file])

[zhangtianhao]

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

Minimum X-grid spacing, DXmin = 1.00000000E+17 km
Maximum X-grid spacing, DXmax = Infinity km
Minimum Y-grid spacing, DYmin = 1.00000000E+17 km
Maximum Y-grid spacing, DYmax = Infinity km
Minimum Z-grid spacing, DZmin = -5.47066471E+03 m
Maximum Z-grid spacing, DZmax = 3.07163748E+03 m

Minimum barotropic Courant Number = 0.00000000E+00
Maximum barotropic Courant Number = 0.00000000E+00
Maximum Coriolis Courant Number = 2.29183476E-07


NLM: GET_STATE - Read state initial conditions, t = 0 00:00:00
(Grid 01, File: init.nc, Rec=0001, Index=1)
- free-surface
(Min = 0.00000000E+00 Max = 0.00000000E+00)
- vertically integrated u-momentum component
(Min = 0.00000000E+00 Max = 0.00000000E+00)
- vertically integrated v-momentum component
(Min = 0.00000000E+00 Max = 0.00000000E+00)
- u-momentum component
(Min = 0.00000000E+00 Max = 0.00000000E+00)
- v-momentum component
(Min = 0.00000000E+00 Max = 0.00000000E+00)
- potential temperature
(Min = 1.00000000E+01 Max = 1.00000000E+01)
- salinity
(Min = 3.50000000E+01 Max = 3.50000000E+01)

Basin information for Grid 01:

Maximum grid stiffness ratios: rx0 = 3.623208E+04 (Beckmann and Haidvogel)
rx1 = 8.566216E+05 (Haney)

Initial basin volumes: TotVolume = NaN m3
MinVolume = -Infinity m3
MaxVolume = Infinity m3
Max/Min = NaN

[CharlesJames]

rx0 can become very large if the depth changes sign within the unmasked part of the grid - I think the topography might not match up quite with your coastline after smoothing and raising it, This may have caused land to be exposed near the coast in unmasked cells which is not good and may be causing the problems you are seeing.

Looking back at previous posts, perhaps I'm misunderstanding your problem, it looks like you just want to make topography behind the land mask equal to -10m but have normal smoothed topography everywhere there are wet cells. This would essentially treat the land/sea boundary as a cliff and wetting and drying would only occur in cells that were originally wet. Is that right?

GridBuilder should be able to do create a suitable grid and topography in the wet cells but doesn't have a feature for setting dry cell depths. If I made the max and min depth edits apply only to the unmasked region then you could invert the mask to set the land values directly (by setting max depth to -10), there might be merit to doing that.

In the absence of that you would have use Matlab to replace depths where mask==0 with -10 after the grid is created.

[zhangtianhao]

rx0 can become very large if the depth changes sign within the unmasked part of the grid - I think the topography might not match up quite with your coastline after smoothing and raising it, This may have caused land to be exposed near the coast in unmasked cells which is not good and may be causing the problems you are seeing.

Looking back at previous posts, perhaps I'm misunderstanding your problem, it looks like you just want to make topography behind the land mask equal to -10m but have normal smoothed topography everywhere there are wet cells. This would essential
ly treat the land/sea boundary as a cliff and wetting and drying would only occur in cells that were originally wet. Is that right?

GridBuilder should be able to do create a suitable grid and topography in the wet cells but doesn't have a feature for setting dry cell depths. If I made the max and min depth edits apply only to the unmasked region then you could invert the mask to set the land values directly (by setting max depth to -10), there might be merit to doing that.
In the absence of that you would have use Matlab to replace depths where mask==0 with -10 after the grid is created.

Thank you so much for your patient answer.

I am so sorry that I didn't make myself clear . In fact, I only want to simulate a tidal area, like the picture shows below.

I use topograpy in the land area(h<0). But if the land/sea mask of land area equals 0, the tide will never rise into the beach (or the land). So I want to change the beach mask from 0 to 1. So the key point is that I want to set mask of a h<0 cell (beach cell) euqals 1.

So your mean :

I should set depth of beach area even land area >0 (for example, add 10 meters to all cells) and set initial zeta(1:xi_rho,1:eta_rho) =-10m in initial.nc file, then sea level=-10m and the depth of land area is positive. So Gridbuilder and ROMS can deal with it？

Thanks

Hoty

[CharlesJames]

OK, I think I get it:
So the steps would be (following what you said):

Using GridBuilder:
Add enough depth to your bathymetry so that everything that you expect to get wet is under water and load it into GridBuilder
Use the new 0 contour as your coast line (let GridBuilder use topography to set mask not coastlines)
Make any adjustments to the mask to remove unwanted cells
Now filter the bathymetry to get your rx0 right and export the grid.

Using ROMS:
Offset your inital zeta by the amount necessary to establish the correct mean sea level for the model run
If you use zeta to directly force tides make sure you offset that as well in your forcing.
set WET_DRY and MASKING and build your ROMS
set suitable DCRIT to determine minimum depth for dry cell to turn on

The model zero is now set at a point higher that you expect zeta to reach and wetting and drying should occur correctly in the intertidal zone.

Because all the bathymetry in the unmasked region is positive the grid stiffness shouldn't be too high now.

That makes sense to me - I think it should work.

[zhangtianhao]

OK, I think I get it:
So the steps would be (following what you said):

Using GridBuilder:
Add enough depth to your bathymetry so that everything that you expect to get wet is under water and load it into GridBuilder
Use the new 0 contour as your coast line (let GridBuilder use topography to set mask not coastlines)
Make any adjustments to the mask to remove unwanted cells
Now filter the bathymetry to get your rx0 right and export the grid.

Using ROMS:
Offset your inital zeta by the amount necessary to establish the correct mean sea level for the model run
If you use zeta to directly force tides make sure you offset that as well in your forcing.
set WET_DRY and MASKING and build your ROMS
set suitable DCRIT to determine minimum depth for dry cell to turn on

The model zero is now set at a point higher that you expect zeta to reach and wetting and drying should occur correctly in the intertidal zone.

Because all the bathymetry in the unmasked region is positive the grid stiffness shouldn't be too high now.

That makes sense to me - I think it should work.

Thank you so much,

I will have a try and I wonder that when using your Gridbuilder, which paper or which website should we cite?

Thanks again and look forward to your new version

[wilkin]

Keep in mind that the rx0 value is guidance only with respect to bathymetry steepness than might introduce truncation errors in the baroclinic pressure gradient. In shallow tidal flats that wet and dry it is unlikely you will see sufficiently large vertical stratification to introduce baroclinicity in the pressure field. The barotropic pressure gradient is always computed accurately from the free surface slope.

I would not obsess about rx0 in this application in the wetting/drying area. Filtering to remove length scales of a few dx in h is probably all you need to do.