Changeset 10499 for NEMO/trunk/tests
- Timestamp:
- 2019-01-10T16:12:24+01:00 (6 years ago)
- Location:
- NEMO/trunk/tests/WAD/MY_DOCS
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/tests/WAD/MY_DOCS/Namelists/nam_wad
r9024 r10499 11 11 rn_wdld = 2.5 ! Land elevation below which wetting/drying is allowed 12 12 nn_wdit = 20 ! Max iterations for W/D limiter 13 rn_wd_sbcdep = 5.0 ! Depth at which to taper sbc fluxes 14 rn_wd_sbcfra = 0.999 ! Fraction of SBC fluxes at taper depth (Must be <1) 13 15 / -
NEMO/trunk/tests/WAD/MY_DOCS/WAD_doc.tex
r9024 r10499 25 25 % Wetting and drying 26 26 % ================================================================ 27 \section{Wetting and drying }28 \label{DYN_wetdry}29 30 There are two main options for wetting and drying code (wd):31 (a) an iterative limiter (il) and (b) a directional limiter (dl).32 The directional limiter is based on the scheme developed by \cite{WarnerEtal13} for ROMS33 which was in turn based on ideas developed for POM by \cite{Oey06}. The iterative34 limiter is a new scheme. The iterative limiter is activated by setting $\mathrm{ln\_wd\_il} = \mathrm{.true.}$35 and $\mathrm{ln\_wd\_dl} = \mathrm{.false.}$. The directional limiter is activated36 by setting $\mathrm{ln\_wd\_dl} = \mathrm{.true.}$ and $\mathrm{ln\_wd\_il} = \mathrm{.false.}$.37 38 \namdisplay{nam_wad}39 40 The following terminology is used. The depth of the topography (positive downwards)41 at each $(i,j)$ point is the quantity stored in array $\mathrm{ht\_wd}$ in the NEMO code.42 The height of the free surface (positive upwards) is denoted by $ \mathrm{ssh}$. Given the sign43 conventions used, the water depth, $h$, is the height of the free surface plus the depth of the44 topography (i.e. $\mathrm{ssh} + \mathrm{ht\_wd}$).45 46 Both wd schemes take all points in the domain below a land elevation of $\mathrm{rn\_wdld}$ to be47 covered by water. They require the topography specified with a model48 configuration to have negative depths at points where the land is higher than the49 topography's reference sea-level. The vertical grid in NEMO is normally computed relative to an50 initial state with zero sea surface height elevation.51 The user can choose to compute the vertical grid and heights in the model relative to52 a non-zero reference height for the free surface. This choice affects the calculation of the metrics and depths53 (i.e. the $\mathrm{e3t\_0, ht\_0}$ etc. arrays).54 55 Points where the water depth is less than $\mathrm{rn\_wdmin1}$ are interpreted as ``dry''.56 $\mathrm{rn\_wdmin1}$ is usually chosen to be of order $0.05$m but extreme topographies57 with very steep slopes require larger values for normal choices of time-step.58 59 Both versions of the code have been tested in six test cases provided in the WAD\_TEST\_CASES configuration60 and in ``realistic'' configurations covering parts of the north-west European shelf.61 All these configurations have used pure sigma coordinates. It is expected that62 the wetting and drying code will work in domains with more general s-coordinates provided63 the coordinates are pure sigma in the region where wetting and drying actually occurs.64 65 The next sub-section descrbies the directional limiter and the following sub-section the iterative limiter.66 The final sub-section covers some additional considerations that are relevant to both schemes.67 68 %-----------------------------------------------------------------------------------------69 % Iterative limiters70 %-----------------------------------------------------------------------------------------71 \subsection [Directional limiter (\textit{wet\_dry})]72 {Directional limiter (\mdl{wet\_dry})}73 \label{DYN_wd_directional_limiter}74 75 The principal idea of the directional limiter is that76 water should not be allowed to flow out of a dry tracer cell (i.e. one whose water depth is less than rn\_wdmin1).77 78 All the changes associated with this option are made to the barotropic solver for the non-linear79 free surface code within dynspg\_ts.80 On each barotropic sub-step the scheme determines the direction of the flow across each face of all the tracer cells81 and sets the flux across the face to zero when the flux is from a dry tracer cell. This prevents cells82 whose depth is rn\_wdmin1 or less from drying out further. The scheme does not force $h$ (the water depth) at tracer cells83 to be at least the minimum depth and hence is able to conserve mass / volume.84 85 The flux across each $u$-face of a tracer cell is multiplied by a factor zuwdmask (an array which depends on ji and jj).86 If the user sets ln\_wd\_dl\_ramp = .False. then zuwdmask is 1 when the87 flux is from a cell with water depth greater than rn\_wdmin1 and 0 otherwise. If the user sets88 ln\_wd\_dl\_ramp = .True. the flux across the face is ramped down as the water depth decreases89 from 2 * rn\_wdmin1 to rn\_wdmin1. The use of this ramp reduced grid-scale noise in idealised test cases.90 91 At the point where the flux across a $u$-face is multiplied by zuwdmask , we have chosen92 also to multiply the corresponding velocity on the ``now'' step at that face by zuwdmask. We could have93 chosen not to do that and to allow fairly large velocities to occur in these ``dry'' cells.94 The rationale for setting the velocity to zero is that it is the momentum equations that are being solved95 and the total momentum of the upstream cell (treating it as a finite volume) should be considered96 to be its depth times its velocity. This depth is considered to be zero at ``dry'' $u$-points consistent with its97 treatment in the calculation of the flux of mass across the cell face.98 99 \cite{WarnerEtal13} state that in their scheme the velocity masks at the cell faces for the baroclinic100 timesteps are set to 0 or 1 depending on whether the average of the masks over the barotropic sub-steps is respectively less than101 or greater than 0.5. That scheme does not conserve tracers in integrations started from constant tracer102 fields (tracers independent of $x$, $y$ and $z$). Our scheme conserves constant tracers because103 the velocities used at the tracer cell faces on the baroclinic timesteps are carefully calculated by dynspg\_ts104 to equal their mean value during the barotropic steps. If the user sets ln\_wd\_dl\_bc = .True., the105 baroclinic velocities are also multiplied by a suitably weighted average of zuwdmask.106 107 %-----------------------------------------------------------------------------------------108 % Iterative limiters109 %-----------------------------------------------------------------------------------------110 \subsection [Iterative limiter (\textit{wet\_dry})]111 {Iterative limiter (\mdl{wet\_dry})}112 \label{DYN_wd_iterative_limiter}113 114 \subsubsection [Iterative flux limiter (\textit{wet\_dry})]115 {Iterative flux limiter (\mdl{wet\_dry})}116 \label{DYN_wd_il_spg_limiter}117 118 The iterative limiter modifies the fluxes across the faces of cells that are either already ``dry''119 or may become dry within the next time-step using an iterative method.120 121 The flux limiter for the barotropic flow (devised by Hedong Liu) can be understood as follows:122 123 The continuity equation for the total water depth in a column124 \begin{equation} \label{dyn_wd_continuity}125 \frac{\partial h}{\partial t} + \mathbf{\nabla.}(h\mathbf{u}) = 0 .126 \end{equation}127 can be written in discrete form as128 129 \begin{align} \label{dyn_wd_continuity_2}130 \frac{e_1 e_2}{\Delta t} ( h_{i,j}(t_{n+1}) - h_{i,j}(t_e) )131 &= - ( \mathrm{flxu}_{i+1,j} - \mathrm{flxu}_{i,j} + \mathrm{flxv}_{i,j+1} - \mathrm{flxv}_{i,j} ) \\132 &= \mathrm{zzflx}_{i,j} .133 \end{align}134 135 In the above $h$ is the depth of the water in the column at point $(i,j)$,136 $\mathrm{flxu}_{i+1,j}$ is the flux out of the ``eastern'' face of the cell and137 $\mathrm{flxv}_{i,j+1}$ the flux out of the ``northern'' face of the cell; $t_{n+1}$ is138 the new timestep, $t_e$ is the old timestep (either $t_b$ or $t_n$) and $ \Delta t =139 t_{n+1} - t_e$; $e_1 e_2$ is the area of the tracer cells centred at $(i,j)$ and140 $\mathrm{zzflx}$ is the sum of the fluxes through all the faces.141 142 The flux limiter splits the flux $\mathrm{zzflx}$ into fluxes that are out of the cell143 (zzflxp) and fluxes that are into the cell (zzflxn). Clearly144 145 \begin{equation} \label{dyn_wd_zzflx_p_n_1}146 \mathrm{zzflx}_{i,j} = \mathrm{zzflxp}_{i,j} + \mathrm{zzflxn}_{i,j} .147 \end{equation}148 149 The flux limiter iteratively adjusts the fluxes $\mathrm{flxu}$ and $\mathrm{flxv}$ until150 none of the cells will ``dry out''. To be precise the fluxes are limited until none of the151 cells has water depth less than $\mathrm{rn\_wdmin1}$ on step $n+1$.152 153 Let the fluxes on the $m$th iteration step be denoted by $\mathrm{flxu}^{(m)}$ and154 $\mathrm{flxv}^{(m)}$. Then the adjustment is achieved by seeking a set of coefficients,155 $\mathrm{zcoef}_{i,j}^{(m)}$ such that:156 157 \begin{equation} \label{dyn_wd_continuity_coef}158 \begin{split}159 \mathrm{zzflxp}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxp}^{(0)}_{i,j} \\160 \mathrm{zzflxn}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxn}^{(0)}_{i,j}161 \end{split}162 \end{equation}163 164 where the coefficients are $1.0$ generally but can vary between $0.0$ and $1.0$ around165 cells that would otherwise dry.166 167 The iteration is initialised by setting168 169 \begin{equation} \label{dyn_wd_zzflx_initial}170 \mathrm{zzflxp^{(0)}}_{i,j} = \mathrm{zzflxp}_{i,j} , \quad \mathrm{zzflxn^{(0)}}_{i,j} = \mathrm{zzflxn}_{i,j} .171 \end{equation}172 173 The fluxes out of cell $(i,j)$ are updated at the $m+1$th iteration if the depth of the174 cell on timestep $t_e$, namely $h_{i,j}(t_e)$, is less than the total flux out of the cell175 times the timestep divided by the cell area. Using (\ref{dyn_wd_continuity_2}) this176 condition is177 178 \begin{equation} \label{dyn_wd_continuity_if}179 h_{i,j}(t_e) - \mathrm{rn\_wdmin1} < \frac{\Delta t}{e_1 e_2} ( \mathrm{zzflxp}^{(m)}_{i,j} + \mathrm{zzflxn}^{(m)}_{i,j} ) .180 \end{equation}181 182 Rearranging (\ref{dyn_wd_continuity_if}) we can obtain an expression for the maximum183 outward flux that can be allowed and still maintain the minimum wet depth:184 185 \begin{equation} \label{dyn_wd_max_flux}186 \begin{split}187 \mathrm{zzflxp}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2}) \frac{e_1 e_2}{\Delta t} \phantom{]} \\188 \phantom{[} & - \mathrm{zzflxn}^{(m)}_{i,j} \Big]189 \end{split}190 \end{equation}191 192 Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\it [Q: Why is193 this necessary/desirable?]}. Substituting from (\ref{dyn_wd_continuity_coef}) gives an194 expression for the coefficient needed to multiply the outward flux at this cell in order195 to avoid drying.196 197 \begin{equation} \label{dyn_wd_continuity_nxtcoef}198 \begin{split}199 \mathrm{zcoef}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2}) \frac{e_1 e_2}{\Delta t} \phantom{]} \\200 \phantom{[} & - \mathrm{zzflxn}^{(m)}_{i,j} \Big] \frac{1}{ \mathrm{zzflxp}^{(0)}_{i,j} }201 \end{split}202 \end{equation}203 204 Only the outward flux components are altered but, of course, outward fluxes from one cell205 are inward fluxes to adjacent cells and the balance in these cells may need subsequent206 adjustment; hence the iterative nature of this scheme. Note, for example, that the flux207 across the ``eastern'' face of the $(i,j)$th cell is only updated at the $m+1$th iteration208 if that flux at the $m$th iteration is out of the $(i,j)$th cell. If that is the case then209 the flux across that face is into the $(i+1,j)$ cell and that flux will not be updated by210 the calculation for the $(i+1,j)$th cell. In this sense the updates to the fluxes across211 the faces of the cells do not ``compete'' (they do not over-write each other) and one212 would expect the scheme to converge relatively quickly. The scheme is flux based so213 conserves mass. It also conserves constant tracers for the same reason that the214 directional limiter does.215 216 217 %----------------------------------------------------------------------------------------218 % Surface pressure gradients219 %----------------------------------------------------------------------------------------220 \subsubsection [Modification of surface pressure gradients (\textit{dynhpg})]221 {Modification of surface pressure gradients (\mdl{dynhpg})}222 \label{DYN_wd_il_spg}223 224 At ``dry'' points the water depth is usually close to $\mathrm{rn\_wdmin1}$. If the225 topography is sloping at these points the sea-surface will have a similar slope and there226 will hence be very large horizontal pressure gradients at these points. The WAD modifies227 the magnitude but not the sign of the surface pressure gradients (zhpi and zhpj) at such228 points by mulitplying them by positive factors (zcpx and zcpy respectively) that lie229 between $0$ and $1$.230 231 We describe how the scheme works for the ``eastward'' pressure gradient, zhpi, calculated232 at the $(i,j)$th $u$-point. The scheme uses the ht\_wd depths and surface heights at the233 neighbouring $(i+1,j)$ and $(i,j)$ tracer points. zcpx is calculated using two logicals234 variables, $\mathrm{ll\_tmp1}$ and $\mathrm{ll\_tmp2}$ which are evaluated for each grid235 column. The three possible combinations are illustrated in figure \ref{Fig_WAD_dynhpg}.236 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>237 \begin{figure}[!ht] \begin{center}238 \includegraphics[width=0.8\textwidth]{Fig_WAD_dynhpg}239 \caption{ \label{Fig_WAD_dynhpg}240 Illustrations of the three possible combinations of the logical variables controlling the241 limiting of the horizontal pressure gradient in wetting and drying regimes}242 \end{center}\end{figure}243 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>244 245 The first logical, $\mathrm{ll\_tmp1}$, is set to true if and only if the water depth at246 both neighbouring points is greater than $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$ and247 the minimum height of the sea surface at the two points is greater than the maximum height248 of the topography at the two points:249 250 \begin{equation} \label{dyn_ll_tmp1}251 \begin{split}252 \mathrm{ll\_tmp1} = & \mathrm{MIN(sshn(ji,jj), sshn(ji+1,jj))} > \\253 & \quad \mathrm{MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj))\ .and.} \\254 & \mathrm{MAX(sshn(ji,jj) + ht\_wd(ji,jj),} \\255 & \mathrm{\phantom{MAX(}sshn(ji+1,jj) + ht\_wd(ji+1,jj))} >\\256 & \quad\quad\mathrm{rn\_wdmin1 + rn\_wdmin2 }257 \end{split}258 \end{equation}259 260 The second logical, $\mathrm{ll\_tmp2}$, is set to true if and only if the maximum height261 of the sea surface at the two points is greater than the maximum height of the topography262 at the two points plus $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$263 264 \begin{equation} \label{dyn_ll_tmp2}265 \begin{split}266 \mathrm{ ll\_tmp2 } = & \mathrm{( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 )\ .AND.}\\267 & \mathrm{( MAX(sshn(ji,jj), sshn(ji+1,jj)) > } \\268 & \mathrm{\phantom{(} MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj)) + rn\_wdmin1 + rn\_wdmin2}) .269 \end{split}270 \end{equation}271 272 If $\mathrm{ll\_tmp1}$ is true then the surface pressure gradient, zhpi at the $(i,j)$273 point is unmodified. If both logicals are false zhpi is set to zero.274 275 If $\mathrm{ll\_tmp1}$ is true and $\mathrm{ll\_tmp2}$ is false then the surface pressure276 gradient is multiplied through by zcpx which is the absolute value of the difference in277 the water depths at the two points divided by the difference in the surface heights at the278 two points. Thus the sign of the sea surface height gradient is retained but the magnitude279 of the pressure force is determined by the difference in water depths rather than the280 difference in surface height between the two points. Note that dividing by the difference281 between the sea surface heights can be problematic if the heights approach parity. An282 additional condition is applied to $\mathrm{ ll\_tmp2 }$ to ensure it is .false. in such283 conditions.284 285 \subsection [Additional considerations (\textit{usrdef\_zgr})]286 {Additional considerations (\mdl{usrdef\_zgr})}287 \label{WAD_additional}288 289 In the very shallow water where wetting and drying occurs the parametrisation of290 bottom drag is clearly very important. In order to promote stability291 it is sometimes useful to calculate the bottom drag using an implicit time-stepping approach.292 293 Suitable specifcation of the surface heat flux in wetting and drying domains in forced and294 coupled simulations needs further consideration. In order to prevent freezing or boiling295 in uncoupled integrations the net surface heat fluxes need to be appropriately limited.296 27 297 28 %---------------------------------------------------------------------------------------- 298 29 % The WAD test cases 299 30 %---------------------------------------------------------------------------------------- 300 \s ubsection [The WAD test cases (\textit{usrdef\_zgr})]31 \section [The WAD test cases (\textit{usrdef\_zgr})] 301 32 {The WAD test cases (\mdl{usrdef\_zgr})} 302 33 \label{WAD_test_cases} … … 327 58 328 59 \clearpage 329 \subs ubsection [WAD test case 1 : A simple linear slope]60 \subsection [WAD test case 1 : A simple linear slope] 330 61 {WAD test case 1 : A simple linear slope} 331 62 \label{WAD_test_case1} … … 351 82 352 83 \clearpage 353 \subs ubsection [WAD test case 2 : A parabolic channel ]84 \subsection [WAD test case 2 : A parabolic channel ] 354 85 {WAD test case 2 : A parabolic channel} 355 86 \label{WAD_test_case2} … … 374 105 375 106 \clearpage 376 \subs ubsection [WAD test case 3 : A parabolic channel (extreme slope) ]107 \subsection [WAD test case 3 : A parabolic channel (extreme slope) ] 377 108 {WAD test case 3 : A parabolic channel (extreme slope)} 378 109 \label{WAD_test_case3} … … 393 124 394 125 \clearpage 395 \subs ubsection [WAD test case 4 : A parabolic bowl ]126 \subsection [WAD test case 4 : A parabolic bowl ] 396 127 {WAD test case 4 : A parabolic bowl} 397 128 \label{WAD_test_case4} … … 415 146 416 147 \clearpage 417 \subs ubsection [WAD test case 5 : A double slope with shelf channel ]148 \subsection [WAD test case 5 : A double slope with shelf channel ] 418 149 {WAD test case 5 : A double slope with shelf channel} 419 150 \label{WAD_test_case5} … … 434 165 435 166 \clearpage 436 \subs ubsection [WAD test case 6 : A parabolic channel with central bar ]167 \subsection [WAD test case 6 : A parabolic channel with central bar ] 437 168 {WAD test case 6 : A parabolic channel with central bar} 438 169 \label{WAD_test_case6} … … 463 194 464 195 \clearpage 465 \subs ubsection [WAD test case 7 : A double slope with shelf, open-ended channel ]196 \subsection [WAD test case 7 : A double slope with shelf, open-ended channel ] 466 197 {WAD test case 7 : A double slope with shelf, open-ended channel} 467 198 \label{WAD_test_case7} … … 511 242 % ================================================================ 512 243 513 \bibliographystyle{wileyqj} 514 \bibliography{references} 244 %\bibliographystyle{wileyqj} 245 %\bibliographystyle{../../../doc/latex/NEMO/main/ametsoc.bst} 246 %\bibliography{references} 515 247 516 248 \end{document} -
NEMO/trunk/tests/WAD/MY_DOCS/references.bib
r9024 r10499 8 8 9 9 @article{Oey06, 10 author = {L-Y Oey}, 11 title = {An OGCM with movable land--sea boundaries}, 12 journal = {Ocean Mod.}, 13 year = {2006}, 14 volume = {13}, 15 pages = {176--195}} 10 title = "An OGCM with movable land-sea boundaries", 11 journal = "Ocean Modelling", 12 volume = "13", 13 number = "2", 14 pages = "176 - 195", 15 year = "2006", 16 issn = "1463-5003", 17 doi = "https://doi.org/10.1016/j.ocemod.2006.01.001", 18 url = "http://www.sciencedirect.com/science/article/pii/S1463500306000084", 19 author = "Lie-Yauw Oey", 20 keywords = "Wetting and drying, Inundations, Ocean general circulation model (OGCM), Princeton Ocean Model (POM), Tides, Tsunamis, Estuarine outflows", 21 abstract = "An ocean general circulation model (OGCM) with wetting and drying (WAD) capabilities removes the vertical-wall coastal assumption and allows simultaneous modeling of open-ocean currents and water run-up (and run-down) across movable land-sea boundaries. This paper implements and tests such a WAD scheme for the Princeton Ocean Model (POM) in its most general three-dimensional setting with stratification, bathymetry and forcing. The scheme can be easily exported to other OGCM's." 22 } 16 23 17 24 @article{WarnerEtal13, 18 author = {J C Warner and Z Defne and K Haas and H G Arango}, 19 title = {{A wetting and drying scheme for ROMS}}, 20 journal = {Computers and Geosc.}, 21 year = {2013}, 22 volume = {58}, 23 pages = {54--61}} 25 title = "A wetting and drying scheme for ROMS", 26 journal = "Computers \& Geosciences", 27 volume = "58", 28 pages = "54 - 61", 29 year = "2013", 30 issn = "0098-3004", 31 doi = "https://doi.org/10.1016/j.cageo.2013.05.004", 32 url = "http://www.sciencedirect.com/science/article/pii/S0098300413001362", 33 author = "John C. Warner and Zafer Defne and Kevin Haas and Hernan G. Arango", 34 keywords = "Wetting and drying, ROMS, Cell-face blocking", 35 abstract = "The processes of wetting and drying have many important physical and biological impacts on shallow water systems. Inundation and dewatering effects on coastal mud flats and beaches occur on various time scales ranging from storm surge, periodic rise and fall of the tide, to infragravity wave motions. To correctly simulate these physical processes with a numerical model requires the capability of the computational cells to become inundated and dewatered. In this paper, we describe a method for wetting and drying based on an approach consistent with a cell-face blocking algorithm. The method allows water to always flow into any cell, but prevents outflow from a cell when the total depth in that cell is less than a user defined critical value. We describe the method, the implementation into the three-dimensional Regional Oceanographic Modeling System (ROMS), and exhibit the new capability under three scenarios: an analytical expression for shallow water flows, a dam break test case, and a realistic application to part of a wetland area along the Georgia Coast, USA." 36 }
Note: See TracChangeset
for help on using the changeset viewer.