/[MITgcm]/manual/s_outp_pkgs/text/pvdiag.tex
ViewVC logotype

Contents of /manual/s_outp_pkgs/text/pvdiag.tex

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (show annotations) (download) (as text)
Tue Jan 15 21:58:33 2008 UTC (16 years, 4 months ago) by cnh
Branch: MAIN
File MIME type: application/x-tex
Adding Guillaume's pvdiag section

1 \section{Potential vorticity Matlab Toolbox}
2 \label{sec:diag:pv}
3 \begin{rawhtml}
4 %<!-- CMIREDIR:pvdiag_matlab: -->
5 \end{rawhtml}
6
7 %%%%%%%%%%%%%%%%%%%%%%%%%%
8 \subsection{Introduction}
9
10 \noindent
11 This section of the documentation describes a Matlab package that aims to provide useful
12 routines to compute vorticity fields (relative, potential and planetary) and its related
13 components. This is an offline computation. It was developed to be used in mode water
14 studies, so that it comes with other related routines, in particular ones computing
15 surface vertical potential vorticity fluxes.
16
17 %%%%%%%%%%%%%%%%%%%%%%%%%%
18 \subsection{Equations}
19
20 \subsubsection{Potential Vorticity}
21 \noindent
22 The package computes the three components of the relative vorticity defined by:
23 \begin{eqnarray}
24 \omega &= \nabla \times {\bf U} = \left( \begin{array}{c}
25 \omega_x\\
26 \omega_y\\
27 \zeta
28 \end{array}\right)
29 \simeq &\left( \begin{array}{c}
30 -\frac{\partial v}{\partial z}\\
31 -\frac{\partial u}{\partial z}\\
32 \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}
33 \end{array}\right) \label{sec:diag:pv:eq1}
34 \end{eqnarray}
35 where we omitted (like all across the package) the vertical velocity component.
36
37 The package then computes the potential vorticity as:
38 \begin{eqnarray}
39 Q &=& -\frac{1}{\rho} \omega\cdot\nabla\sigma_\theta\\
40 Q &=& -\frac{1}{\rho}\left(\omega_x \frac{\partial \sigma_\theta}{\partial x} +
41 \omega_y \frac{\partial \sigma_\theta}{\partial y} +
42 \left(f+\zeta\right) \frac{\partial \sigma_\theta}{\partial z}\right) \label{sec:diag:pv:eq2}
43 \end{eqnarray}
44 where $\rho$ is the density, $\sigma_\theta$ is the potential density (both eventually
45 computed by the package) and $f$ is the Coriolis parameter.
46
47 The package is also able to compute the simpler planetary vorticity as:
48 \begin{eqnarray}
49 splQ &=& -\frac{f}{\rho}\frac{\sigma_\theta}{\partial z} \label{sec:diag:pv:eq3}
50 \end{eqnarray}
51
52 \subsubsection{Surface vertical potential vorticity fluxes}
53 These quantities are useful in mode water studies because of the impermeability theorem
54 which states that for a given potential density layer (embedding a mode water), the
55 integrated PV only changes through surface input/output.
56
57 Vertical PV fluxes due to frictional and diabatic processes are given by:
58 \begin{eqnarray}
59 J^B_z &=& -\frac{f}{h}\left( \frac{\alpha Q_{net}}{C_w}-\rho_0 \beta S_{net}\right) \label{sec:diag:pv:eq14a}\\
60 J^F_z &=& \frac{1}{\rho\delta_e} \vec{k}\times\tau\cdot\nabla\sigma_m \label{sec:diag:pv:eq15a}
61 \end{eqnarray}
62 These components can be computed with the package. Details on the variables definition and
63 the way these fluxes are derived can be found in section
64 \ref{sec:diag:pv:notes-flux-form}.
65
66 \vspace{.5cm}
67 \noindent
68 We now give some simple explanations about these fluxes and how they can reduce the PV
69 level of an oceanic potential density layer.
70
71 \paragraph{Diabatic process}
72 Let's take the PV flux due to surface buoyancy forcing from Eq.\ref{sec:diag:pv:eq14a} and
73 simplify it as:
74 \begin{eqnarray}
75 J^B_z &\simeq& -\frac{\alpha f}{hC_w} Q_{net}
76 \end{eqnarray}
77 When the net surface heat flux $Q_{net}$ is upward i.e. negative and cooling the ocean
78 (buoyancy loss), surface density will increase, triggering mixing which reduces the
79 stratification and then the PV.
80 \begin{eqnarray}
81 Q_{net} &<& 0 \,\,\,\hbox{(upward, cooling)} \nonumber \\
82 J^B_z &>& 0 \,\,\,\hbox{(upward)} \nonumber \\
83 -\rho^{-1}\nabla\cdot J^B_z &<& 0 \,\,\, \hbox{(PV flux divergence)} \nonumber \\
84 PV &\searrow& \hbox{where $Q_{net}<0$}\nonumber
85 \end{eqnarray}
86
87
88 \paragraph{Frictional process: "Down-front" wind-stress}
89 Now let's take the PV flux due to the "wind-driven buoyancy flux" from
90 Eq.\ref{sec:diag:pv:eq15a} and simplify it as:
91 \begin{eqnarray}
92 J^F_z &=& \frac{1}{\rho\delta_e} \left( \tau_x\frac{\partial \sigma}{\partial y} - \tau_y\frac{\partial \sigma}{\partial x} \right) \\
93 J^F_z &\simeq& \frac{1}{\rho\delta_e} \tau_x\frac{\partial \sigma}{\partial y} \nonumber
94 \end{eqnarray}
95 When the wind is blowing from the east above the Gulf Stream (a region of high meridional
96 density gradient), it induces an advection of dense water from the northern side of the GS
97 to the southern side through Ekman currents. Then, it induces a "wind-driven" buoyancy
98 lost and mixing which reduces the stratification and the PV.
99 \begin{eqnarray}
100 \vec{k}\times\tau\cdot\nabla\sigma &>& 0 \,\,\, \hbox{("Down-front" wind)} \nonumber \\
101 J^F_z &>& 0 \,\,\, \hbox{(upward)} \nonumber \\
102 -\rho^{-1}\nabla\cdot J^F_z &<& 0 \,\,\, \hbox{(PV flux divergence)} \nonumber \\
103 PV &\searrow& \hbox{where $\vec{k}\times\tau\cdot\nabla\sigma>0$}\nonumber
104 \end{eqnarray}
105
106
107 \paragraph{Diabatic versus frictional processes}
108 A recent debate in the community arose about the relative role of these processes. Taking
109 the ratio of Eq.\ref{sec:diag:pv:eq14a} and Eq.\ref{sec:diag:pv:eq15a} leads to:
110 \begin{eqnarray}
111 \frac{J^F_z}{J^B_Z} &=& \frac{ \frac{1}{\rho\delta_e} \vec{k}\times\tau\cdot\nabla\sigma }
112 {-\frac{f}{h}\left( \frac{\alpha Q_{net}}{C_w}-\rho_0 \beta S_{net}\right)} \\
113 &\simeq& \frac{Q_{Ek}/\delta_e}{Q_{net}/h} \nonumber
114 \end{eqnarray}
115 where appears the lateral heat flux induced by Ekman currents:
116 \begin{eqnarray}
117 Q_{Ek} &=& -\frac{C_w}{\alpha\rho f}\vec{k}\times\tau\cdot\nabla\sigma
118 \nonumber \\
119 &=& \frac{C_w}{\alpha}\delta_e\vec{u_{Ek}}\cdot\nabla\sigma
120 \end{eqnarray}
121 which can be computed with the package. In the aim of comparing both processes, it will be
122 useful to plot surface net and lateral Ekman-induced heat fluxes together with PV fluxes.
123
124
125
126
127 %%%%%%%%%%%%%%%%%%%%%%%%%%
128 \subsection{Key routines}
129
130 \begin{itemize}
131 \item {\bf {\ttfamily A\_compute\_potential\_density.m}}: Compute the potential density
132 field. Requires the potential temperature and salinity (either total or anomalous) and
133 produces one output file with the potential density field (file prefix is {\ttfamily
134 SIGMATHETA}). The routine uses {\ttfamily densjmd95.m} a Matlab counterpart of the
135 MITgcm built-in function to compute the density.
136 \item {\bf {\ttfamily B\_compute\_relative\_vorticity.m}}: Compute the three components of
137 the relative vorticity defined in Eq.~(\ref{sec:diag:pv:eq1}). Requires the two
138 horizontal velocity components and produces three output files with the three components
139 (files prefix are {\ttfamily OMEGAX}, {\ttfamily OMEGAY} and {\ttfamily ZETA}).
140 \item {\bf {\ttfamily C\_compute\_potential\_vorticity.m}}: Compute the potential
141 vorticity without the negative ratio by the density. Two options are possible in order
142 to compute either the full component (term into parenthesis in
143 Eq.~\ref{sec:diag:pv:eq2}) or the planetary component ($f\partial_z\sigma_\theta$ in
144 Eq.~\ref{sec:diag:pv:eq3}). Requires the relative vorticity components and the potential
145 density, and produces one output file with the potential vorticity (file prefix is
146 {\ttfamily PV} for the full term and {\ttfamily splPV} for the planetary component).
147 \item {\bf {\ttfamily D\_compute\_potential\_vorticity.m}}: Load the field computed with
148 \mbox{{\ttfamily C\_comp...}} and divide it by $-\rho$ to obtain the
149 correct potential vorticity. Require the density field and after loading, overwrite the
150 file with prefix {\ttfamily PV} or {\ttfamily splPV}.
151 \item {\bf {\ttfamily compute\_density.m}}: Compute the density $\rho$ from the potential
152 temperature and the salinity fields.
153 \item {\bf {\ttfamily compute\_JFz.m}}: Compute the surface vertical PV flux due to
154 frictional processes. Requires the wind stress components, density, potential density
155 and Ekman layer depth (all of them, except the wind stress, may be computed with the
156 package), and produces one output file with the PV flux $J^F_z$ (see
157 Eq.~\ref{sec:diag:pv:eq15a}) and with {\ttfamily JFz} as a prefix.
158 \item {\bf {\ttfamily compute\_JBz.m}}: Compute the surface vertical PV flux due to
159 diabatic processes as:
160 \begin{eqnarray}
161 J^B_z &=& -\frac{f}{h}\frac{\alpha Q_{net}}{C_w} \nonumber
162 \end{eqnarray}
163 which is a simplified version of the full expression given in
164 Eq.~(\ref{sec:diag:pv:eq14a}). Requires the net surface heat flux and the mixed layer depth
165 (of which an estimation can be computed with the package), and produces one output file
166 with the PV flux $J^B_z$ and with {\ttfamily JBz} as a prefix.
167 \item {\bf {\ttfamily compute\_QEk.m}}: Compute the horizontal heat flux due to Ekman
168 currents from the PV flux induced by frictional forces as:
169 \begin{eqnarray}
170 Q_{Ek} &=& - \frac{C_w \delta_e}{\alpha f}J^F_z\nonumber
171 \end{eqnarray}
172 Requires the PV flux due to frictional forces and the Ekman layer depth, and produces one
173 output with the heat flux and with {\ttfamily QEk} as a prefix.
174 \item {\bf {\ttfamily eg\_main\_getPV}}: A complete example of how to set up a master
175 routine able to compute everything from the package.
176 \end{itemize}
177
178 %%%%%%%%%%%%%%%%%%%%%%%%%%
179 \subsection{Technical details}
180
181 \subsubsection{File name}
182 \noindent
183 A file name is formed by three parameters which need to be set up as global variables in Matlab
184 before running any routines. They are:
185 \begin{itemize}
186 \item the prefix, ie the variable name ({\ttfamily netcdf\_UVEL} for example). This
187 parameter is specified in the help section of all diagnostic routines.
188 \item {\ttfamily netcdf\_domain}: the geographical domain.
189 \item {\ttfamily netcdf\_suff}: the netcdf extension (nc or cdf for example).
190 \end{itemize}
191 Then, for example, if the calling Matlab routine had set up:
192 \begin{verbatim}
193 global netcdf_THETA netcdf_SALTanom netcdf_domain netcdf_suff
194 netcdf_THETA = 'THETA';
195 netcdf_SALTanom = 'SALT';
196 netcdf_domain = 'north_atlantic';
197 netcdf_suff = 'nc';
198 \end{verbatim}
199 the routine {\ttfamily A\_compute\_potential\_density.m} to compute the potential density
200 field, will look for the files:
201 \begin{verbatim}
202 THETA.north_atlantic.nc
203 SALT.north_atlantic.nc
204 \end{verbatim}
205 and the output file will automatically be: {\ttfamily SIGMATHETA.north\_atlantic.nc}.
206
207 Otherwise indicated, output file prefix cannot be changed.
208
209 \subsubsection{Path to file}
210 \noindent
211 All diagnostic routines look for input files in a subdirectory (relative to the Matlab
212 routine directory) called {\ttfamily ./netcdf-files} which in turn, is supposed to contain
213 subdirectories for each set of fields. For example, computing the potential density for
214 the timestep 12H00 02/03/2005 will require a
215 subdirectory with the potential temperature and salinity files like:
216 \begin{verbatim}
217 ./netcdf-files/200501031200/THETA.north_atlantic.nc
218 ./netcdf-files/200501031200/SALT.north_atlantic.nc
219 \end{verbatim}
220 \noindent
221 The output file {\ttfamily SIGMATHETA.north\_atlantic.nc} will be created in {\ttfamily
222 ./netcdf-files/200501031200/}.
223 \noindent
224 All diagnostic routines take as argument the name of the timestep subdirectory into
225 {\ttfamily ./netcdf-files}.
226
227 \subsubsection{Grids}
228 \noindent
229 With MITgcm numerical outputs, velocity and tracer fields may not be defined on the same
230 grid. Usually, UVEL and VVEL are defined on a C-grid but when interpolated from a
231 cube-sphere simulation they are defined on a A-grid. When it is needed, routines allow
232 to set up a global variable which define the grid to use.
233
234
235 %%%%%%%%%%%%%%%%%%%%%%%%%%
236 \subsection{Notes on the flux form of the PV equation and vertical PV fluxes}
237 \label{sec:diag:pv:notes-flux-form}
238
239 \subsubsection{Flux form of the PV equation}
240 The conservative flux form of the potential vorticity equation is:
241 \begin{eqnarray}
242 \frac{\partial \rho Q}{\partial t} + \nabla \cdot \vec{J} &=& 0 \label{sec:diag:pv:eq4}
243 \end{eqnarray}
244 where the potential vorticity $Q$ is given by the Eq.\ref{sec:diag:pv:eq2}.
245
246 The generalized flux vector of potential vorticity is:
247 \begin{eqnarray}
248 \vec{J} &=& \rho Q \vec{u} + \vec{N_Q}
249 \end{eqnarray}
250 which allows to rewrite Eq.\ref{sec:diag:pv:eq4} as:
251 \begin{eqnarray}
252 \frac{DQ}{dt} &=& -\frac{1}{\rho}\nabla\cdot\vec{N_Q} \label{sec:diag:pv:eq5}
253 \end{eqnarray}
254 where the nonadvective PV flux $\vec{N_Q}$ is given by:
255 \begin{eqnarray}
256 \vec{N_Q} &=& -\frac{\rho_0}{g}B\vec{\omega_a} + \vec{F}\times\nabla\sigma_\theta
257 \label{sec:diag:pv:eq6}
258 \end{eqnarray}
259
260 Its first component is linked to the buoyancy forcing\footnote{
261 Note that introducing B into Eq.\ref{sec:diag:pv:eq6} yields to:
262 \begin{eqnarray}
263 \vec{N_Q} &=& \omega_a \frac{D \sigma_\theta}{dt} + \vec{F}\times\nabla\sigma_\theta
264 \end{eqnarray}}:
265 \begin{eqnarray}
266 B &=& -\frac{g}{\rho_o}\frac{D \sigma_\theta}{dt}
267 \end{eqnarray}
268 and the second one to the nonconservative body forces per unit mass:
269 \begin{eqnarray}
270 \vec{F} &=& \frac{D \vec{u}}{dt} + 2\Omega\times\vec{u} + \nabla p
271 \end{eqnarray}
272
273 \subsubsection{Determining the PV flux at the ocean's surface}
274 In the context of mode water study, we're particularly interested in how the PV may be
275 reduced by surface PV fluxes because a mode water is characterised by a low PV level.
276 Considering the volume limited by two $iso-\sigma_\theta$, PV
277 flux is limited to surface processes and then vertical component of $\vec{N_Q}$. It is
278 supposed that $B$ and $\vec{F}$ will only be nonzero in the mixed layer (of depth $h$ and
279 variable density $\sigma_m$) exposed to mechanical forcing by the wind and buoyancy fluxes
280 through the ocean's surface.
281
282 Given the assumption of a mechanical forcing confined to a thin surface Ekman layer (of
283 depth $\delta_e$, eventually computed by the package) and of hydrostatic and geostrophic
284 balances, we can write:
285 \begin{eqnarray}
286 \vec{u_g} &=& \frac{1}{\rho f} \vec{k}\times\nabla p \\
287 \frac{\partial p_m}{\partial z} &=& -\sigma_m g \\
288 \frac{\partial \sigma_m}{\partial t} + \vec{u}_m\cdot\nabla\sigma_m &=& -\frac{\rho_0}{g}B \label{sec:diag:pv:eq7}
289 \end{eqnarray}
290 where:
291 \begin{eqnarray}
292 \vec{u}_m &=& \vec{u}_g + \vec{u}_{Ek} + o(R_o) \label{sec:diag:pv:eq8}
293 \end{eqnarray}
294 is the full velocity field composed by the geostrophic current $\vec{u}_g$ and the Ekman
295 drift:
296 \begin{eqnarray}
297 \vec{u}_{Ek} &=& -\frac{1}{\rho f}\vec{k}\times\frac{\partial \tau}{\partial z} \label{sec:diag:pv:eq9}
298 \end{eqnarray}
299 (where $\tau$ is the wind stress) and last by other ageostrophic components of $o(R_o)$
300 which are neglected.
301
302 Partitioning the buoyancy forcing as:
303 \begin{eqnarray}
304 B &=& B_g + B_{Ek} \label{sec:diag:pv:eq10}
305 \end{eqnarray}
306 and using Eq.\ref{sec:diag:pv:eq8} and Eq.\ref{sec:diag:pv:eq9}, the Eq.\ref{sec:diag:pv:eq7} becomes:
307 \begin{eqnarray}
308 \frac{\partial \sigma_m}{\partial t} + \vec{u}_g\cdot\nabla\sigma_m &=& -\frac{\rho_0}{g} B_g
309 \end{eqnarray}
310 revealing the "wind-driven buoyancy forcing":
311 \begin{eqnarray}
312 B_{Ek} &=& \frac{g}{\rho_0}\frac{1}{\rho f}\left(\vec{k}\times\frac{\partial \tau}{\partial z}\right)\cdot\nabla\sigma_m
313 \end{eqnarray}
314 Note that since:
315 \begin{eqnarray}
316 \frac{\partial B_g}{\partial z} &=& \frac{\partial}{\partial z}\left(-\frac{g}{\rho_0}\vec{u_g}\cdot\nabla\sigma_m\right)
317 = -\frac{g}{\rho_0}\frac{\partial \vec{u_g}}{\partial z}\cdot\nabla\sigma_m
318 = 0
319 \end{eqnarray}
320 $B_g$ must be uniform throughout the depth of the mixed layer and then being related to
321 the surface buoyancy flux by integrating Eq.\ref{sec:diag:pv:eq10} through the mixed layer:
322 \begin{eqnarray}
323 \int_{-h}^0B\,dz &=\, hB_g + \int_{-h}^0B_{Ek}\,dz \,=& \mathcal{B}_{in} \label{sec:diag:pv:eq11}
324 \end{eqnarray}
325 where $\mathcal{B}_{in}$ is the vertically integrated surface buoyancy (in)flux:
326 \begin{eqnarray}
327 \mathcal{B}_{in} &=& \frac{g}{\rho_o}\left( \frac{\alpha Q_{net}}{C_w} - \rho_0\beta S_{net}\right)
328 \label{sec:diag:pv:eq12}
329 \end{eqnarray}
330 with $\alpha\simeq 2.5\times10^{-4}\, K^{-1}$ the thermal expansion coefficient (computed
331 by the package otherwise), $C_w=4187J.kg^{-1}.K^{-1}$ the specific heat of seawater,
332 $Q_{net}[W.m^{-2}]$ the net heat surface flux (positive downward, warming the ocean),
333 $\beta[PSU^{-1}]$ the saline contraction coefficient, and
334 $S_{net}=S*(E-P)[PSU.m.s^{-1}]$ the net freshwater surface flux with $S[PSU]$
335 the surface salinity and $(E-P)[m.s^{-1}]$ the fresh water flux.
336
337 Introducing the body force in the Ekman layer:
338 \begin{eqnarray}
339 F_z &=& \frac{1}{\rho}\frac{\partial \tau}{\partial z}
340 \end{eqnarray}
341 the vertical component of Eq.\ref{sec:diag:pv:eq6} is:
342 \begin{eqnarray}
343 \vec{N_Q}_z &=& -\frac{\rho_0}{g}(B_g+B_{Ek})\omega_z
344 + \frac{1}{\rho}
345 \left( \frac{\partial \tau}{\partial z}\times\nabla\sigma_\theta \right)\cdot\vec{k}
346 \nonumber \\
347 &=& -\frac{\rho_0}{g}B_g\omega_z
348 -\frac{\rho_0}{g}
349 \left(\frac{g}{\rho_0}\frac{1}{\rho f}\vec{k}\times\frac{\partial \tau}{\partial z}
350 \cdot\nabla\sigma_m\right)\omega_z
351 + \frac{1}{\rho}
352 \left( \frac{\partial \tau}{\partial z}\times\nabla\sigma_\theta \right)\cdot\vec{k}
353 \nonumber \\
354 &=& -\frac{\rho_0}{g}B_g\omega_z
355 + \left(1-\frac{\omega_z}{f}\right)\left(\frac{1}{\rho}\frac{\partial \tau}{\partial z}
356 \times\nabla\sigma_\theta \right)\cdot\vec{k}
357 \end{eqnarray}
358 and given the assumption that $\omega_z\simeq f$, the second term vanishes and we obtain:
359 \begin{eqnarray}
360 \vec{N_Q}_z &=& -\frac{\rho_0}{g}f B_g \label{sec:diag:pv:eq12}
361 \end{eqnarray}
362 Note that the wind-stress forcing does not appear explicitly here but is implicit in $B_g$
363 through Eq.\ref{sec:diag:pv:eq11}: the buoyancy forcing $B_g$ is determined by the
364 difference between the integrated surface buoyancy flux $\mathcal{B}_{in}$ and the
365 integrated "wind-driven buoyancy forcing":
366 \begin{eqnarray}
367 B_g &=& \frac{1}{h}\left( \mathcal{B}_{in} - \int_{-h}^0B_{Ek}dz \right) \nonumber \\
368 &=& \frac{1}{h}\frac{g}{\rho_0}\left( \frac{\alpha Q_{net}}{C_w} - \rho_0 \beta S_{net}\right)
369 - \frac{1}{h}\int_{-h}^0
370 \frac{g}{\rho_0}\frac{1}{\rho f}\vec{k}\times \frac{\partial \tau}{\partial z} \cdot\nabla\sigma_m dz
371 \nonumber \\
372 &=& \frac{1}{h}\frac{g}{\rho_0}\left( \frac{\alpha Q_{net}}{C_w} - \rho_0 \beta S_{net}\right)
373 - \frac{g}{\rho_0}\frac{1}{\rho f \delta_e}\vec{k}\times\tau\cdot\nabla\sigma_m
374 \end{eqnarray}
375 Finally, from Eq.\ref{sec:diag:pv:eq6}, the vertical surface flux of PV may be written as:
376 \begin{eqnarray}
377 \vec{N_Q}_z &=& J^B_z + J^F_z \label{sec:diag:pv:eq13} \\
378 J^B_z &=& -\frac{f}{h}\left( \frac{\alpha Q_{net}}{C_w}-\rho_0 \beta S_{net}\right) \label{sec:diag:pv:eq14}\\
379 J^F_z &=& \frac{1}{\rho\delta_e} \vec{k}\times\tau\cdot\nabla\sigma_m \label{sec:diag:pv:eq15}
380 \end{eqnarray}

  ViewVC Help
Powered by ViewVC 1.1.22