New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 2268 for branches/DEV_r1826_DOC – NEMO

Ignore:
Timestamp:
2010-10-14T17:48:21+02:00 (14 years ago)
Author:
sga
Message:

NEMO branch DEV_r1826_DOC
Edit of first half of TRA chapter.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_TRA.tex

    r2195 r2268  
    4848In the present chapter we also describe the diagnostic equations used to compute  
    4949the sea-water properties (density, Brunt-Vais\"{a}l\"{a} frequency, specific heat and  
    50 freezing point) although the associated modules ($i.e.$ \mdl{eosbn2}, \mdl{ocfzpt}  
    51 and \mdl{phycst}) are (temporarily) located in the NEMO/OPA directory. 
    52  
    53 The different options available to the user are managed by namelist logical or  
     50freezing point with associated modules \mdl{eosbn2}, \mdl{ocfzpt} and \mdl{phycst}). 
     51 
     52The different options available to the user are managed by namelist logicals or  
    5453CPP keys. For each equation term \textit{ttt}, the namelist logicals are \textit{ln\_trattt\_xxx},  
    55 where \textit{xxx} is a 3 or 4 letter acronym accounting for each optional scheme.  
    56 The CPP key (when it exists) is \textbf{key\_trattt}. The corresponding code can be  
     54where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme.  
     55The CPP key (when it exists) is \textbf{key\_trattt}. The equivalent code can be  
    5756found in the \textit{trattt} or \textit{trattt\_xxx} module, in the NEMO/OPA/TRA directory. 
    5857 
     
    8887since the vertical scale factors are functions of $k$ only, and thus  
    8988$e_{3u} =e_{3v} =e_{3t} $. The flux form in \eqref{Eq_tra_adv}  
    90 requires implicitly the use of the continuity equation. Indeed, it is obtained 
     89implicitly requires the use of the continuity equation. Indeed, it is obtained 
    9190by using the following equality : $\nabla \cdot \left( \vect{U}\,T \right)=\vect{U} \cdot \nabla T$  
    9291which results from the use of the continuity equation, $\nabla \cdot \vect{U}=0$ or  
    93 $\partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$ in constant (default option)  
    94 or variable (\key{vvl} defined) volume case, respectively.  
     92$\partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$ in constant volume (default option)  
     93or variable volume (\key{vvl} defined) case, respectively.  
    9594Therefore it is of paramount importance to design the discrete analogue of the  
    9695advection tendency so that it is consistent with the continuity equation in order to  
    9796enforce the conservation properties of the continuous equations. In other words,  
    98 by substituting $\tau$ by 1 in (\ref{Eq_tra_adv}) we recover the discrete form of  
     97by replacing $\tau$ by the number 1 in (\ref{Eq_tra_adv}) we recover the discrete form of  
    9998the continuity equation which is used to calculate the vertical velocity. 
    10099%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    113112%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    114113 
    115 The key difference between the advection schemes used in \NEMO is the choice  
     114The key difference between the advection schemes available in \NEMO is the choice  
    116115made in space and time interpolation to define the value of the tracer at the  
    117116velocity points (Fig.~\ref{Fig_adv_scheme}).  
    118117 
    119 Along solid lateral and bottom boundaries a zero tracer flux is naturally  
     118Along solid lateral and bottom boundaries a zero tracer flux is automatically  
    120119specified, since the normal velocity is zero there. At the sea surface the  
    121120boundary condition depends on the type of sea surface chosen:  
     
    167166that are pure numerical artefacts.  
    168167 
     168\sgacomment{not sure whether 3 is still relevant after TRA-TRC branch is merged in?} 
    169169% ------------------------------------------------------------------------------------------------------------- 
    170170%        2nd order centred scheme   
     
    188188(\ref{Eq_tra_adv_cen2}) is the \textit{now} tracer value. The centered second  
    189189order advection is computed in the \mdl{traadv\_cen2} module. In this module, 
    190 it is also proposed to combine the \textit{cen2} scheme with an upstream scheme 
    191 in specific areas which requires a strong diffusion in order to avoid the generation  
     190it is advantageous to combine the \textit{cen2} scheme with an upstream scheme 
     191in specific areas which require a strong diffusion in order to avoid the generation  
    192192of false extrema. These areas are the vicinity of large river mouths, some straits  
    193193with coarse resolution, and the vicinity of ice cover area ($i.e.$ when the ocean  
    194194temperature is close to the freezing point). 
     195This combined scheme has been included for specific grid points in the ORCA2 and ORCA4 configurations only. 
    195196 
    196197Note that using the cen2 scheme, the overall tracer advection is of second  
    197198order accuracy since both (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_cen2})  
    198 have this order of accuracy. Note also that  
     199have this order of accuracy. \gmcomment{Note also that ... blah, blah} 
    199200 
    200201% ------------------------------------------------------------------------------------------------------------- 
     
    206207 
    207208In the $4^{th}$ order formulation (to be implemented), tracer values are  
    208 evaluated at velocity points as a $4^{th}$ order interpolation, and thus uses  
     209evaluated at velocity points as a $4^{th}$ order interpolation, and thus depend on  
    209210the four neighbouring $T$-points. For example, in the $i$-direction: 
    210211\begin{equation} \label{Eq_tra_adv_cen4} 
     
    257258\end{equation} 
    258259where $c_u$ is a flux limiter function taking values between 0 and 1.  
    259 There exist many ways to define $c_u$, each correcponding to a different  
     260There exist many ways to define $c_u$, each corresponding to a different  
    260261total variance decreasing scheme. The one chosen in \NEMO is described in  
    261262\citet{Zalesak_JCP79}. $c_u$ only departs from $1$ when the advective term  
     
    264265This scheme is tested and compared with MUSCL and the MPDATA scheme in  
    265266\citet{Levy_al_GRL01}; note that in this paper it is referred to as "FCT" (Flux corrected  
    266 transport) rather than TVD. The TVD scheme is computed in the \mdl{traadv\_tvd} module. 
    267  
    268 For stability reasons (see \S\ref{DOM_nxt}), in (\ref{Eq_tra_adv_tvd})  
    269 $\tau _u^{cen2}$ is evaluated using the \textit{now} tracer while $\tau _u^{ups}$  
     267transport) rather than TVD. The TVD scheme is implemented in the \mdl{traadv\_tvd} module. 
     268 
     269For stability reasons (see \S\ref{DOM_nxt}), 
     270$\tau _u^{cen2}$ is evaluated  in (\ref{Eq_tra_adv_tvd}) using the \textit{now} tracer while $\tau _u^{ups}$  
    270271is evaluated using the \textit{before} tracer. In other words, the advective part of  
    271272the scheme is time stepped with a leap-frog scheme while a forward scheme is  
     
    302303(\np{ln\_traadv\_muscl2}=true). Note that the latter choice does not ensure  
    303304the \textit{positive} character of the scheme. Only the former can be used  
    304 on both active and passive tracers. The two MUSCL schemes are computed  
     305on both active and passive tracers. The two MUSCL schemes are implemented  
    305306in the \mdl{traadv\_tvd} and \mdl{traadv\_tvd2} modules. 
    306307 
     
    339340\np{ln\_traadv\_ubs}=true. 
    340341 
    341 For stability reasons  (see \S\ref{DOM_nxt}), in \eqref{Eq_tra_adv_ubs},  
    342 the first term (which corresponds to a second order centred scheme)  
     342For stability reasons  (see \S\ref{DOM_nxt}), 
     343the first term  in \eqref{Eq_tra_adv_ubs} (which corresponds to a second order centred scheme)  
    343344is evaluated using the \textit{now} tracer (centred in time) while the  
    344345second term (which is the diffusive part of the scheme), is  
     
    352353substitution in the \mdl{traadv\_ubs} module and obtain a QUICK scheme. 
    353354 
     355Four different options are possible for the vertical  
     356component used in the UBS scheme. $\tau _w^{ubs}$ can be evaluated  
     357using either \textit{(a)} a centred $2^{nd}$ order scheme, or  \textit{(b)}  
     358a TVD scheme, or  \textit{(c)} an interpolation based on conservative  
     359parabolic splines following the \citet{Shchepetkin_McWilliams_OM05}  
     360implementation of UBS in ROMS, or  \textit{(d)} a UBS. The $3^{rd}$ case  
     361has dispersion properties similar to an eighth-order accurate conventional scheme. 
     362The current reference version uses method b) 
     363 
    354364Note that : 
    355365 
     
    361371in the current reference version.  
    362372 
    363 (2) In a forthcoming release four options will be available for the vertical  
    364 component used in the UBS scheme. $\tau _w^{ubs}$ will be evaluated  
    365 using either \textit{(a)} a centred $2^{nd}$ order scheme, or  \textit{(b)}  
    366 a TVD scheme, or  \textit{(c)} an interpolation based on conservative  
    367 parabolic splines following the \citet{Shchepetkin_McWilliams_OM05}  
    368 implementation of UBS in ROMS, or  \textit{(d)} a UBS. The $3^{rd}$ case  
    369 has dispersion properties similar to an eighth-order accurate conventional scheme. 
    370  
    371 (3) It is straightforward to rewrite \eqref{Eq_tra_adv_ubs} as follows: 
     373(2) It is straightforward to rewrite \eqref{Eq_tra_adv_ubs} as follows: 
    372374\begin{equation} \label{Eq_traadv_ubs2} 
    373375\tau _u^{ubs} = \tau _u^{cen4} + \frac{1}{12} \left\{   
     
    391393Thirdly, the diffusion term is in fact a biharmonic operator with an eddy  
    392394coefficient which is simply proportional to the velocity: 
    393  $A_u^{lm}= - \frac{1}{12}\,{e_{1u}}^3\,|u|$. Note that NEMO v2.3 still uses  
    394  \eqref{Eq_tra_adv_ubs}, not \eqref{Eq_traadv_ubs2}. This should be  
    395  changed in forthcoming release. 
     395 $A_u^{lm}= - \frac{1}{12}\,{e_{1u}}^3\,|u|$. Note that NEMO v3.3 still uses  
     396 \eqref{Eq_tra_adv_ubs}, not \eqref{Eq_traadv_ubs2}. 
    396397 %%% 
    397398 \gmcomment{the change in UBS scheme has to be done} 
     
    412413The resulting scheme is quite expensive but \emph{positive}.  
    413414It can be used on both active and passive tracers.  
    414 Nevertheless, the intrinsic diffusion of QCK makes its use risky in the vertical  
     415However, the intrinsic diffusion of QCK makes its use risky in the vertical  
    415416direction where the control of artificial diapycnal fluxes is of paramount importance.  
    416417Therefore the vertical flux is evaluated using the CEN2 scheme.  
    417 This no more ensure the positivity of the scheme. The use of TVD in the vertical  
    418 direction as for the UBS case should be implemented to maintain the property. 
     418This no longer guarantees the positivity of the scheme. The use of TVD in the vertical  
     419direction (as for the UBS case) should be implemented to restore this property. 
    419420 
    420421 
     
    427428 
    428429The Piecewise Parabolic Method (PPM) proposed by Colella and Woodward (1984)  
    429 is based on a quadradic piecewise rebuilding. Like the QCK scheme, it is associated  
     430\sgacomment{reference?} 
     431is based on a quadradic piecewise construction. Like the QCK scheme, it is associated  
    430432with the ULTIMATE QUICKEST limiter \citep{Leonard1991}. It has been implemented  
    431433in \NEMO by G. Reffray (MERCATOR-ocean) but is not yet offered in the reference  
    432 version 3.0. 
     434version 3.3. 
    433435 
    434436% ================================================================ 
     
    469471\end{equation} 
    470472where  $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$  is the volume of $T$-cells.  
    471 It can be found in the \mdl{traadv\_lap} module. 
     473It is implemented in the \mdl{traadv\_lap} module. 
    472474 
    473475This lateral operator is computed in \mdl{traldf\_lap}. It is a \emph{horizontal}  
    474476operator ($i.e.$ acting along geopotential surfaces) in the $z$-coordinate with  
    475 or without partial step, but is simply an iso-level operator in the $s$-coordinate.  
     477or without partial steps, but is simply an iso-level operator in the $s$-coordinate.  
    476478It is thus used when, in addition to \np{ln\_traldf\_lap}=true, we have  
    477 \np{ln\_traldf\_level}=true, or \np{ln\_traldf\_hor}=\np{ln\_zco}=true.  
     479\np{ln\_traldf\_level}=true or \np{ln\_traldf\_hor}=\np{ln\_zco}=true.  
    478480In both cases, it significantly contributes to diapycnal mixing.  
    479481It is therefore not recommended. 
    480482 
    481483Note that  
    482 (a) In pure $z$-coordinate (\key{zco} is defined), $e_{3u}$=$e_{3v}$=$e_{3t}$,  
     484(a) In the pure $z$-coordinate (\key{zco} is defined), $e_{3u}$=$e_{3v}$=$e_{3t}$,  
    483485so that the vertical scale factors disappear from (\ref{Eq_tra_ldf_lap}) ;  
    484 (b) In partial step $z$-coordinate (\np{ln\_zps}=true), tracers in horizontally  
     486(b) In the partial step $z$-coordinate (\np{ln\_zps}=true), tracers in horizontally  
    485487adjacent cells are located at different depths in the vicinity of the bottom.  
    486488In this case, horizontal derivatives in (\ref{Eq_tra_ldf_lap}) at the bottom level  
     
    541543background horizontal diffusion \citep{Guilyardi_al_CD01}. An alternative scheme  
    542544developed by \cite{Griffies_al_JPO98} which preserves both tracer and its variance  
    543 is currently been tested in \NEMO. It should be available in a forthcoming 
    544 release. 
     545is also available in \NEMO (\np{ln\_traldf_grif}=true). 
    545546 
    546547Note that in the partial step $z$-coordinate (\np{ln\_zps}=true), the horizontal  
     
    572573where $a=\sqrt{|A|}$ and $A<0$. This was a mistake: both formulations  
    573574ensure the total variance decrease, but the former requires a larger  
    574 number of code-lines. It will be corrected in a forthcoming release. 
     575number of code-lines. 
    575576 
    576577% ------------------------------------------------------------------------------------------------------------- 
     
    584585applying (\ref{Eq_tra_ldf_iso}) twice. It requires an additional assumption  
    585586on boundary conditions: first and third derivative terms normal to the  
    586 coast, the bottom and the surface are set to zero. It can be found in the 
     587coast, normal to the bottom and normal to the surface are set to zero. It can be found in the 
    587588\mdl{traldf\_bilapg}. 
    588589 
    589590It is used when, in addition to \np{ln\_traldf\_bilap}=true, we have  
    590591\np{ln\_traldf\_iso}= .true, or both \np{ln\_traldf\_hor}=true and \np{ln\_zco}=true.  
    591 Nevertheless, this rotated bilaplacian operator has never been seriously  
    592 tested. No warranties that it is neither free of bugs or correctly formulated.  
     592This rotated bilaplacian operator has never been seriously  
     593tested. There are no guarantees that it is either free of bugs or correctly formulated.  
    593594Moreover, the stability range of such an operator will be probably quite  
    594 narrow, requiring a significantly smaller time-step than the one used on  
     595narrow, requiring a significantly smaller time-step than the one used with an 
    595596unrotated operator. 
    596597 
     
    662663the thickness of the top model layer.  
    663664 
    664 Due to interactions and mass exchange with other media ($i.e.$ atmosphere, sea-ice, land), 
     665Due to interactions and mass exchange of water ($F_{mass}$) with other Earth system components ($i.e.$ atmosphere, sea-ice, land), 
    665666the change in the heat and salt content of the surface layer of the ocean is due both  
    666 to the heat and salt fluxes crossing the sea surface and not linked with $F_{mass}$, the water  
    667 exchange with the other media, and to the heat and salt content of this water exchange. 
    668 In a forcoming release, these two parts, computed in the surface module (SBC), will included directly 
     667to the heat and salt fluxes crossing the sea surface (not linked with $F_{mass}$) 
     668 and to the heat and salt content of the mass exchange. 
     669\sgacomment{ the following does not apply to the release to which this documentation is  
     670attached and so should not be included .... 
     671In a forthcoming release, these two parts, computed in the surface module (SBC), will be included directly 
    669672in $Q_{ns}$, the surface heat flux and $F_{salt}$, the surface salt flux. 
    670 \gmcomment{  The specification of those fluxes is further detail in the SBC chapter (see \S\ref{SBC}).  } 
    671 This change will provide a same forcing formulation for any tracers (including temperature and salinity). 
     673The specification of these fluxes is further detailed in the SBC chapter (see \S\ref{SBC}).  
     674This change will provide a forcing formulation which is the same for any tracer (including temperature and salinity). 
    672675  
    673 In the current version, the situation is a little bit more complicated.  
     676In the current version, the situation is a little bit more complicated. } 
    674677The surface module (\mdl{sbcmod}, see \S\ref{SBC}) provides the following  
    675678forcing fields (used on tracers): 
    676679 
    677 $\bullet$ $Q_{ns}$,the non solar part of the net surface heat flux that cross the sea surface  
    678 (difference between the total surface heat flux and the fraction of the short wave flux that  
     680$\bullet$ $Q_{ns}$, the non-solar part of the net surface heat flux that crosses the sea surface  
     681(i.e. the difference between the total surface heat flux and the fraction of the short wave flux that  
    679682penetrates into the water column, see \S\ref{TRA_qsr}) 
    680683 
    681684$\bullet$ \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation) 
    682685 
    683 $\bullet$ $\textit{emp}_S$, an equivalent mass flux taking into account the effect of ice-ocean mass exchanged 
     686$\bullet$ $\textit{emp}_S$, an equivalent mass flux taking into account the effect of ice-ocean mass exchange 
    684687 
    685688$\bullet$ \textit{rnf}, the mass flux associated with runoff (see \S\ref{SBC_rnf} for further detail of how it acts on temperature and salinity tendencies) 
    686689 
    687 The $\textit{emp}_S$ field is not simply the budget evaporation-precipitation+freezing-melting because  
    688 the sea-ice is not currently embedded in the ocean but levitates above it. There is not mass 
     690The $\textit{emp}_S$ field is not simply the budget of evaporation-precipitation+freezing-melting because  
     691the sea-ice is not currently embedded in the ocean but levitates above it. There is no mass 
    689692exchanged between the sea-ice and the ocean. Instead we only take into account the salt 
    690 flux link to the fact that sea-ice has a non-sero salinity, and the concentration/dilution effect 
     693flux associated with the non-zero salinity of sea-ice, and the concentration/dilution effect 
    691694due to the freezing/melting (F/M) process. These two parts of the forcing are then converted into  
    692 a equivalent mass flux given by $\textit{emp}_S - \textit{emp}$. As a result of this mess,  
     695an equivalent mass flux given by $\textit{emp}_S - \textit{emp}$. As a result of this mess,  
    693696the surface boundary condition on temperature and salinity is applied as follows: 
    694697 
    695 In the nonlinear free surface case (\key{vvl} is defined, \jp{lk\_vvl}=true): 
     698In the nonlinear free surface case (\key{vvl} is defined): 
    696699\begin{equation} \label{Eq_tra_sbc} 
    697700\begin{aligned} 
     
    704707\end{equation}  
    705708 
    706 In the linear free surface case (\key{vvl} not defined, , \jp{lk\_vvl}=false): 
     709In the linear free surface case (\key{vvl} not defined): 
    707710\begin{equation} \label{Eq_tra_sbc_lin} 
    708711\begin{aligned} 
     
    713716 \end{aligned} 
    714717\end{equation}  
    715 where $\overline{x }^t$ means that $x$ is averaged over two consecutive time step  
    716 ($t-\rdt/2$ and $t+\rdt/2$). Such a time averaged prevents the excitation of the  
     718where $\overline{x }^t$ means that $x$ is averaged over two consecutive time steps  
     719($t-\rdt/2$ and $t+\rdt/2$). Such time averaging prevents the  
    717720divergence of odd and even time step (see \S\ref{STP}). 
    718721 
    719722The two set of equations, \eqref{Eq_tra_sbc} and \eqref{Eq_tra_sbc_lin}, are obtained  
    720723by assuming that the temperature of precipitation and evaporation are equal to 
    721 the ocean surface temperature while their salinity is zero. Therefore, the heat content 
    722 of \textit{emp} budget must be added to the temperature equation in variable volume case,  
    723 while it does not appear in constant volume. Similarly, the \textit{emp} budget affects  
    724 the ocean surface salinity in constant volume case (through the concentration dilution effect) 
    725 while it does not appears explicitly in variable volume as salinity change will be 
    726 induced by volume change. In both constant and variable volume, surface salinity  
    727 will change with ice-ocean salt flux and F/M flux without mass exchanges  
    728 ($\textit{emp}_S - \textit{emp}$). 
    729  
    730 Note that concentration/dilution effect due to F/M is computed using 
     724the ocean surface temperature and that their salinity is zero. Therefore, the heat content 
     725of the \textit{emp} budget must be added to the temperature equation in the variable volume case,  
     726while it does not appear in the constant volume case. Similarly, the \textit{emp} budget affects  
     727the ocean surface salinity in the constant volume case (through the concentration dilution effect) 
     728while it does not appears explicitly in the variable volume case since salinity change will be 
     729induced by volume change. In both constant and variable volume cases, surface salinity  
     730will change with ice-ocean salt flux and F/M flux (both contained in $\textit{emp}_S - \textit{emp}$) without mass exchanges. 
     731 
     732Note that the concentration/dilution effect due to F/M is computed using 
    731733a constant ice salinity as well as a constant ocean salinity.  
    732734This approximation suppresses the correlation between \textit{SSS}  
     
    734736Indeed, if this approximation is not made, even if the F/M budget is zero  
    735737on average over the whole ocean domain and over the seasonal cycle,  
    736 the associated salt flux is not, since sea-surface salinity and F/M flux are  
     738the associated salt flux is not zero, since sea-surface salinity and F/M flux are  
    737739intrinsically correlated (high \textit{SSS} are found where freezing is  
    738 strong whilst low \textit{SSS} is usually associated with high melting areas. 
     740strong whilst low \textit{SSS} is usually associated with high melting areas). 
    739741 
    740742Even using this approximation, an exact conservation of heat and salt content  
    741743is only achieved in the variable volume case. In the constant volume case,  
    742 there is a small unbalance associated with the product $(\partial_t\eta - \textit{emp}) * \textit{SSS}$. 
     744there is a small imbalance associated with the product $(\partial_t\eta - \textit{emp}) * \textit{SSS}$. 
    743745Nevertheless, the salt content variation is quite small and will not induce 
    744 a long term drift as there is no physical reason that $(\partial_t\eta - \textit{emp})$  
    745 and \textit{SSS} are correlated \citep{Roullet_Madec_JGR00}.  
    746 Note that, while quite small, the unbalance in constant volume case is larger  
    747 than the unbalance associated with the Asselin time filter \citep{Leclair_Madec_OM09}.  
    748 This is the reason why the modified filter is not applied in constant volume case. 
     746a long term drift as there is no physical reason for $(\partial_t\eta - \textit{emp})$  
     747and \textit{SSS} to be correlated \citep{Roullet_Madec_JGR00}.  
     748Note that, while quite small, the imbalance in the constant volume case is larger  
     749than the imbalance associated with the Asselin time filter \citep{Leclair_Madec_OM09}.  
     750This is the reason why the modified filter is not applied in the constant volume case. 
    749751 
    750752% ------------------------------------------------------------------------------------------------------------- 
     
    759761 
    760762When the penetrative solar radiation option is used (\np{ln\_flxqsr}=true),  
    761 the solar radiation penetrates the top few 10 meters of the ocean, otherwise  
    762 all the heat flux is absorbed in the first ocean level (\np{ln\_flxqsr}=false).  
     763the solar radiation penetrates the top few tens of meters of the ocean. If it is not used  
     764(\np{ln\_flxqsr}=false) all the heat flux is absorbed in the first ocean level.  
    763765Thus, in the former case a term is added to the time evolution equation of  
    764 temperature \eqref{Eq_PE_tra_T} whilst the surface boundary condition is  
     766temperature \eqref{Eq_PE_tra_T} and the surface boundary condition is  
    765767modified to take into account only the non-penetrative part of the surface  
    766768heat flux: 
     
    780782The shortwave radiation,  $Q_{sr}$, consists of energy distributed across a wide spectral range.  
    781783The ocean is strongly absorbing for wavelengths longer than 700~nm and these  
    782 wavelengths contribute to heating the upper few 10 centimetres. The fraction of $Q_{sr}$  
     784wavelengths contribute to heating the upper few tens of centimetres. The fraction of $Q_{sr}$  
    783785that resides in these almost non-penetrative wavebands, $R$, is $\sim 58\%$ (specified  
    784786through namelist parameter \np{rn\_abs}).  It is assumed to penetrate the ocean  
    785 following a decreasing exponential profile, with an e-folding depth scale, $\xi_0$,  
    786 of a few 10 centimetres (typically $\xi_0=0.35~m$ set as \np{rn\_si0} in the namtra\_qsr namlist). 
     787with a decreasing exponential profile, with an e-folding depth scale, $\xi_0$,  
     788of a few tens of centimetres (typically $\xi_0=0.35~m$ set as \np{rn\_si0} in the namtra\_qsr namelist). 
    787789For shorter wavelengths (400-700~nm), the ocean is more transparent, and solar energy  
    788 propagates to depths where it contributes to a penetrating flux of solar energy and thus  
    789 to local heating below the surface.  
    790 The way this second part of the solar energy penetrates in the ocean depends on  
    791 which formulation is chosen. In the simple 2-wavebands light penetration  (\np{ln\_qsr\_2bd}=true)  
    792 a chlorophyll-independent monochromatic formulation is also chosen for the shorter wavelengths,  
     790propagates to larger depths where it contributes to  
     791local heating.  
     792The way this second part of the solar energy penetrates into the ocean depends on  
     793which formulation is chosen. In the simple 2-waveband light penetration scheme  (\np{ln\_qsr\_2bd}=true)  
     794a chlorophyll-independent monochromatic formulation is chosen for the shorter wavelengths,  
    793795leading to the following expression  \citep{Paulson1977}: 
    794796\begin{equation} \label{Eq_traqsr_iradiance} 
    795797I(z) = Q_{sr} \left[Re^{-z / \xi_0} + \left( 1-R\right) e^{-z / \xi_1} \right] 
    796798\end{equation} 
    797 where $\xi_1$ is the second extinction length scales associated with the shorter wavebands.   
    798 It is usually chosen to be 23~m through \np{rn\_si0} namelist parameter.  
     799where $\xi_1$ is the second extinction length scale associated with the shorter wavelengths.   
     800It is usually chosen to be 23~m by setting the \np{rn\_si0} namelist parameter.  
    799801The set of default values ($\xi_0$, $\xi_1$, $R$) corresponds to a Type I water in  
    800802Jerlov's (1968) classification (oligotrophic waters). 
     
    802804Such assumptions have been shown to provide a very crude and simplistic  
    803805representation of observed light penetration profiles (\cite{Morel_JGR88}, see also  
    804 Fig.\ref{Fig_traqsr_irradiance}). Light absorption in the ocean depends on the  
    805 particules concentration and it is spectrally selective. \cite{Morel_JGR88} has shown  
     806Fig.\ref{Fig_traqsr_irradiance}). Light absorption in the ocean depends on  
     807particle concentration and is spectrally selective. \cite{Morel_JGR88} has shown  
    806808that an accurate representation of light penetration can be provided by a 61 waveband  
    807809formulation. Unfortunately, such a model is very computationally expensive.  
    808810Thus, \cite{Lengaigne_al_CD07} have constructed a simplified version of this  
    809 formulation in which visible light is splitted into three wavebands: blue (400-500 nm),  
    810 green (500-600 nm) and red (600-700nm). For each wave-band, the chlorophyll-dependant  
     811formulation in which visible light is split into three wavebands: blue (400-500 nm),  
     812green (500-600 nm) and red (600-700nm). For each wave-band, the chlorophyll-dependent  
    811813attenuation coefficient is fitted to the coefficients computed from the full spectral model  
    812 of \cite{Morel_JGR88} (as modified by \cite{Morel_Maritorena_JGR01}) assuming  
    813 the same power-law expression. As shown on Fig.\ref{Fig_traqsr_irradiance},  
    814 this formulation, called RGB (Reed-Green-Blue), reproduces quite closely  
    815 the light penetration profiles predicted by the full spectal model with much faster  
    816 computing efficiently, in contrast with the 2-bands formulation.  
     814of \cite{Morel_JGR88} (as modified by \cite{Morel_Maritorena_JGR01}), assuming  
     815the same power-law relationship. As shown in Fig.\ref{Fig_traqsr_irradiance},  
     816this formulation, called RGB (Red-Green-Blue), reproduces quite closely  
     817the light penetration profiles predicted by the full spectal model, but with much greater  
     818computational efficiency. The 2-bands formulation does not reproduce the full model very well.  
    817819 
    818820The RGB formulation is used when \np{ln\_qsr\_rgb}=true. The RGB attenuation coefficients 
    819 ($i.e.$ the inverse of the extinction length scales) are tabulated over 61 nonuniform  
     821($i.e.$ the inverses of the extinction length scales) are tabulated over 61 nonuniform  
    820822chlorophyll classes ranging from 0.01 to 10 g.Chl/L (see the routine \rou{trc\_oce\_rgb}  
    821 in \mdl{trc\_oce} module). Three type of chlorophyll can be used in the RGB formulation: 
    822 (1) a constant 0.05 g.Chl/L value everywhere (\np{nn\_chdta}=0) ; (2) observed  
    823 time varying chlorophyll (\np{nn\_chdta}=0) ; (3) simulated time varying chlorophyll 
    824 by TOP biogeochemical model (\np{ln\_qsr\_bio}=true). In the later case, the RGB  
    825 formulation is used to calculated both the phytoplankton light limitation in PISCES  
     823in \mdl{trc\_oce} module). Three types of chlorophyll can be chosen in the RGB formulation: 
     824(1) a constant 0.05 g.Chl/L value everywhere (\np{nn\_chdta}=0) ; (2) an observed  
     825time varying chlorophyll (\np{nn\_chdta}=1) ; (3) simulated time varying chlorophyll 
     826by TOP biogeochemical model (\np{ln\_qsr\_bio}=true). In the latter case, the RGB  
     827formulation is used to calculate both the phytoplankton light limitation in PISCES  
    826828or LOBSTER and the oceanic heating rate.  
    827829 
     
    829831is added to the temperature trend, and the surface heat flux is modified in routine \mdl{traqsr}.  
    830832 
    831 When $z$-coordinate is preferred to $s$-coordinate, the depth of $w-$levels does  
     833When the $z$-coordinate is preferred to the $s$-coordinate, the depth of $w-$levels does  
    832834not significantly vary with location. The level at which the light has been totally  
    833835absorbed ($i.e.$ it is less than the computer precision) is computed once,  
    834 and the trend associated with the penetration of the solar radiation is only added until that level.  
     836and the trend associated with the penetration of the solar radiation is only added down to that level.  
    835837Finally, note that when the ocean is shallow ($<$ 200~m), part of the  
    836838solar radiation can reach the ocean floor. In this case, we have  
     
    841843\begin{figure}[!t] \label{Fig_traqsr_irradiance}  \begin{center} 
    842844\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_TRA_Irradiance.pdf} 
    843 \caption{Penetration profile of the Downward solar irradiance  
    844 calculated by four models. Two wavebands chlorophyll-independant formulation (blue),  
    845 a chlorophyll-dependant monochromatic formulation (green), 4 waveband RGB formulation (red),  
     845\caption{Penetration profile of the downward solar irradiance  
     846calculated by four models. Two waveband chlorophyll-independent formulation (blue),  
     847a chlorophyll-dependent monochromatic formulation (green), 4 waveband RGB formulation (red),  
    84684861 waveband Morel (1988) formulation (black) for a chlorophyll concentration of  
    847849(a) Chl=0.05 mg/m$^3$ and (b) Chl=0.5 mg/m$^3$. From \citet{Lengaigne_al_CD07}.} 
     
    873875earth cooling. This flux is weak compared to surface fluxes (a mean  
    874876global value of $\sim0.1\;W/m^2$ \citep{Stein_Stein_Nat92}), but it is  
    875 systematically positive and acts on the densest water masses.  
     877systematically positive \sgacomment{positive definite?} and acts on the densest water masses.  
    876878Taking this flux into account in a global ocean model increases 
    877879the deepest overturning cell ($i.e.$ the one associated with the Antarctic  
    878880Bottom Water) by a few Sverdrups  \citep{Emile-Geay_Madec_OS09}.  
    879881 
    880 The presence or not of geothermal heating is controlled by the namelist  
     882The presence of geothermal heating is controlled by the namelist  
    881883parameter  \np{nn\_geoflx}. When this parameter is set to 1, a constant  
    882884geothermal heating is introduced whose value is given by the  
     
    904906and volume flux of the densest waters of the ocean, such as Antarctic Bottom Water,  
    905907or North Atlantic Deep Water. $z$-coordinate models tend to overestimate the  
    906 entrainment, because the gravity flow is mixed down vertically by convection  
     908entrainment, because the gravity flow is mixed vertically by convection  
    907909as it goes ''downstairs'' following the step topography, sometimes over a thickness  
    908910much larger than the thickness of the observed gravity plume. A similar problem  
    909 occurs in the $s$-coordinate when the thickness of the bottom level varies in large  
    910 proportions downstream of a sill \citep{Willebrand_al_PO01}, and the thickness  
     911occurs in the $s$-coordinate when the thickness of the bottom level varies rapidly  
     912downstream of a sill \citep{Willebrand_al_PO01}, and the thickness  
    911913of the plume is not resolved.  
    912914 
     
    914916\citet{Beckmann_Doscher1997}, is to allow a direct communication between  
    915917two adjacent bottom cells at different levels, whenever the densest water is  
    916 located above the less dense water. The communication can be by a diffusive  
    917 (diffusive BBL), advective fluxes (advective BBL), or both. In the current  
     918located above the less dense water. The communication can be by a diffusive flux 
     919(diffusive BBL), an advective flux (advective BBL), or both. In the current  
    918920implementation of the BBL, only the tracers are modified, not the velocities.  
    919921Furthermore, it only connects ocean bottom cells, and therefore does not include  
    920 one of the improvment introduced by \citet{Campin_Goosse_Tel99}. 
     922the improvements introduced by \citet{Campin_Goosse_Tel99}. 
    921923 
    922924% ------------------------------------------------------------------------------------------------------------- 
     
    927929 
    928930When applying sigma-diffusion (\key{trabbl} defined and \np{nn\_bbl\_ldf} set to 1),  
    929 the diffusive flux between two adjacent cells living at the ocean bottom is given by  
     931the diffusive flux between two adjacent cells at the ocean floor is given by  
    930932\begin{equation} \label{Eq_tra_bbl_diff} 
    931933{\rm {\bf F}}_\sigma=A_l^\sigma \; \nabla_\sigma T 
     
    933935with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells,  
    934936and  $A_l^\sigma$ the lateral diffusivity in the BBL. Following \citet{Beckmann_Doscher1997},  
    935 the latter is prescribed with a spatial dependence, $e.g.$ in the conditional form 
     937the latter is prescribed with a spatial dependence, $i.e.$ in the conditional form 
    936938\begin{equation} \label{Eq_tra_bbl_coef} 
    937939A_l^\sigma (i,j,t)=\left\{ {\begin{array}{l} 
     
    943945where $A_{bbl}$ is the BBL diffusivity coefficient, given by the namelist  
    944946parameter \np{rn\_ahtbbl} and usually set to a value much larger  
    945 than the one used on lateral mixing in open ocean. The constraint in \eqref{Eq_tra_bbl_coef}  
    946 implies that sigma-like diffusion only occurs when density above the sea floor, at the top of  
     947than the one used for lateral mixing in the open ocean. The constraint in \eqref{Eq_tra_bbl_coef}  
     948implies that sigma-like diffusion only occurs when the density above the sea floor, at the top of  
    947949the slope, is larger than in the deeper ocean (see green arrow in Fig.\ref{Fig_bbl}).  
    948950In practice, this constraint is applied separately in the two horizontal directions,  
     
    951953   \nabla_\sigma \rho / \rho = \alpha \,\nabla_\sigma T + \beta   \,\nabla_\sigma S 
    952954\end{equation}  
    953 where $\rho$, $\alpha$ and $\beta$ are function of $\overline{T}^\sigma$,  
    954 $\overline{S}^\sigma$, $\overline{H}^\sigma$, the along bottom mean temperature,  
     955where $\rho$, $\alpha$ and $\beta$ are functions of $\overline{T}^\sigma$,  
     956$\overline{S}^\sigma$ and $\overline{H}^\sigma$, the along bottom mean temperature,  
    955957salinity and depth, respectively. 
    956958 
     
    961963\label{TRA_bbl_adv} 
    962964 
     965\sgacomment{"downsloping flow" has been replaced by "downslope flow" in the following 
     966if this is not what is meant then "downwards sloping flow" is also a possibility"} 
    963967 
    964968%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    965969\begin{figure}[!t] \label{Fig_bbl}  \begin{center} 
    966970\includegraphics[width=0.7\textwidth]{./TexFiles/Figures/Fig_BBL_adv.pdf} 
    967 \caption{Advective/diffusive Bottom Boundary Layer. The BBl parameterisation is  
     971\caption{Advective/diffusive Bottom Boundary Layer. The BBL parameterisation is  
    968972activated when $\rho^i_{kup}$ is larger than $\rho^{i+1}_{kdnw}$.  
    969973Red arrows indicate the additional overturning circulation due to the advective BBL.  
    970 The transport of the downsloping flow is defined either as the transport of the bottom  
     974The transport of the downslope flow is defined either as the transport of the bottom  
    971975ocean cell (black arrow), or as a function of the along slope density gradient.  
    972 The green arrow indicates the diffusive BBL flux connecting directly $kup$ and $kdwn$ 
     976The green arrow indicates the diffusive BBL flux directly connecting $kup$ and $kdwn$ 
    973977ocean bottom cells. 
    974978connection} 
     
    988992water to move down the slope.  
    989993 
    990 \np{nn\_bbl\_adv} = 1 : the downsloping velocity is chosen to be the Eulerian 
     994\np{nn\_bbl\_adv} = 1 : the downslope velocity is chosen to be the Eulerian 
    991995ocean velocity just above the topographic step (see black arrow in Fig.\ref{Fig_bbl})  
    992996\citep{Beckmann_Doscher1997}. It is a \textit{conditional advection}, that is, advection 
     
    995999greater depth ($i.e.$ $\vect{U}  \cdot  \nabla H>0$). 
    9961000 
    997 \np{nn\_bbl\_adv} = 2 : the downsloping velocity is chosen to be proportional to $\Delta \rho$, 
    998 the density difference between the top and down cell densities \citep{Campin_Goosse_Tel99}. 
     1001\np{nn\_bbl\_adv} = 2 : the downslope velocity is chosen to be proportional to $\Delta \rho$, 
     1002the density difference between the higher cell and lower cell densities \citep{Campin_Goosse_Tel99}. 
    9991003The advection is allowed only  if dense water overlies less dense water on the slope ($i.e.$  
    10001004$\nabla_\sigma \rho  \cdot  \nabla H<0$). For example, the resulting transport of the  
    1001 downsloping flow, here in the $i$-direction (Fig.\ref{Fig_bbl}), is simply given by the  
     1005downslope flow, here in the $i$-direction (Fig.\ref{Fig_bbl}), is simply given by the  
    10021006following expression: 
    10031007\begin{equation} \label{Eq_bbl_Utr} 
    10041008 u^{tr}_{bbl} = \gamma \, g \frac{\Delta \rho}{\rho_o}  e_{1u} \; min \left( {e_{3u}}_{kup},{e_{3u}}_{kdwn} \right) 
    10051009\end{equation} 
    1006 where $\gamma$, expressed in second, is the coefficient of proportionality  
     1010where $\gamma$, expressed in seconds, is the coefficient of proportionality  
    10071011provided as \np{rn\_gambbl}, a namelist parameter, and \textit{kup} and \textit{kdwn}  
    1008 are the vertical index of the top and bottom cells, respectively. 
     1012are the vertical index of the higher and lower cells, respectively. 
    10091013The parameter $\gamma$ should take a different value for each bathymetric  
    1010 step. But, for simplicity, and because no direct estimation of this parameter is  
    1011 available, a uniform value has been retained. The possible values for $\gamma$  
     1014step, but for simplicity, and because no direct estimation of this parameter is  
     1015available, a uniform value has been assumed. The possible values for $\gamma$  
    10121016range between 1 and $10~s$ \citep{Campin_Goosse_Tel99}.   
    10131017 
    1014 The scalar properties are advected by this additional transport $( u^{tr}_{bbl}, v^{tr}_{bbl} )$  
     1018Scalar properties are advected by this additional transport $( u^{tr}_{bbl}, v^{tr}_{bbl} )$  
    10151019using the upwind scheme. Such a diffusive advective scheme has been chosen  
    1016 to mimic the entrainment between the downsloping plume and the surrounding  
    1017 water at intermediate depth. The entrainment is replaced by the vertical mixing  
    1018 included in the advection scheme. Let us consider as an example the  
    1019 case display in Fig.\ref{Fig_bbl} where the density at level $(i,kup)$ is  
     1020to mimic the entrainment between the downslope plume and the surrounding  
     1021water at intermediate depths. The entrainment is replaced by the vertical mixing  
     1022implicit in the advection scheme. Let us consider as an example the  
     1023case displayed in Fig.\ref{Fig_bbl} where the density at level $(i,kup)$ is  
    10201024larger than the one at level $(i,kdwn)$. The advective BBL scheme 
    10211025modifies the tracer time tendency of the ocean cells near the  
    1022 topographic step by the downsloping flow \eqref{Eq_bbl_dw},  
     1026topographic step by the downslope flow \eqref{Eq_bbl_dw},  
    10231027the horizontal \eqref{Eq_bbl_hor}  and the upward \eqref{Eq_bbl_up}  
    10241028return flows as follows:  
     
    10621066are given temperature and salinity fields (usually a climatology).  
    10631067The restoring term is added when \key{tradmp} is defined.  
    1064 It also requires that both \key{temdta} and \key{saldta} are defined  
     1068It also requires that both \key{dtatem} and \key{dtasal} are defined  
    10651069($i.e.$ that $T_o$ and $S_o$ are read). The restoring coefficient  
    1066 $S_o$ is a three-dimensional array initialized by the user in routine  
     1070$\gamma$ is a three-dimensional array initialized by the user in routine  
    10671071\rou{dtacof} also located in module \mdl{tradmp}.  
    10681072 
     
    10741078field for a passive tracer study). The first case applies to regional  
    10751079models that have artificial walls instead of open boundaries.  
    1076 In the vicinity of these walls, $S_o$ takes large values (equivalent to  
     1080In the vicinity of these walls, $\gamma$ takes large values (equivalent to  
    10771081a time scale of a few days) whereas it is zero in the interior of the  
    10781082model domain. The second case corresponds to the use of the robust  
    10791083diagnostic method \citep{Sarmiento1982}. It allows us to find the velocity  
    1080 field consistent with the model dynamics whilst having a $T$-$S$ field  
    1081 close to a given climatological field ($T_o -S_o$). The time scale  
     1084field consistent with the model dynamics whilst having a $T$, $S$ field  
     1085close to a given climatological field ($T_o$, $S_o$). The time scale  
    10821086associated with $S_o$ is generally not a constant but spatially varying  
    10831087in order to respect other properties. For example, it is usually set to zero  
     
    10851089\citep{Madec_al_JPO96} and in the equatorial region  
    10861090\citep{Reverdin1991, Fujio1991, Marti_PhD92} since these two regions  
    1087 have a short time scale of adjustment; while smaller $S_o$ are used  
     1091have a short time scale of adjustment; while smaller $\gamma$ are used  
    10881092in the deep ocean where the typical time scale is long \citep{Sarmiento1982}.  
    10891093In addition the time scale is reduced (even to zero) along the western  
     
    10911095structure in equilibrium with its physics.  
    10921096The choice of the shape of the Newtonian damping is controlled by two  
    1093 namelist parameters \np{nn\_zdmp}. The former allows to specified the  
    1094 width of the equatorial band in which no damping is applied as well as a decrease  
    1095 in vicinity of the coast and a damping everywhere in the Red and Med Seas,  
    1096 whereas the latter set a damping acting in the mixed layer or not.  
     1097namelist parameters \np{??} and \np{nn\_zdmp}. The former allows us to specify: the  
     1098width of the equatorial band in which no damping is applied; a decrease  
     1099in the vicinity of the coast; and a damping everywhere in the Red and Med Seas. 
     1100The latter sets whether damping should act in the mixed layer or not.  
    10971101The time scale associated with the damping depends on the depth as 
    10981102a hyperbolic tangent, with \np{rn\_surf} as surface value, \np{rn\_bot} as  
Note: See TracChangeset for help on using the changeset viewer.