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 6625 for branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters/Chap_LBC.tex – NEMO

Ignore:
Timestamp:
2016-05-26T11:08:07+02:00 (8 years ago)
Author:
kingr
Message:

Rolled back to r6613

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/DOC/TexFiles/Chapters/Chap_LBC.tex

    r6617 r6625  
    11% ================================================================ 
    2 % Chapter Lateral Boundary Condition (LBC)  
     2% Chapter Lateral Boundary Condition (LBC)  
    33% ================================================================ 
    44\chapter{Lateral Boundary Condition (LBC) } 
     
    204204%        North fold (\textit{jperio = 3 }to $6)$  
    205205% ------------------------------------------------------------------------------------------------------------- 
    206 \subsection{North-fold (\textit{jperio = 3 }to $6$) } 
     206\subsection{North-fold (\textit{jperio = 3 }to $6)$ } 
    207207\label{LBC_north_fold} 
    208208 
    209209The north fold boundary condition has been introduced in order to handle the north  
    210 boundary of a three-polar ORCA grid. Such a grid has two poles in the northern hemisphere  
    211 (Fig.\ref{Fig_MISC_ORCA_msh}, and thus requires a specific treatment illustrated in Fig.\ref{Fig_North_Fold_T}.  
    212 Further information can be found in \mdl{lbcnfd} module which applies the north fold boundary condition. 
     210boundary of a three-polar ORCA grid. Such a grid has two poles in the northern hemisphere.  
     211\colorbox{yellow}{to be completed...} 
    213212 
    214213%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    251250ocean model. Second order finite difference schemes lead to local discrete  
    252251operators that depend at the very most on one neighbouring point. The only  
    253 non-local computations concern the vertical physics (implicit diffusion,  
     252non-local computations concern the vertical physics (implicit diffusion, 1.5  
    254253turbulent closure scheme, ...) (delocalization over the whole water column),  
    255254and the solving of the elliptic equation associated with the surface pressure  
    256255gradient computation (delocalization over the whole horizontal domain).  
    257256Therefore, a pencil strategy is used for the data sub-structuration  
     257\gmcomment{no idea what this means!} 
    258258: the 3D initial domain is laid out on local processor  
    259259memories following a 2D horizontal topological splitting. Each sub-domain  
     
    264264phase starts: each processor sends to its neighbouring processors the update  
    265265values of the points corresponding to the interior overlapping area to its  
    266 neighbouring sub-domain ($i.e.$ the innermost of the two overlapping rows).  
    267 The communication is done through the Message Passing Interface (MPI).  
    268 The data exchanges between processors are required at the very  
     266neighbouring sub-domain (i.e. the innermost of the two overlapping rows).  
     267The communication is done through message passing. Usually the parallel virtual  
     268language, PVM, is used as it is a standard language available on  nearly  all  
     269MPP computers. More specific languages (i.e. computer dependant languages)  
     270can be easily used to speed up the communication, such as SHEM on a T3E  
     271computer. The data exchanges between processors are required at the very  
    269272place where lateral domain boundary conditions are set in the mono-domain  
    270 computation : the \rou{lbc\_lnk} routine (found in \mdl{lbclnk} module)  
    271 which manages such conditions is interfaced with routines found in \mdl{lib\_mpp} module  
    272 when running on an MPP computer ($i.e.$ when \key{mpp\_mpi} defined).  
    273 It has to be pointed out that when using the MPP version of the model,  
    274 the east-west cyclic boundary condition is done implicitly,  
    275 whilst the south-symmetric boundary condition option is not available. 
     273computation (\S III.10-c): the lbc\_lnk routine which manages such conditions  
     274is substituted by mpplnk.F or mpplnk2.F routine when running on an MPP  
     275computer (\key{mpp\_mpi} defined). It has to be pointed out that when using  
     276the MPP version of the model, the east-west cyclic boundary condition is done  
     277implicitly, whilst the south-symmetric boundary condition option is not available. 
    276278 
    277279%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    283285%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    284286 
    285 In the standard version of \NEMO, the splitting is regular and arithmetic. 
    286 The i-axis is divided by \jp{jpni} and the j-axis by \jp{jpnj} for a number of processors  
    287 \jp{jpnij} most often equal to $jpni \times jpnj$ (parameters set in  
    288  \ngn{nammpp} namelist). Each processor is independent and without message passing  
    289  or synchronous process, programs run alone and access just its own local memory.  
    290  For this reason, the main model dimensions are now the local dimensions of the subdomain (pencil)  
     287In the standard version of the OPA model, the splitting is regular and arithmetic. 
     288 the i-axis is divided by \jp{jpni} and the j-axis by \jp{jpnj} for a number of processors  
     289 \jp{jpnij} most often equal to $jpni \times jpnj$ (model parameters set in  
     290 \mdl{par\_oce}). Each processor is independent and without message passing  
     291 or synchronous process  
     292 \gmcomment{how does a synchronous process relate to this?},  
     293 programs run alone and access just its own local memory. For this reason, the  
     294 main model dimensions are now the local dimensions of the subdomain (pencil)  
    291295 that are named \jp{jpi}, \jp{jpj}, \jp{jpk}. These dimensions include the internal  
    292296 domain and the overlapping rows. The number of rows to exchange (known as  
     
    300304where \jp{jpni}, \jp{jpnj} are the number of processors following the i- and j-axis. 
    301305 
    302 One also defines variables nldi and nlei which correspond to the internal domain bounds,  
    303 and the variables nimpp and njmpp which are the position of the (1,1) grid-point in the global domain.  
    304 An element of $T_{l}$, a local array (subdomain) corresponds to an element of $T_{g}$,  
    305 a global array (whole domain) by the relationship:  
     306\colorbox{yellow}{Figure IV.3: example of a domain splitting with 9 processors and  
     307no east-west cyclic boundary conditions.} 
     308 
     309One also defines variables nldi and nlei which correspond to the internal  
     310domain bounds, and the variables nimpp and njmpp which are the position  
     311of the (1,1) grid-point in the global domain. An element of $T_{l}$, a local array  
     312(subdomain) corresponds to an element of $T_{g}$, a global array  
     313(whole domain) by the relationship:  
    306314\begin{equation} \label{Eq_lbc_nimpp} 
    307315T_{g} (i+nimpp-1,j+njmpp-1,k) = T_{l} (i,j,k), 
     
    312320nproc. In the standard version, a processor has no more than four neighbouring  
    313321processors named nono (for north), noea (east), noso (south) and nowe (west)  
    314 and two variables, nbondi and nbondj, indicate the relative position of the processor : 
     322and two variables, nbondi and nbondj, indicate the relative position of the processor  
     323\colorbox{yellow}{(see Fig.IV.3)}: 
    315324\begin{itemize} 
    316325\item       nbondi = -1    an east neighbour, no west processor, 
     
    323332processor on its overlapping row, and sends the data issued from internal  
    324333domain corresponding to the overlapping row of the other processor. 
     334        
     335\colorbox{yellow}{Figure IV.4: pencil splitting with the additional outer halos } 
    325336 
    326337 
     
    332343global ocean where more than 50 \% of points are land points. For this reason, a  
    333344pre-processing tool can be used to choose the mpp domain decomposition with a  
    334 maximum number of only land points processors, which can then be eliminated (Fig. \ref{Fig_mppini2}) 
    335 (For example, the mpp\_optimiz tools, available from the DRAKKAR web site).  
     345maximum number of only land points processors, which can then be eliminated.  
     346(For example, the mpp\_optimiz tools, available from the DRAKKAR web site.)  
    336347This optimisation is dependent on the specific bathymetry employed. The user  
    337348then chooses optimal parameters \jp{jpni}, \jp{jpnj} and \jp{jpnij} with  
    338349$jpnij < jpni \times jpnj$, leading to the elimination of $jpni \times jpnj - jpnij$  
    339 land processors. When those parameters are specified in \ngn{nammpp} namelist,  
     350land processors. When those parameters are specified in module \mdl{par\_oce},  
    340351the algorithm in the \rou{inimpp2} routine sets each processor's parameters (nbound,  
    341352nono, noea,...) so that the land-only processors are not taken into account.  
    342353 
    343 \gmcomment{Note that the inimpp2 routine is general so that the original inimpp  
     354\colorbox{yellow}{Note that the inimpp2 routine is general so that the original inimpp  
    344355routine should be suppressed from the code.} 
    345356 
    346357When land processors are eliminated, the value corresponding to these locations in  
    347 the model output files is undefined. Note that this is a problem for the meshmask file  
    348 which requires to be defined over the whole domain. Therefore, user should not eliminate  
    349 land processors when creating a meshmask file ($i.e.$ when setting a non-zero value to \np{nn\_msh}). 
     358the model output files is zero. Note that this is a problem for a mesh output file written  
     359by such a model configuration, because model users often divide by the scale factors  
     360($e1t$, $e2t$, etc) and do not expect the grid size to be zero, even on land. It may be  
     361best not to eliminate land processors when running the model especially to write the  
     362mesh files as outputs (when \np{nn\_msh} namelist parameter differs from 0). 
     363%% 
     364\gmcomment{Steven : dont understand this, no land processor means no output file  
     365covering this part of globe; its only when files are stitched together into one that you  
     366can leave a hole} 
     367%% 
    350368 
    351369%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    362380%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    363381 
     382 
     383% ================================================================ 
     384% Open Boundary Conditions  
     385% ================================================================ 
     386\section{Open Boundary Conditions (\key{obc}) (OBC)} 
     387\label{LBC_obc} 
     388%-----------------------------------------nam_obc  ------------------------------------------- 
     389%-    nobc_dta    =    0     !  = 0 the obc data are equal to the initial state 
     390%-                           !  = 1 the obc data are read in 'obc   .dta' files 
     391%-    rn_dpein      =    1.    !  ??? 
     392%-    rn_dpwin      =    1.    !  ??? 
     393%-    rn_dpnin      =   30.    !  ??? 
     394%-    rn_dpsin      =    1.    !  ??? 
     395%-    rn_dpeob      = 1500.    !  time relaxation (days) for the east  open boundary 
     396%-    rn_dpwob      =   15.    !    "        "           for the west  open boundary 
     397%-    rn_dpnob      =  150.    !    "        "           for the north open boundary 
     398%-    rn_dpsob      =   15.    !    "        "           for the south open boundary 
     399%-    ln_obc_clim = .true.   !  climatological obc data files (default T) 
     400%-    ln_vol_cst  = .true.   !  total volume conserved 
     401\namdisplay{namobc}  
     402 
     403It is often necessary to implement a model configuration limited to an oceanic  
     404region or a basin, which communicates with the rest of the global ocean through  
     405''open boundaries''. As stated by \citet{Roed1986}, an open boundary is a  
     406computational border where the aim of the calculations is to allow the perturbations  
     407generated inside the computational domain to leave it without deterioration of the  
     408inner model solution. However, an open boundary also has to let information from  
     409the outer ocean enter the model and should support inflow and outflow conditions.  
     410 
     411The open boundary package OBC is the first open boundary option developed in  
     412NEMO (originally in OPA8.2). It allows the user to  
     413\begin{itemize} 
     414\item tell the model that a boundary is ''open'' and not closed by a wall, for example  
     415by modifying the calculation of the divergence of velocity there; 
     416\item impose values of tracers and velocities at that boundary (values which may  
     417be taken from a climatology): this is the``fixed OBC'' option.  
     418\item calculate boundary values by a sophisticated algorithm combining radiation  
     419and relaxation (``radiative OBC'' option) 
     420\end{itemize} 
     421 
     422Options are defined through the \ngn{namobc} namelist variables. 
     423The package resides in the OBC directory. It is described here in four parts: the  
     424boundary geometry (parameters to be set in \mdl{obc\_par}), the forcing data at  
     425the boundaries (module \mdl{obcdta}),  the radiation algorithm involving the  
     426namelist and module \mdl{obcrad}, and a brief presentation of boundary update  
     427and restart files. 
     428 
     429%---------------------------------------------- 
     430\subsection{Boundary geometry} 
     431\label{OBC_geom} 
     432% 
     433First one has to realize that open boundaries may not necessarily be located  
     434at the extremities of the computational domain. They may exist in the middle  
     435of the domain, for example at Gibraltar Straits if one wants to avoid including  
     436the Mediterranean in an Atlantic domain. This flexibility has been found necessary  
     437for the CLIPPER project \citep{Treguier_al_JGR01}. Because of the complexity of the  
     438geometry of ocean basins, it may even be necessary to have more than one  
     439''west'' open boundary, more than one ''north'', etc. This is not possible with  
     440the OBC option: only one open boundary of each kind, west, east, south and  
     441north is allowed; these names refer to the grid geometry (not to the direction  
     442of the geographical ''west'', ''east'', etc). 
     443 
     444The open boundary geometry is set by a series of parameters in the module  
     445\mdl{obc\_par}. For an eastern open boundary, parameters are \jp{lp\_obc\_east}  
     446(true if an east open boundary exists), \jp{jpieob} the $i$-index along which  
     447the eastern open boundary (eob) is located, \jp{jpjed} the $j$-index at which  
     448it starts, and \jp{jpjef} the $j$-index where it ends (note $d$ is for ''d\'{e}but''  
     449and $f$ for ''fin'' in French). Similar parameters exist for the west, south and  
     450north cases (Table~\ref{Tab_obc_param}). 
     451 
     452 
     453%--------------------------------------------------TABLE-------------------------------------------------- 
     454\begin{table}[htbp]     \begin{center}    \begin{tabular}{|l|c|c|c|} 
     455\hline 
     456Boundary and  & Constant index  & Starting index (d\'{e}but) & Ending index (fin) \\ 
     457Logical flag  &                 &                            &                     \\ 
     458\hline 
     459West          & \jp{jpiwob} $>= 2$         &  \jp{jpjwd}$>= 2$          &  \jp{jpjwf}<= \np{jpjglo}-1 \\ 
     460lp\_obc\_west & $i$-index of a $u$ point   & $j$ of a $T$ point   &$j$ of a $T$ point \\ 
     461\hline 
     462East            & \jp{jpieob}$<=$\np{jpiglo}-2&\jp{jpjed} $>= 2$         & \jp{jpjef}$<=$ \np{jpjglo}-1 \\ 
     463 lp\_obc\_east  & $i$-index of a $u$ point    & $j$ of a $T$ point & $j$ of a $T$ point \\ 
     464\hline 
     465South           & \jp{jpjsob} $>= 2$         & \jp{jpisd} $>= 2$          & \jp{jpisf}$<=$\np{jpiglo}-1 \\ 
     466lp\_obc\_south  & $j$-index of a $v$ point   & $i$ of a $T$ point   & $i$ of a $T$ point \\ 
     467\hline 
     468North           & \jp{jpjnob} $<=$ \np{jpjglo}-2& \jp{jpind} $>= 2$        & \jp{jpinf}$<=$\np{jpiglo}-1 \\ 
     469lp\_obc\_north  & $j$-index of a $v$ point      & $i$  of a $T$ point & $i$ of a $T$ point \\ 
     470\hline 
     471\end{tabular}    \end{center} 
     472\caption{     \label{Tab_obc_param} 
     473Names of different indices relating to the open boundaries. In the case  
     474of a completely open ocean domain with four ocean boundaries, the parameters  
     475take exactly the values indicated.} 
     476\end{table} 
     477%------------------------------------------------------------------------------------------------------------ 
     478 
     479The open boundaries must be along coordinate lines. On the C-grid, the boundary  
     480itself is along a line of normal velocity points: $v$ points for a zonal open boundary  
     481(the south or north one), and $u$ points for a meridional open boundary (the west  
     482or east one). Another constraint is that there still must be a row of masked points  
     483all around the domain, as if the domain were a closed basin (unless periodic conditions  
     484are used together with open boundary conditions). Therefore, an open boundary  
     485cannot be located at the first/last index, namely, 1, \jp{jpiglo} or \jp{jpjglo}. Also,  
     486the open boundary algorithm involves calculating the normal velocity points situated  
     487just on the boundary, as well as the tangential velocity and temperature and salinity  
     488just outside the boundary. This means that for a west/south boundary, normal  
     489velocities and temperature are calculated at the same index \jp{jpiwob} and  
     490\jp{jpjsob}, respectively. For an east/north boundary, the normal velocity is  
     491calculated at index \jp{jpieob} and \jp{jpjnob}, but the ``outside'' temperature is  
     492at index \jp{jpieob}+1 and \jp{jpjnob}+1. This means that \jp{jpieob}, \jp{jpjnob}  
     493cannot be bigger than \jp{jpiglo}-2, \jp{jpjglo}-2.  
     494 
     495 
     496The starting and ending indices are to be thought of as $T$ point indices: in many  
     497cases they indicate the first land $T$-point, at the extremity of an open boundary  
     498(the coast line follows the $f$ grid points, see Fig.~\ref{Fig_obc_north} for an example  
     499of a northern open boundary). All indices are relative to the global domain. In the  
     500free surface case it is possible to have ``ocean corners'', that is, an open boundary  
     501starting and ending in the ocean. 
     502 
     503%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     504\begin{figure}[!t]     \begin{center} 
     505\includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_obc_north.pdf} 
     506\caption{    \label{Fig_obc_north} 
     507Localization of the North open boundary points.} 
     508\end{center}     \end{figure} 
     509%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     510 
     511Although not compulsory, it is highly recommended that the bathymetry in the  
     512vicinity of an open boundary follows the following rule: in the direction perpendicular  
     513to the open line, the water depth should be constant for 4 grid points. This is in  
     514order to ensure that the radiation condition, which involves model variables next  
     515to the boundary, is calculated in a consistent way. On Fig.\ref{Fig_obc_north} we  
     516indicate by an $=$ symbol, the points which should have the same depth. It means  
     517that at the 4 points near the boundary, the bathymetry is cylindrical \gmcomment{not sure  
     518why cylindrical}. The line behind the open $T$-line must be 0 in the bathymetry file  
     519(as shown on Fig.\ref{Fig_obc_north} for example). 
     520 
     521%---------------------------------------------- 
     522\subsection{Boundary data} 
     523\label{OBC_data} 
     524 
     525It is necessary to provide information at the boundaries. The simplest case is  
     526when this information does not change in time and is equal to the initial conditions  
     527(namelist variable \np{nn\_obcdta}=0). This is the case for the standard configuration  
     528EEL5 with open boundaries. When (\np{nn\_obcdta}=1), open boundary information  
     529is read from netcdf files. For convenience the input files are supposed to be similar  
     530to the ''history'' NEMO output files, for dimension names and variable names.  
     531Open boundary arrays must be dimensioned according to the parameters of table~ 
     532\ref{Tab_obc_param}: for example, at the western boundary, arrays have a  
     533dimension of \jp{jpwf}-\jp{jpwd}+1 in the horizontal and \jp{jpk} in the vertical.  
     534 
     535When ocean observations are used to generate the boundary data (a hydrographic  
     536section for example, as in \citet{Treguier_al_JGR01}) it happens often that only the velocity  
     537normal to the boundary is known, which is the reason why the initial OBC code  
     538assumes that only $T$, $S$, and the normal velocity ($u$ or $v$) needs to be  
     539specified. As more and more global model solutions and ocean analysis products  
     540become available, it will be possible to provide information about all the variables  
     541(including the tangential velocity) so that the specification of four variables at each  
     542boundaries will become standard. For the sea surface height, one must distinguish  
     543between the filtered free surface case and the time-splitting or explicit treatment of  
     544the free surface. 
     545 In the first case, it is assumed that the user does not wish to represent high  
     546 frequency motions such as tides. The boundary condition is thus one of zero  
     547 normal gradient of sea surface height at the open boundaries, following \citet{Marchesiello2001}.  
     548No information other than the total velocity needs to be provided at the open  
     549boundaries in that case. In the other two cases (time splitting or explicit free surface),  
     550the user must provide barotropic information (sea surface height and barotropic  
     551velocities) and the use of the Flather algorithm for barotropic variables is  
     552recommanded. However, this algorithm has not yet been fully tested and bugs  
     553remain in NEMO v2.3. Users should read the code carefully before using it. Finally,  
     554in the case of the rigid lid approximation the barotropic streamfunction must be  
     555provided, as documented in \citet{Treguier_al_JGR01}). This option is no longer  
     556recommended but remains in NEMO V2.3. 
     557 
     558One frequently encountered case is when an open boundary domain is constructed  
     559from a global or larger scale NEMO configuration. Assuming the domain corresponds  
     560to indices $ib:ie$, $jb:je$ of the global domain, the bathymetry and forcing of the  
     561small domain can be created by using the following netcdf utility on the global files:  
     562ncks -F $-d\;x,ib,ie$ $-d\;y,jb,je$ (part of the nco series of utilities,  
     563see their \href{http://nco.sourceforge.net}{website}).  
     564The open boundary files can be constructed using ncks  
     565commands, following table~\ref{Tab_obc_ind}.  
     566 
     567%--------------------------------------------------TABLE-------------------------------------------------- 
     568\begin{table}[htbp]     \begin{center}      \begin{tabular}{|l|c|c|c|c|c|} 
     569\hline 
     570OBC  & Variable   & file name      & Index  & Start  & end  \\ 
     571West &  T,S       &   obcwest\_TS.nc &  $ib$+1     &   $jb$+1 &  $je-1$  \\ 
     572     &    U       &   obcwest\_U.nc  &  $ib$+1     &   $jb$+1 &  $je-1$  \\  
     573     &    V       &   obcwest\_V.nc  &  $ib$+1     &   $jb$+1 &  $je-1$  \\        
     574\hline 
     575East &  T,S       &   obceast\_TS.nc &  $ie$-1     &   $jb$+1 &  $je-1$  \\ 
     576     &    U       &   obceast\_U.nc  &  $ie$-2     &   $jb$+1 &  $je-1$  \\  
     577     &    V       &   obceast\_V.nc  &  $ie$-1     &   $jb$+1 &  $je-1$  \\        
     578\hline          
     579South &  T,S      &   obcsouth\_TS.nc &  $jb$+1     &  $ib$+1 &  $ie-1$  \\ 
     580      &    U      &   obcsouth\_U.nc  &  $jb$+1     &  $ib$+1 &  $ie-1$  \\  
     581      &    V      &   obcsouth\_V.nc  &  $jb$+1     &  $ib$+1 &  $ie-1$  \\     
     582\hline 
     583North &  T,S      &   obcnorth\_TS.nc &  $je$-1     &  $ib$+1 &  $ie-1$  \\ 
     584      &    U      &   obcnorth\_U.nc  &  $je$-1     &  $ib$+1 &  $ie-1$  \\  
     585      &    V      &   obcnorth\_V.nc  &  $je$-2     &  $ib$+1 &  $ie-1$  \\   
     586\hline 
     587\end{tabular}     \end{center} 
     588\caption{    \label{Tab_obc_ind} 
     589Requirements for creating open boundary files from a global configuration,  
     590appropriate for the subdomain of indices $ib:ie$, $jb:je$. ``Index'' designates the  
     591$i$ or $j$ index along which the $u$ of $v$ boundary point is situated in the global  
     592configuration, starting and ending with the $j$ or $i$ indices indicated.  
     593For example, to generate file obcnorth\_V.nc, use the command ncks  
     594$-F$ $-d\;y,je-2$  $-d\;x,ib+1,ie-1$ }  
     595\end{table} 
     596%----------------------------------------------------------------------------------------------------------- 
     597 
     598It is assumed that the open boundary files contain the variables for the period of  
     599the model integration. If the boundary files contain one time frame, the boundary  
     600data is held fixed in time. If the files contain 12 values, it is assumed that the input  
     601is a climatology for a repeated annual cycle (corresponding to the case \np{ln\_obc\_clim}  
     602=true). The case of an arbitrary number of time frames is not yet implemented  
     603correctly; the user is required to write his own code in the module \mdl{obc\_dta}  
     604to deal with this situation.  
     605 
     606\subsection{Radiation algorithm} 
     607\label{OBC_rad} 
     608 
     609The art of open boundary management consists in applying a constraint strong  
     610enough that the inner domain "feels" the rest of the ocean, but weak enough 
     611that perturbations are allowed to leave the domain with minimum false reflections  
     612of energy. The constraints are specified separately at each boundary as time  
     613scales for ''inflow'' and ''outflow'' as defined below. The time scales are set (in days)  
     614by namelist parameters such as \np{rn\_dpein}, \np{rn\_dpeob} for the eastern open  
     615boundary for example. When both time scales are zero for a given boundary  
     616($e.g.$ for the western boundary, \jp{lp\_obc\_west}=true, \np{rn\_dpwob}=0 and  
     617\np{rn\_dpwin}=0) this means that the boundary in question is a ''fixed '' boundary  
     618where the solution is set exactly by the boundary data. This is not recommended,  
     619except in combination with increased viscosity in a ''sponge'' layer next to the  
     620boundary in order to avoid spurious reflections.   
     621 
     622 
     623The radiation\/relaxation \gmcomment{the / doesnt seem to appear in the output}  
     624algorithm is applied when either relaxation time (for ''inflow'' or ''outflow'') is  
     625non-zero. It has been developed and tested in the SPEM model and its  
     626successor ROMS \citep{Barnier1996, Marchesiello2001}, which is an  
     627$s$-coordinate model on an Arakawa C-grid. Although the algorithm has  
     628been numerically successful in the CLIPPER Atlantic models, the physics  
     629do not work as expected \citep{Treguier_al_JGR01}. Users are invited to consider  
     630open boundary conditions (OBC hereafter) with some scepticism  
     631\citep{Durran2001, Blayo2005}. 
     632 
     633The first part of the algorithm calculates a phase velocity to determine  
     634whether perturbations tend to propagate toward, or away from, the  
     635boundary. Let us consider a model variable $\phi$.  
     636The phase velocities ($C_{\phi x}$,$C_{\phi y}$) for the variable $\phi$,  
     637in the directions normal and tangential to the boundary are 
     638\begin{equation} \label{Eq_obc_cphi} 
     639C_{\phi x} = \frac{ -\phi_{t} }{ ( \phi_{x}^{2} + \phi_{y}^{2}) } \phi_{x}  
     640\;\;\;\;\; \;\;\;  
     641C_{\phi y} = \frac{ -\phi_{t} }{ ( \phi_{x}^{2} + \phi_{y}^{2}) } \phi_{y}.  
     642\end{equation} 
     643Following \citet{Treguier_al_JGR01} and \citet{Marchesiello2001} we retain only  
     644the normal component of the velocity, $C_{\phi x}$, setting $C_{\phi y} =0$  
     645(but unlike the original Orlanski radiation algorithm we retain $\phi_{y}$ in  
     646the expression for $C_{\phi x}$).   
     647 
     648The discrete form of (\ref{Eq_obc_cphi}), described by \citet{Barnier1998}, 
     649takes into account the two rows of grid points situated inside the domain  
     650next to the boundary, and the three previous time steps ($n$, $n-1$, 
     651and $n-2$). The same equation can then be discretized at the boundary at 
     652time steps $n-1$, $n$ and $n+1$ \gmcomment{since the original was three time-level}  
     653in order to extrapolate for the new boundary value $\phi^{n+1}$.  
     654 
     655In the open boundary algorithm as implemented in NEMO v2.3, the new boundary  
     656values are updated differently depending on the sign of $C_{\phi x}$. Let us take  
     657an eastern boundary as an example. The solution for variable $\phi$ at the  
     658boundary is given by a generalized wave equation with phase velocity $C_{\phi}$,  
     659with the addition of a relaxation term, as: 
     660\begin{eqnarray} 
     661\phi_{t} &  =  & -C_{\phi x} \phi_{x} + \frac{1}{\tau_{o}} (\phi_{c}-\phi)  
     662                        \;\;\; \;\;\; \;\;\; (C_{\phi x} > 0), \label{Eq_obc_rado} \\ 
     663\phi_{t} &  =  & \frac{1}{\tau_{i}} (\phi_{c}-\phi)  
     664\;\;\; \;\;\; \;\;\;\;\;\; (C_{\phi x} < 0), \label{Eq_obc_radi} 
     665\end{eqnarray} 
     666where $\phi_{c}$ is the estimate of $\phi$ at the boundary, provided as boundary  
     667data. Note that in (\ref{Eq_obc_rado}), $C_{\phi x}$ is bounded by the ratio  
     668$\delta x/\delta t$ for stability reasons. When $C_{\phi x}$ is eastward (outward  
     669propagation), the radiation condition (\ref{Eq_obc_rado}) is used.  
     670When  $C_{\phi x}$ is westward (inward propagation), (\ref{Eq_obc_radi}) is  
     671used with a strong relaxation to climatology (usually $\tau_{i}=\np{rn\_dpein}=$1~day). 
     672Equation (\ref{Eq_obc_radi}) is solved with a Euler time-stepping scheme. As a  
     673consequence, setting $\tau_{i}$ smaller than, or equal to the time step is equivalent  
     674to a fixed boundary condition. A time scale of one day is usually a good compromise  
     675which guarantees that the inflow conditions remain close to climatology while ensuring  
     676numerical stability.  
     677 
     678In  the case of a western boundary located in the Eastern Atlantic, \citet{Penduff_al_JGR00}  
     679have been able to implement the radiation algorithm without any boundary data,  
     680using persistence from the previous time step instead. This solution has not worked  
     681in other cases \citep{Treguier_al_JGR01}, so that the use of boundary data is recommended.  
     682Even in the outflow condition (\ref{Eq_obc_rado}), we have found it desirable to  
     683maintain a weak relaxation to climatology. The time step is usually chosen so as to  
     684be larger than typical turbulent scales (of order 1000~days \gmcomment{or maybe seconds?}). 
     685 
     686The radiation condition is applied to the model variables: temperature, salinity,  
     687tangential and normal velocities. For normal and tangential velocities, $u$ and $v$,  
     688radiation is applied with phase velocities calculated from $u$ and $v$ respectively.   
     689For the radiation of tracers, we use the phase velocity calculated from the tangential  
     690velocity in order to avoid calculating too many independent radiation velocities and  
     691because tangential velocities and tracers have the same position along the boundary  
     692on a C-grid.   
     693 
     694\subsection{Domain decomposition (\key{mpp\_mpi})} 
     695\label{OBC_mpp} 
     696When \key{mpp\_mpi} is active in the code, the computational domain is divided  
     697into rectangles that are attributed each to a different processor. The open boundary  
     698code is ``mpp-compatible'' up to a certain point. The radiation algorithm will not  
     699work if there is an mpp subdomain boundary parallel to the open boundary at the  
     700index of the boundary, or the grid point after (outside), or three grid points before  
     701(inside). On the other hand, there is no problem if an mpp subdomain boundary  
     702cuts the open boundary perpendicularly. These geometrical limitations must be  
     703checked for by the user (there is no safeguard in the code).   
     704The general principle for the open boundary mpp code is that loops over the open  
     705boundaries {not sure what this means} are performed on local indices (nie0,  
     706nie1, nje0, nje1 for an eastern boundary for instance) that are initialized in module  
     707\mdl{obc\_ini}. Those indices have relevant values on the processors that contain  
     708a segment of an open boundary. For processors that do not include an open  
     709boundary segment, the indices are such that the calculations within the loops are  
     710not performed. 
     711\gmcomment{I dont understand most of the last few sentences} 
     712  
     713Arrays of climatological data that are read from files are seen by all processors  
     714and have the same dimensions for all (for instance, for the eastern boundary,  
     715uedta(jpjglo,jpk,2)). On the other hand, the arrays for the calculation of radiation  
     716are local to each processor (uebnd(jpj,jpk,3,3) for instance).  This allowed the  
     717CLIPPER model for example, to save on memory where the eastern boundary  
     718crossed 8 processors so that \jp{jpj} was much smaller than (\jp{jpjef}-\jp{jpjed}+1).  
     719 
     720\subsection{Volume conservation} 
     721\label{OBC_vol} 
     722 
     723It is necessary to control the volume inside a domain when using open boundaries.  
     724With fixed boundaries, it is enough to ensure that the total inflow/outflow has  
     725reasonable values (either zero or a value compatible with an observed volume  
     726balance). When using radiative boundary conditions it is necessary to have a  
     727volume constraint because each open boundary works independently from the  
     728others. The methodology used to control this volume is identical to the one  
     729coded in the ROMS model \citep{Marchesiello2001}. 
     730 
     731 
     732%---------------------------------------- EXTRAS 
     733\colorbox{yellow}{Explain obc\_vol{\ldots}} 
     734 
     735\colorbox{yellow}{OBC algorithm for update, OBC restart, list of routines where obc key appears{\ldots}} 
     736 
     737\colorbox{yellow}{OBC rigid lid? {\ldots}} 
    364738 
    365739% ==================================================================== 
Note: See TracChangeset for help on using the changeset viewer.