 Timestamp:
 20190108T17:17:38+01:00 (5 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

NEMO/branches/UKMO/dev_10448_WAD_SBC_BUGFIX/tests/WAD/MY_DOCS/WAD_doc.tex
r10467 r10472 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 sealevel. 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 nonzero 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 timestep. Surface fluxes58 are also switched off for dry cells to prevent freezing, boiling etc. of very thin water layers.59 The fluxes are tappered down using a $\mathrm{tanh}$ weighting function60 to no flux as the dry limit $\mathrm{rn\_wdmin1}$ is approached. Even wet cells can be very shallow.61 The depth at which to start tapering is controlled by the user by setting $\mathrm{rn\_wd\_sbcdep}$.62 The fraction $(<1)$ of sufrace fluxes to use at this depth is set by $\mathrm{rn\_wd\_sbcfra}$.63 64 Both versions of the code have been tested in six test cases provided in the WAD\_TEST\_CASES configuration65 and in ``realistic'' configurations covering parts of the northwest European shelf.66 All these configurations have used pure sigma coordinates. It is expected that67 the wetting and drying code will work in domains with more general scoordinates provided68 the coordinates are pure sigma in the region where wetting and drying actually occurs.69 70 The next subsection descrbies the directional limiter and the following subsection the iterative limiter.71 The final subsection covers some additional considerations that are relevant to both schemes.72 73 %74 % Iterative limiters75 %76 \subsection [Directional limiter (\textit{wet\_dry})]77 {Directional limiter (\mdl{wet\_dry})}78 \label{DYN_wd_directional_limiter}79 80 The principal idea of the directional limiter is that81 water should not be allowed to flow out of a dry tracer cell (i.e. one whose water depth is less than rn\_wdmin1).82 83 All the changes associated with this option are made to the barotropic solver for the nonlinear84 free surface code within dynspg\_ts.85 On each barotropic substep the scheme determines the direction of the flow across each face of all the tracer cells86 and sets the flux across the face to zero when the flux is from a dry tracer cell. This prevents cells87 whose depth is rn\_wdmin1 or less from drying out further. The scheme does not force $h$ (the water depth) at tracer cells88 to be at least the minimum depth and hence is able to conserve mass / volume.89 90 The flux across each $u$face of a tracer cell is multiplied by a factor zuwdmask (an array which depends on ji and jj).91 If the user sets ln\_wd\_dl\_ramp = .False. then zuwdmask is 1 when the92 flux is from a cell with water depth greater than rn\_wdmin1 and 0 otherwise. If the user sets93 ln\_wd\_dl\_ramp = .True. the flux across the face is ramped down as the water depth decreases94 from 2 * rn\_wdmin1 to rn\_wdmin1. The use of this ramp reduced gridscale noise in idealised test cases.95 96 At the point where the flux across a $u$face is multiplied by zuwdmask , we have chosen97 also to multiply the corresponding velocity on the ``now'' step at that face by zuwdmask. We could have98 chosen not to do that and to allow fairly large velocities to occur in these ``dry'' cells.99 The rationale for setting the velocity to zero is that it is the momentum equations that are being solved100 and the total momentum of the upstream cell (treating it as a finite volume) should be considered101 to be its depth times its velocity. This depth is considered to be zero at ``dry'' $u$points consistent with its102 treatment in the calculation of the flux of mass across the cell face.103 104 \cite{WarnerEtal13} state that in their scheme the velocity masks at the cell faces for the baroclinic105 timesteps are set to 0 or 1 depending on whether the average of the masks over the barotropic substeps is respectively less than106 or greater than 0.5. That scheme does not conserve tracers in integrations started from constant tracer107 fields (tracers independent of $x$, $y$ and $z$). Our scheme conserves constant tracers because108 the velocities used at the tracer cell faces on the baroclinic timesteps are carefully calculated by dynspg\_ts109 to equal their mean value during the barotropic steps. If the user sets ln\_wd\_dl\_bc = .True., the110 baroclinic velocities are also multiplied by a suitably weighted average of zuwdmask.111 112 %113 % Iterative limiters114 %115 \subsection [Iterative limiter (\textit{wet\_dry})]116 {Iterative limiter (\mdl{wet\_dry})}117 \label{DYN_wd_iterative_limiter}118 119 \subsubsection [Iterative flux limiter (\textit{wet\_dry})]120 {Iterative flux limiter (\mdl{wet\_dry})}121 \label{DYN_wd_il_spg_limiter}122 123 The iterative limiter modifies the fluxes across the faces of cells that are either already ``dry''124 or may become dry within the next timestep using an iterative method.125 126 The flux limiter for the barotropic flow (devised by Hedong Liu) can be understood as follows:127 128 The continuity equation for the total water depth in a column129 \begin{equation} \label{dyn_wd_continuity}130 \frac{\partial h}{\partial t} + \mathbf{\nabla.}(h\mathbf{u}) = 0 .131 \end{equation}132 can be written in discrete form as133 134 \begin{align} \label{dyn_wd_continuity_2}135 \frac{e_1 e_2}{\Delta t} ( h_{i,j}(t_{n+1})  h_{i,j}(t_e) )136 &=  ( \mathrm{flxu}_{i+1,j}  \mathrm{flxu}_{i,j} + \mathrm{flxv}_{i,j+1}  \mathrm{flxv}_{i,j} ) \\137 &= \mathrm{zzflx}_{i,j} .138 \end{align}139 140 In the above $h$ is the depth of the water in the column at point $(i,j)$,141 $\mathrm{flxu}_{i+1,j}$ is the flux out of the ``eastern'' face of the cell and142 $\mathrm{flxv}_{i,j+1}$ the flux out of the ``northern'' face of the cell; $t_{n+1}$ is143 the new timestep, $t_e$ is the old timestep (either $t_b$ or $t_n$) and $ \Delta t =144 t_{n+1}  t_e$; $e_1 e_2$ is the area of the tracer cells centred at $(i,j)$ and145 $\mathrm{zzflx}$ is the sum of the fluxes through all the faces.146 147 The flux limiter splits the flux $\mathrm{zzflx}$ into fluxes that are out of the cell148 (zzflxp) and fluxes that are into the cell (zzflxn). Clearly149 150 \begin{equation} \label{dyn_wd_zzflx_p_n_1}151 \mathrm{zzflx}_{i,j} = \mathrm{zzflxp}_{i,j} + \mathrm{zzflxn}_{i,j} .152 \end{equation}153 154 The flux limiter iteratively adjusts the fluxes $\mathrm{flxu}$ and $\mathrm{flxv}$ until155 none of the cells will ``dry out''. To be precise the fluxes are limited until none of the156 cells has water depth less than $\mathrm{rn\_wdmin1}$ on step $n+1$.157 158 Let the fluxes on the $m$th iteration step be denoted by $\mathrm{flxu}^{(m)}$ and159 $\mathrm{flxv}^{(m)}$. Then the adjustment is achieved by seeking a set of coefficients,160 $\mathrm{zcoef}_{i,j}^{(m)}$ such that:161 162 \begin{equation} \label{dyn_wd_continuity_coef}163 \begin{split}164 \mathrm{zzflxp}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxp}^{(0)}_{i,j} \\165 \mathrm{zzflxn}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxn}^{(0)}_{i,j}166 \end{split}167 \end{equation}168 169 where the coefficients are $1.0$ generally but can vary between $0.0$ and $1.0$ around170 cells that would otherwise dry.171 172 The iteration is initialised by setting173 174 \begin{equation} \label{dyn_wd_zzflx_initial}175 \mathrm{zzflxp^{(0)}}_{i,j} = \mathrm{zzflxp}_{i,j} , \quad \mathrm{zzflxn^{(0)}}_{i,j} = \mathrm{zzflxn}_{i,j} .176 \end{equation}177 178 The fluxes out of cell $(i,j)$ are updated at the $m+1$th iteration if the depth of the179 cell on timestep $t_e$, namely $h_{i,j}(t_e)$, is less than the total flux out of the cell180 times the timestep divided by the cell area. Using (\ref{dyn_wd_continuity_2}) this181 condition is182 183 \begin{equation} \label{dyn_wd_continuity_if}184 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} ) .185 \end{equation}186 187 Rearranging (\ref{dyn_wd_continuity_if}) we can obtain an expression for the maximum188 outward flux that can be allowed and still maintain the minimum wet depth:189 190 \begin{equation} \label{dyn_wd_max_flux}191 \begin{split}192 \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{]} \\193 \phantom{[} &  \mathrm{zzflxn}^{(m)}_{i,j} \Big]194 \end{split}195 \end{equation}196 197 Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\it [Q: Why is198 this necessary/desirable?]}. Substituting from (\ref{dyn_wd_continuity_coef}) gives an199 expression for the coefficient needed to multiply the outward flux at this cell in order200 to avoid drying.201 202 \begin{equation} \label{dyn_wd_continuity_nxtcoef}203 \begin{split}204 \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{]} \\205 \phantom{[} &  \mathrm{zzflxn}^{(m)}_{i,j} \Big] \frac{1}{ \mathrm{zzflxp}^{(0)}_{i,j} }206 \end{split}207 \end{equation}208 209 Only the outward flux components are altered but, of course, outward fluxes from one cell210 are inward fluxes to adjacent cells and the balance in these cells may need subsequent211 adjustment; hence the iterative nature of this scheme. Note, for example, that the flux212 across the ``eastern'' face of the $(i,j)$th cell is only updated at the $m+1$th iteration213 if that flux at the $m$th iteration is out of the $(i,j)$th cell. If that is the case then214 the flux across that face is into the $(i+1,j)$ cell and that flux will not be updated by215 the calculation for the $(i+1,j)$th cell. In this sense the updates to the fluxes across216 the faces of the cells do not ``compete'' (they do not overwrite each other) and one217 would expect the scheme to converge relatively quickly. The scheme is flux based so218 conserves mass. It also conserves constant tracers for the same reason that the219 directional limiter does.220 221 222 %223 % Surface pressure gradients224 %225 \subsubsection [Modification of surface pressure gradients (\textit{dynhpg})]226 {Modification of surface pressure gradients (\mdl{dynhpg})}227 \label{DYN_wd_il_spg}228 229 At ``dry'' points the water depth is usually close to $\mathrm{rn\_wdmin1}$. If the230 topography is sloping at these points the seasurface will have a similar slope and there231 will hence be very large horizontal pressure gradients at these points. The WAD modifies232 the magnitude but not the sign of the surface pressure gradients (zhpi and zhpj) at such233 points by mulitplying them by positive factors (zcpx and zcpy respectively) that lie234 between $0$ and $1$.235 236 We describe how the scheme works for the ``eastward'' pressure gradient, zhpi, calculated237 at the $(i,j)$th $u$point. The scheme uses the ht\_wd depths and surface heights at the238 neighbouring $(i+1,j)$ and $(i,j)$ tracer points. zcpx is calculated using two logicals239 variables, $\mathrm{ll\_tmp1}$ and $\mathrm{ll\_tmp2}$ which are evaluated for each grid240 column. The three possible combinations are illustrated in figure \ref{Fig_WAD_dynhpg}.241 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>242 \begin{figure}[!ht] \begin{center}243 \includegraphics[width=0.8\textwidth]{Fig_WAD_dynhpg}244 \caption{ \label{Fig_WAD_dynhpg}245 Illustrations of the three possible combinations of the logical variables controlling the246 limiting of the horizontal pressure gradient in wetting and drying regimes}247 \end{center}\end{figure}248 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>249 250 The first logical, $\mathrm{ll\_tmp1}$, is set to true if and only if the water depth at251 both neighbouring points is greater than $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$ and252 the minimum height of the sea surface at the two points is greater than the maximum height253 of the topography at the two points:254 255 \begin{equation} \label{dyn_ll_tmp1}256 \begin{split}257 \mathrm{ll\_tmp1} = & \mathrm{MIN(sshn(ji,jj), sshn(ji+1,jj))} > \\258 & \quad \mathrm{MAX(ht\_wd(ji,jj), ht\_wd(ji+1,jj))\ .and.} \\259 & \mathrm{MAX(sshn(ji,jj) + ht\_wd(ji,jj),} \\260 & \mathrm{\phantom{MAX(}sshn(ji+1,jj) + ht\_wd(ji+1,jj))} >\\261 & \quad\quad\mathrm{rn\_wdmin1 + rn\_wdmin2 }262 \end{split}263 \end{equation}264 265 The second logical, $\mathrm{ll\_tmp2}$, is set to true if and only if the maximum height266 of the sea surface at the two points is greater than the maximum height of the topography267 at the two points plus $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$268 269 \begin{equation} \label{dyn_ll_tmp2}270 \begin{split}271 \mathrm{ ll\_tmp2 } = & \mathrm{( ABS( sshn(ji,jj)  sshn(ji+1,jj) ) > 1.E12 )\ .AND.}\\272 & \mathrm{( MAX(sshn(ji,jj), sshn(ji+1,jj)) > } \\273 & \mathrm{\phantom{(} MAX(ht\_wd(ji,jj), ht\_wd(ji+1,jj)) + rn\_wdmin1 + rn\_wdmin2}) .274 \end{split}275 \end{equation}276 277 If $\mathrm{ll\_tmp1}$ is true then the surface pressure gradient, zhpi at the $(i,j)$278 point is unmodified. If both logicals are false zhpi is set to zero.279 280 If $\mathrm{ll\_tmp1}$ is true and $\mathrm{ll\_tmp2}$ is false then the surface pressure281 gradient is multiplied through by zcpx which is the absolute value of the difference in282 the water depths at the two points divided by the difference in the surface heights at the283 two points. Thus the sign of the sea surface height gradient is retained but the magnitude284 of the pressure force is determined by the difference in water depths rather than the285 difference in surface height between the two points. Note that dividing by the difference286 between the sea surface heights can be problematic if the heights approach parity. An287 additional condition is applied to $\mathrm{ ll\_tmp2 }$ to ensure it is .false. in such288 conditions.289 290 \subsection [Additional considerations (\textit{usrdef\_zgr})]291 {Additional considerations (\mdl{usrdef\_zgr})}292 \label{WAD_additional}293 294 In the very shallow water where wetting and drying occurs the parametrisation of295 bottom drag is clearly very important. In order to promote stability296 it is sometimes useful to calculate the bottom drag using an implicit timestepping approach.297 298 Suitable specifcation of the surface heat flux in wetting and drying domains in forced and299 coupled simulations needs further consideration. In order to prevent freezing or boiling300 in uncoupled integrations the net surface heat fluxes need to be appropriately limited.301 27 302 28 % 303 29 % The WAD test cases 304 30 % 305 \s ubsection [The WAD test cases (\textit{usrdef\_zgr})]31 \section [The WAD test cases (\textit{usrdef\_zgr})] 306 32 {The WAD test cases (\mdl{usrdef\_zgr})} 307 33 \label{WAD_test_cases} … … 332 58 333 59 \clearpage 334 \subs ubsection [WAD test case 1 : A simple linear slope]60 \subsection [WAD test case 1 : A simple linear slope] 335 61 {WAD test case 1 : A simple linear slope} 336 62 \label{WAD_test_case1} … … 356 82 357 83 \clearpage 358 \subs ubsection [WAD test case 2 : A parabolic channel ]84 \subsection [WAD test case 2 : A parabolic channel ] 359 85 {WAD test case 2 : A parabolic channel} 360 86 \label{WAD_test_case2} … … 379 105 380 106 \clearpage 381 \subs ubsection [WAD test case 3 : A parabolic channel (extreme slope) ]107 \subsection [WAD test case 3 : A parabolic channel (extreme slope) ] 382 108 {WAD test case 3 : A parabolic channel (extreme slope)} 383 109 \label{WAD_test_case3} … … 398 124 399 125 \clearpage 400 \subs ubsection [WAD test case 4 : A parabolic bowl ]126 \subsection [WAD test case 4 : A parabolic bowl ] 401 127 {WAD test case 4 : A parabolic bowl} 402 128 \label{WAD_test_case4} … … 420 146 421 147 \clearpage 422 \subs ubsection [WAD test case 5 : A double slope with shelf channel ]148 \subsection [WAD test case 5 : A double slope with shelf channel ] 423 149 {WAD test case 5 : A double slope with shelf channel} 424 150 \label{WAD_test_case5} … … 439 165 440 166 \clearpage 441 \subs ubsection [WAD test case 6 : A parabolic channel with central bar ]167 \subsection [WAD test case 6 : A parabolic channel with central bar ] 442 168 {WAD test case 6 : A parabolic channel with central bar} 443 169 \label{WAD_test_case6} … … 468 194 469 195 \clearpage 470 \subs ubsection [WAD test case 7 : A double slope with shelf, openended channel ]196 \subsection [WAD test case 7 : A double slope with shelf, openended channel ] 471 197 {WAD test case 7 : A double slope with shelf, openended channel} 472 198 \label{WAD_test_case7} … … 517 243 518 244 %\bibliographystyle{wileyqj} 519 \bibliographystyle{../../../doc/latex/NEMO/main/ametsoc.bst}520 \bibliography{references}245 %\bibliographystyle{../../../doc/latex/NEMO/main/ametsoc.bst} 246 %\bibliography{references} 521 247 522 248 \end{document}
Note: See TracChangeset
for help on using the changeset viewer.