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

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

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


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

1 cnh 1.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