Numerical modeling of Atmosphere and Ocean

[vinu]

Many of us are familiar in using Re-analysis data sets from
various sources. A few may be aware of how these are modelled
in real time. A lot of effort has been made since 1940 (the early
stages of computer) to model the atmospher and ocean. Since then,
thoughts about 'Numerical Schemes with stable solutions'
were also involved.

Here you can download a technical report by Global
Atmospheric Research Programme- (GARP) published in early 1970's.
Any of todays models can not exclude these basic laws of modeling
structured by Messinger and Arakawa (1970's). Most of us may be aware
of NCEP/NCAR, ECMWF, GFDL models and many others. All these models
are designed according to these 'basic' rules
developed in early 1970's.


Numerical Methods used in Atmospheric Models (Mesinger and Arakawa,1976),
GARP publication series, No.17. (69 pages, 14 MB) (http://wwwoa.ees.hokudai.ac.jp/people/vinu/share_pal/book.pdf)

PS: If the file or link does not work properly, please feed back.

----------------------------------------------------------------------------------------------------

[vinu]

This is a short note on Numerical Ocean Model developement.
A Good and Short note to understand the basics about
Ocean Modeling like,

What is an Ocean Model
What are the governing equations of a simple ocean model
What are the various ocean models usually in practice
What is a grid and how they are designed
How to arrange these governing equations in a grid
How to integrate it in time for a forecast.

read more in Ocean Model (PDF file ) -74 pages (1.4 MB) (http://wwwoa.ees.hokudai.ac.jp/people/vinu/ocean_model.pdf)

Source: http://csep1.phy.ornl.gov/om/om.html
Copyright (C) 1991, 1992, 1993, 1994, 1995 by the
Computational Science Education Project

-----------------------------------------------------------------------------------------------------

[Suprit]

I have browsed through the modelling book and its informative indeed, thanks.

A valuable book for the fellows interesting in modelling ( although not available free on the net and not so popular in oceanography student circles but certainly a must for mets ) is "Numerical prediction and dynamic meteorology" by Haltiner & Williams published by wiley. The concepts are dealt from the basics, and in straightforward manner.

[vinu]

Books which are worth refering to understand
about the basics of Ocean models-

Numerical Models of Oceans and Oceanic Processes
by L.H. Kantha, C.A. Clayson, International Geophysics Series,
Vol.66, 2000.
(-This book combines the fundementals of GFD and
a good start point for a beginner. It include examples for simple
models and expanded it in the finite difference form which will be
useful in understanding the development of a model-)

Numerical Modeling of Ocean Dynamics
by Z. Kowalik, T. S. Murty, Advanced Series on Ocean Engineering,
World Scientific, 1993
(-This book points to the essentials of advanced Numerical
methods in Ocean Modeling-).

Numerical Ocean Circulation Modeling
by Dale B. Haidvogel, Aike Beckmann, Series on Environmental
Science and Management. Vol.2, Imperial College Press, 1999.
(-This is a short book wich introduces Numerical methods in
Ocean Modeling and discuss the various OGCM usually under practice-).


Small Scale Processess in Geophysical Fluid Flows
by L. H. Kantha, C. A. Clayson, International Geophysics Series,
Vol.67, 2000.
(-This book is usefull to understand more about upper layer
(boundary layer) processes. It encorporate the modeling of Chemical
and Biological component in the Ocean-)

[vinu]

Here in this post, let us see a simple Hydro-dynamical Ocean Model.

What is a model?

Roughly speaking, it is a 'structure' to represent something. For example, a model of a building which is made up of 'card-board'- is a physical model. Now in our discipline, an ocean process (like its motion, heating, mixing etc.) can be modeled using some governing equations', under which these processes are happening. For example, Miss. A has 'X' income and 'Y' expenditure. Then her savings is, 'S = X-Y'. This is a simple 'model' for 'savings of Miss. A'. Likewise, every process in the ocean is also under some basic 'physical balance' and it can be reasonably represented using some sets of equations. This is our key.

Here let us see a simple hydrodynamical ocean model which is usually called, under practice, as 1&1/2 Layer model .

[b]What is a 1&1/2 layer model?

Here the ocean is assumed into two layer- an upper layer and lower layer with slightly different density. However only the upper layer is 'under motion'. The 'lower layer' is very deep so that it is assumed to be at 'rest'. This is called 1&1/2 layer (means 'one' active layer and the rest 'quiescent layer'). These layers represent the 'well-mixed' upper layer of the real ocean and relatively slow and inactive 'deep-ocean'. This means that 1&1/2 layer version of ocean model can well represent a real ocean and that is its importance. (An upper layer thickness of 100m is a good practice).

The figure represents the structure of a 1&1/2 layer model with some necessary dimensional representations.

http://wwwoa.ees.hokudai.ac.jp/people/vinu/share_pal/model/model.fig.gif

The upper layer density is 'rho-1' and lower layer density is'rho-2' the upper layer velocities are 'U'(zonal) and 'V'(meridional) and lower layer velocities are zero . The 'free-surface' elevation is 'eta' and 'inter-facial' elevation is 'h'. The upper layer thickness is 'H' and lower layer thickness is 'Hd'.

If the lower layer is under 'rest', why is it important? Why the entire bottom layer is not taken as an 'ocean floor?

The lower layer interacts with the upper layer through the density stratification at the interface. Eventhough the lower layer is under rest, the difference in density at the interface causes internal waves to develop. And thus it affects the 'upper-layer'.

What are the governing equations in such situation?

The main process to model is its motion. Let us use the shallow-water
wave equation for the upper layer
http://wwwoa.ees.hokudai.ac.jp/people/vinu/share_pal/model/model.gif

The equation (1) and (2) are the 'momentum equations'. The L.H.S represents the 'evolution of velocity ' in time. The first term in the R.H.S represents the 'coriolis acceleration'. The second term in the R.H.S is the 'pressure-gradient' term. Equation (3) is the 'continuity equation', which means that the 'evolution' of the layer thickness (L.H.S) is balanced by the 'horizontal divergence' (R.H.S) of the layer.

All the terms in the equation (1) and (2) are known for the upper layer. However in (3) we need 'h', which is the 'inter-facial' displacement due to the 'density stratification'.

How can we model 'h' (the inter-facial displacement)?

Since the lower layer is in rest, it implies that, there is no pressure gradient between two points in the 'lower layer' in order to drive the motion. Let us use this tool and try to find 'h' consider the pressure at 'A' (shown in Figure-1). Using the simple hydrostatic balance (p=h.rho.g), the pressure at A is the sum of the pressure of the water column above it. From a careful examination of the figure we can estimate it as shown in Figure-2 as 'Pa'. Now take another point Pb on the same level, and subtract Pa-Pb. According to above constraint (i.e. pressure gradient is zero) Pa-Pb=0. This will lead us to get the equation (4).
Now we can find all the terms in the R.H.S of equation 1 to 3.

Now look at the equation (4). Usually in the ocean, the lower layer is heavier(denser) than the upper layer, resulting a negative value of (rho-1 - rho-2)in equation (4). This implies that, the 'free surface' elevation (eta) is in opposite 'direction' with the interface-displacement (h).

How do these equations work (or how does it model the ocean)?

As mentioned earlier, the savings of Miss. A is 'S = X-Y'. Well, this is a useless equation, because it doesnt tells us what is the savings of Miss. A last month, or next month. Suppose we have the 'rate-of-savings' of Miss. A, (something like, dS/dt = (X-Y), where 't' is the time.) Let us 'integrate this equation', with the knowledge of todays 'X and Y'. [b]PS: S(t2) denotes S_suffix_t2, means value of S at time t2 i.e. (S(t2) - S(t1) )/(t2-t1) = X-Y (in finite difference form).

i.e. S (t2) = S (t1) + (X-Y)*(t2-t1) , where 't1'= today, 't2' = tomorrow.

Now we end-up in S (t2) , that is Miss. A's tomorrows of savings,(provided todays rate is 'dS/dt'). This is how 'a forecast' is being done. The L.H.S of equation 1 to 3 is a 'time-derivative'. We know that if we 'integrate' a time derivative in the present time, we will get values for the future-time' (depending upon the time-interval we used in the integration).

So what?

Now we have to think about the way-of-integration . At first we need a 'space' to do the integration. Simply, we need to design our ocean in a regular 'grid' which represents some 'points (location)' in the real ocean. And then we have to express our equations in a 'finite difference form' for integrating'numerically'. (Just like what we have seen in Miss. A's example). (A good note about stable numerical schemes is here (http://www.oceanographers.net/forums/showthread.php?p=98)).

What is a Grid?

Roughly speaking, it is a group of points representing the ocean. We use these points to 'integrate' the above sets of equations. We can not take as many points as we want, because it needs laborious computational power. So let us restrict our 'points' necessary enough to represent the physical processes of the ocean in which we are interested in. Likewise we can-not restrict our points 'to-a-very-minimum' number to save the computational time because we may lose useful physical details of the ocean. There are well defined criteria, though they are not the concern of this post.

What type of grid we use?

In equation 1 to 3, we have three 'prognostic variables', i.e. variables to find out for 'future time'. They are 'eta or h, U and V. There are certain rules to define these variables in certain points of a grid in the ocean. Depending on various architecture, there are a wide variety of grids that are in practice. Let us see three of them.

A-Grid.
http://wwwoa.ees.hokudai.ac.jp/people/vinu/share_pal/model/a-grid.gif

It is a simple grid. The variable 'eta or h, U, V' are defined at the 'same location' in the ocean. (see figure above)

B-Grid.
http://wwwoa.ees.hokudai.ac.jp/people/vinu/share_pal/model/b-grid.gif

It is a 'staggered grid'. The 'U, V locations are at the corner points of the grid, whereas 'eta or h' is defined at the center of the grid.(see figure above)

C-Grid
http://wwwoa.ees.hokudai.ac.jp/people/vinu/share_pal/model/c-grid.gif

It is also a 'staggered grid'. The 'U' locations are 'half-grid size away' to the 'west of 'eta'. The 'V' locations are 'half-grid size away' to the south of 'eta'. (See figure above).

Why these jargons?

The staggering of the grids is related with the 'computational stability of various numerical schemes'. For designing appropriate finite difference scheme in space we need to 'locate' these variables in suitable positions so that the resulting scheme yields maximum stability in calculation. A note about this is included here. (http://www.oceanographers.net/forums/showthread.php?p=98)



Which grid is appropriate?

This depends upon the problem we are interested in. For the above set of 3 equations, any of the above grid can be used. Each has its own advantages and disadvantages. For example, if we want to solve 'U' component of velocity at a 'grid-location', we need 'dP/dX and 'f.V'. If we set a C-grid (suppose), we have two P (or eta) points on the right and left side of 'U'. So it is easy to take dP/dx in this case. However, we have no 'V' values at the same position of 'U' (see figure C-grid for understanding the situation). So we have to average the 'V' ('f' as well) from the neighboring corners to find the 'f.V'. This is not precise, but theres no other choice. This example tells the various technical details of a particular grid. If we use A-grid instead, all the values are available at the same location. However the resulting finite difference scheme is not necessarily computationally stable. Concluding these technical side (not discussing much in this post), let us set for a C-grid which is more appropriate.

Now, we have a Grid-, then?

In the Miss.A's saving example, we have seen a step of 'integration'. If we use the time-step of integration as 'one-month', we will get the savings of Miss. A for one month at once. (In a single step calculation). However, from the common logic we can follow that it is not-accurate, because the saving rate 'dS/dt ' may be changing from 'day to day'. Likewise in numerical modeling of geophysical flows, the selection of 'time-step' for the integration needs much more attention.

CFL limit

In numerical integration of Geophysical Fluids, it is necessary to follow a constraint in the time-step of the integration. The C.F.L condition says that (Courant-Friedrichs-Lewy condition) the time step chosen, dt < dx/C, where 'dx' is the grid size and 'C' is the phase speed of 'fastest wave' in the domain. In Geophysical Fluid flow, the fastest wave is the 'gravity wave' with speed = square-root(g.H), where g=acceleration due to gravity and 'H' is the fluid depth. For a 'barotropic' mode (the entire fluid column is assumed to be in uniform speed) H is usually 4000m. This yields the 'gravity wave' speed c = 200m/sec, which is very large. In this case for example, for a grid-size that is of 1 degree, the CFL limit says dt < (dx/C) , which comes around 600 seconds. So this is the time-step we can grant for our integration. We have to forecast all the variable in a short-future of 600 sec (or less). From there we do the next integration to get the next future values and so on.

However in 1&1/2 Layer version of ocean models, there is no 'barotropic' mode because we are concerned only with the motion of the'upper-layer' (lower layer is in rest). As well as the depth of the fluid which will be in some order of 100's, so the time-step can be effectively increased.

What is next?

Now we have a set of equations, we have grid, we have finite-difference form, we have chosen a 'time-step' for integration. Now comes the initial-condition as well as the boundary-conditions. In simple words, initial conditions means, the state of the ocean we refer to the model initially. The consequent solutions in time (or evolution in time) will be depending on this initial state of the ocean. Thus an accurate initial condition is necessary to yield a reliable evolution of the model in time.

For example, we can assume, the ocean has a particular Temperature and Salinity values initially. This may be selected from the available climatological values (Data). For the momentum equations, we can assume, for example, Ocean is initially at rest. That means zero velocity. In our above 1&1/2 layer, we are concerned only with the momentum equations. Sea-Surface elevation (eta) can be assumed to be flat initially so that there is no pressure-gradient across the surface and hence no-motion.

Now let us see what is meant by a boundary condition. The boundary means the portion, where the model is separated from its surroundings. Simply speaking, we can expect boundaries at the surface, bottom and the lateral walls. Sometimes the lateral walls may be open to some other ocean basin. From our simple logic we can expect that the water is not allowed (or should not) to flow out of these boundaries (except for open boundaries where the ocean is connected to another ocean basins or channels).

Thus generally we can classify the boundary condition as Closed boundary conditions and Open boundary conditions. For example the lateral wall as well as the bottom of the sea-floor is the closed boundaries whereas the connection of a semi-enclosed bay to the open ocean is an open boundary. At the closed boundary (i.e. for lateral walls) the velocity that is normal to the wall is assumed as zero (flow cannot penetrate the wall). The velocity parallel to the wall can be either set to non-zero or set to zero. In the non-zero case it is called 'free-slip' boundary condition in-which the water does not feel any frictional force of the wall. The zero-velocity case is called as'no-slip means that the friction within the grid close to the lateral wall is important. These two are ideal cases whereas a good practice for general exercises in ocean modeling. (The description about open boundaries are not mentioning here).


(...Continued) (http://www.oceanographers.net/forums/showthread.php?p=329)

If any of the above statements found misleading PLEASE FEED BACK.
- - - - - - - - - - - - - - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - -- - - - - - - - - - - -

[vinu]

The following is Praveen's commend over the above post, I quote it for him here,
-----------------------------------------------------------------------------------------------
under C grid picture, it's description

there "east" of eta is making problem, is it "west" ??

is it like this?
The 'U' locations are 'half grid size away'


g'=g.(rho2-rho1)/rho2

what is g', is it reduced gravity something like that?

the thread is much informative. try to include different cordinates
systems like sigma cordinates,......

with regards
-------------------------------------------------------------------------------------------------

Dear Praveen

Thanks for the commend and type mistake you noticed in the
above post. I have corrected those.

yes, the term g' is the 'reduced gravity'.
However, for the 1&1/2 layer version, g'
need not necessarily find during the calculation, because if you
substitute g' in 'eta' equation, you will reach equation (4) by
eliminating the 'reduced gravity'.

Now, If we apply a Rigid-Lid (i.e. ocean surface is believed to as rest,
or just like a solid slab, then we can eliminate 'eta' from the momentum
equation, and can be replaced by the 'interfacial displacement h'
from the equation (4). As well as, the 'eta' term in equation (3) can
be neglected comparing with the 'interfacial displacement 'h' '. Then
these equations will give the motion under the 'reduced gravity' .Actualy
this is what being done in a 1&1/2 layer model. This way, we can avoid the
fastest moving gravity wave with speed sqrt(g.H) and reduces it to sqrt(g'.H).
This will provide an effective way to increase the time-step for
integration as with the C.F.L limit.


thanks for your feed back

PS: this time i quote your commend, you can post your replies and notices
in the thread itslef. we meant for this 'Forum Structure' to open
discussions, If in future if the postings are overflowing, let us merge it
with appropriate ones.

[vinu]

Shall we begin the modeling?

Thus we have looked at the necessary aspects to begin the simple Numerical Ocean Model. Now a question remains, well, we have a set of equation in 1 to 3. And we mentioned the initial model velocity as zero and sea-surface is assumed to be flat so that there is no pressure gradient either. In such a situation, how does the Ocean begin to move?. Thus it is clear that somehow the motion should initiate. It is a good practice to give some particular 'sea-surface elevation' as the initial condition. This will start the motion in our Ocean Model. However there are other important 'forcing' in the atmosphere which immediately impart momentum to the Ocean model (perhaps most of the large scale ocean circulation are driven by this 'forcing', and that is the 'wind-forcing'. Thus surface wind stress can be added to the momentum equations (1) and (2).

How does these wind go into the ocean model?

According to the Navier-Stokes equation, the stress component in the momentum equation is represented as X/rho, where X is the body force. Here in the real Ocean case the body force is the wind stress and can be added to the equation 1 and 2 as Tx/rho and Ty/rho, where 'Tx and 'Ty' are wind stress in 'x' (zonal) and 'y'(meridional) direction.

'To what depth is this wind stress going to impart its momentum to the ocean?.

Well, we know that, wind stress causes the 'Ekman spiral' in the ocean applicable to an 'Ekman depth'. In our simple 1&1/2 layer version we have only one active layer which has a 'unique' velocity. So we are going to assume that the surface wind is going to affect the entire layer. This can be simply derived using the following concept. In the 'active layer' the velocity is unique. Thus if we integrate the equation 1 and 2 vertically over the depth of the layer 'H', we getH.du/dt = H.(fv) - H.dP/dx + (1/rho).(Ts-Tb), where Ts is the surface stress (of course due to the wind) andTb is the bottom stress- due to the bottom layer friction.

However we can safely neglect the term Tb, so that the above equation becomes du/dt = f.v - dP/dx + (1/rho).Ts/H. So the wind stress is ordered as (1/rho).Ts/H, where 'H' is the depth of the 'upper-layer'.

The above all points are some necessary notes for modeling the Ocean circulation. However, lets consider a situation. Suppose the ocean is forced with a wind stress for some period of time and suddenly the wind is switched off. Now the ocean will reach stability according to equation 1 to 3. However, from the common knowledge we know that after a long time the system eventually comes back to rest (in the real ocean). Why? . It is because the ocean has friction. So we have to account this amount of friction in the ocean so that it accounts the 'dissipation of the motion'. Of course the ocean water has molecular viscosity' as every other fluid. However the scale of such viscous force (frictional force) is very few and not appropriate for large scale motion. However we have the enlarged version of this friction for the ocean, and we refer it as 'Diffusion' or 'eddy diffusion'. This term is represented as Am.d^2u/dx^2., where Am. is called the 'eddy viscos
ity coefficient'. Perhaps this kind of frictional term is necessary for computational stability.

So finally we have to add two more term to the R.H.S of equation 1 and 2.

They are + Tx/(rho.H) + Am.d^2u/dx^2 . Reference

The pioneer works in ocean modeling started from 60's and 70's. Most of the fundamental works in reduced gravity layered model are initiated by several scientist. Among them, the team led by J.J. O'Brien is outstanding. He and his students made many contributions to the development of simple layered models to the state-of-the-art Ocean General Circulation models (OGCM). A number of classical papers can be obtained from here (scanned copy)
(J.J. O'Brien. (http://www.coaps.fsu.edu/bios/obrien.html)



Jensen T. G., 1991: Modeling of the Seasonal Undercurrents in the Somali currentSystem, J. Geophysical Research, Vol. 96, No. C12, 22,151-22,167.

Jensen T. G., 1996: Artificial retardation of barotropic waves in layered ocean models. Monthly Weather Review, 124, 1272-1283.

Busalacchi, A. J., and J. J. O'Brien, 1980: The seasonal variability in a model of the tropical Pacific. J. Physical Oceanography, 10, 1929-1951.

McCreary, J. P., and P. K. Kundu, 1988: Numerical investigation of Somali current during southwest monsoon. J. Marine Research, 46, 25-58.

McCreary, J. P., P. K. Kundu and R. L. Molinary, 1993: A numerical investigationof dynamics and thermodynamics and mixed layer process in Indian Ocean, Progress in Oceanography, 31, 181-244




PS: If any of the above snetence found miss-leading PLEASE FEED BACK,
Excuse is requested for type-mistakes.
-------------------------------------------------------------------------------------------------------

[praveen]

hi vinu
The above thread is becoming more attractive. I have a doubt about open boundaries in models. If we set an open boundary in a semi enclosed region, like if we model Indian ocean basin. we have to set southern part as open boundary. then how the system become stable. Do we need any aditional information in that region?? . In some model I saw inclusion of a sponge in that region. Can you explain this in your thread?

[vinu]

Hi Praveen,
Yes, open boundaries need some more explanation. I hope to add it
with the above thread as time comes.

There are standard methods to deal with 'Open boundaries', In short,
some of them can be listed as (1) Dirichlet B.C. (2) Neumann B.C. (3)
Newtonian damping (4) Sponges (5) Radiation B.C (6) Cyclic B.C etc...

In the Dirichlet, prescribed values are specified at the boundary (can be
from observational data or some other model) and 'clamped' to the
model at each time-step.

In Neumann, it assumes that, there is 'no gradient' of flux across the
boundary, i.e. the flux are freely allowed to flow out/in to the domain.

In Newtonian damping, the variable (eg. Temp or Salinity )are damped to
a reference value at the open boundaries. This is like a 'relaxation to the
pre-defined values'.

Sponges are 'like a cushion' at the boundaries. This is accomplished by
making the fluid 'high viscous' (like oil) (increasing the values of Am) at
the boundaries so that high fluctiations 'impinging' on the boundary is damped quickly.

Radiation boundaries, as the name implies, it allows the 'waves to radiate' out of the domain.

Cyclic boundary let the waves exiting from one side of the domain to
enter through the another side. For example,the waves leaving the western
boundary will enter into the model domain through the easter boundary.
This is accompanied by restoring values of the 'boundaries' to the values
of its opposite side.

well, this is a short discription, each of the above method comes with
advantage and penalty.

The first book listed here (http://www.oceanographers.net/forums/showthread.php?p=131) explains these in detail.
------------------------------------------------------------------------------------------------------

[vinu]

Detailed Physical and Numerical aspects of the Ocean and Atmospheric
Modeling is described in the following text book

Fundamentals of Ocean Climate Models, (wwwoa.ees.hokudai.ac.jp/people/vinu/share_pal/griffies_ocean_modeling_book.pdf)

by Stephen Griffies, GFDL, Princeton University (http://www.gfdl.noaa.gov/~smg/).

Here is a work book for Ocean modeling practices, (wwwoa.ees.hokudai.ac.jp/people/vinu/share_pal/workbook.pdf)

by Zygmunt Kowalik, (http://www.sfos.uaf.edu/directory/faculty/kowalik/) which comes with necessary FORTRAN subroutines.

PS: These two books are contributed to our web-site by our member Ramon Garcia (http://www.oceanographers.net/forums/member.php?u=282).
Thank you Ramon, your contribution is highly appreciated.


Thank you.
------------------------------------------------------------------------------------------------

[sreenesh]

Thanks vinu to present the article, it is very useful for new comers in the sence that it will give better idea about modeling. hope that more such interesting and useful articles will be coming from you and other members.

[rajeshkumar]

Hi Vinu
The materials you had posted is really worthful. Hope you will repeat this kind of things
Thank You
Rajeshkumar

[kkvarma]

Dear Dr.vinu,

The post appears to be very useful. What are the input formats for the code . Or is the model TOW that is the only input. Also, can this model work for coastal seas.
K.K.Varma

[vinu]

Thank you for the interest and feedbacks on this post


The above model is the simplest representation of the deep ocean, with a single active (shallow) layer over a (deep) quiescent layer. Also the motion concerned is only about a particular baroclinic mod, which is apparent in large scale oceanic adjustments to the external forcing. Still, the model is only hydrodynamic and not included the buoyancy effect (the evolution of density).

Coastal seas (shallower than 300m) are very different in processes and circulation compared to abyssal part (>4000m) of the ocean. Thus the model of the above kind WILL NOT be a suitable choice for coastal seas. The important difference is the geography (depth wise) of the coastal region compared to the deep ocean basins. The coast begins with near zero depth at the shore and a gradual steep (still shallower than 200m) within a 200km off shore, which then steeps into 4000m abyssal basin in a further few kilometer (2km to 10km). This turns out an important feature that should be included in the coastal modeling is the bottom topography. This additional necessity truncate our idealisation of an one-and-half-layer representation of the ocean.

The inclusion of bottom topography (a simple approach) can be done for a barotropic coastal model (i.e. neglecting the vertical and horizontal density gradient and assuming the entire depth is uniform density). This is not a serious limitation, because the coastal dynamics are largely barotropic in nature, responding to the alongshore wind forcing. The simplest approach of this kind can be established by considering the term 'H' in continuity equation (equation ---(3)) as a variable off shore. Thus the continuity equation modifies and additional terms appear in equation --(3), namely dH/dx or dH/dy. Assuming, the depth is varying only offshore (not along the coast) and fixing the coastal line North-South, the term dH/dy can be neglected. The remaining part (dH/dx) is indeed the bottom slope 'S'. Thus it turns out that, one important scale governing the coastal model solutions will be depending upon the bottom slope.

The same principle can be extended to a baroclinic ocean with different density layer with depth. For example, a simplest approach will be a two-layer ocean with shallow upper layer and deep bottom layer, but still the bottom layer will be under motion, because the fluid depth varies off shore. Thus the one-and-half-layer has to be modified to a two-layer system.

Also the precesses involved in the coastal dynamics are different from small scale adjustments of deep oceans. The coastal dynamics involves a spectrum of waves, those ranging from sub-inertial scale (period comparable to the inertial period of rotation of the earth) to large time scale but less compared to quasi-geostrophic motions (very slow adjustments). These are called continental shelf waves. These generate either by local wind forcing or remotely propagate into the region from elsewhere, because the coastal region acts as a wave-guide for energies to propagate. A more detailed description of coastal modeling and continental shelf waves can be found in text book listed here (item:1) (http://www.oceanographers.net/forums/showthread.php?p=131).

Now, coming back to the one-and-half-layer system described in the beginning of this post, the forcing is only the wind stress. In the real ocean, (with different density layers and bottom friction), other forces come into play. For example the bottom stress (friction) and the friction at the interface of different density layers also has to included. A fortran-code for the above model is given here (http://www.oceanographers.net/forums/showthread.php?p=1254).



Thanks
---------------------------------------------------------------------------------------------

[kkvarma]

Dear Vinu,
I tried to run the rgmodel. The out put is recorded in some unknown format. Looking at the programe, it shoud write in some format readable with notepad. .While runing the progrogramme on the screen 'writing output for day -> ' is shown as expected. Woud you kindly clarify what might have gone wrong?

Thanking you
K.K.Varma

[praveen]

output is in binary format (check lines 77-78)

Grads will be better for viewing the out put, using the *.ctl file given in post (http://www.oceanographers.net/forums/showthread.php?p=1254)

[vinu]

Output Format

yes, what Praveen said is right. The output is in binary format. The
print appear in the monitor is only to indicate that the model is running.
The output is written in 'out.dat' file.

If required, output can be produced in ASCII format (readable in
notepad), by replacing the line -77,78 and 93 of the code with,

open (1,file='out.asc')
write (1,*) ((h(i,j,3),i=1,lon),j=1,lat)

respectively.


Binary format can be read and plot by using GrADS software.
A cotrol file for GrADS is included in the original post
(as Praveen pointed out).

To find more information about the software-GrADS, see this post (http://www.oceanographers.net/forums/showthread.php?t=8). It is
freely available for both Windows and Unix platforms.



-Thankyou for the feedback on this post

-

[kkvarma]

Dear Dr.Vinu,
Thank u very much for ur promt response.
I have encounterd some problems in using out put in Grads. The out.dat file does not open. It gives the message 'can't open description file' . Of course, consequently rgmodel.ctl file does not run. Would yu kindly enlighten on this

Regards
kk varma

[praveen]

I haven't faced any problem in grads, to open this out put file.

In Grads, you should open the rgmodel.ctl instead of out.dat file

[vinu]

Ploting outputs in Grads

Let me list out the necessary steps for ploting the output.

1) I suppose, the model run, and produced the output file 'out.dat'.

2) From the same directory (or location), where the model
output (out.dat) is written, open grads.

3) from the "Grads command prompt", give the command
open rgmodel.ctl
make sure, both the files "rgmodel.ctl and out.dat" are in the same
directory or location. The rgmodel.ctl is a 'control file' where the
necessary information about the data has been coded. Grads receives
these informations and assign it to the data.
Also please refer the grads manual pages. (http://www.oceanographers.net/forums/showthread.php?t=8)

make sure that the directory (where the Grads is opened) contains
the "out.dat and rgmodel.ctl". If your defualt Grads location is
different from model directory, (most probably this would be the case
with Windows) please edit the file name in the 'rgmodel.ctl' with exact
path. For example,

DSET C:\My_work_area\Model\out.dat (i hope this is the path
specification in Windows, please correct it if wrong).

DSET /data1/My_work_area/Model/out.dat (for linux/unix)

4) on a successful opening, Grads might list out like following,

Scanning description file: rgmodel.ctl
Data file out.dat is open as file 1
LON set to 0 49
LAT set to 0 24
LEV set to 1 1
Time values set: 1900:1:1:0 1900:1:1:0

5) the command 'query file' will list out the dimensions and
variable names of the loaded data.

6) the command 'set t 100' will set the time level to 100.

7) the command 'display h' will plot the variable "h".


(refer grads manual pages (http://www.oceanographers.net/forums/showthread.php?t=8) for more information on ploting)


If none of the above steps work, please respond to this post.


-Thank you for the feedback on this post.

-

[kkvarma]

Dear Dr. Vinu,
I could plot the height. The figure shows height of about 200 near the equator and near the northern boundaries. Could this be the expected result? I am not very sure whether the model run was correct? Is there any modification to be done in the rgmodel?
regards
kkvarma

[vinu]

Solution with the given forcing.

If the model run with the given forcing provided in the code,
(a simple zonal wind stress profile represented by a sin function
with easterlies (trade winds) in the equatorial region and
westerlies on the northern domain, the solution will be a
gyre circulation, with strong and narrow western boundary
currents (parallel to the western coast)
and broad and slow eastern boundary currents
(parallel to the eastern coast) and zonal flow everywhere else.
The corresponding thermocline (or thickness)
shape will be, larger at the center of the
domain and gradient to the domain boundaries, with strongest gradient
in the western boundary (represents the strong western boundary currents).

An example plot is given below.

(Optional)
Also, while plotting, the one grid border at the
domain boundaries can be omitted,
because that is used for applying boundary condition and is treated as
land cell. This border grid can be avoided in the plot by limiting
the output to a dimension (lon-2)*(lat-2). Following alterations
in the code may help to avoid this land cell from being printed in the output.

(with reference to the code given in the post).

line-78 & recl=(lon-2)*(lat-2)*4)
line-93 write(1,rec=lrec)((h(i,j,3),i=2,lon-1),j=2,lat-1)

Also, in the corresponding rgmodel.ctl file, the dimension values
50 and 25 has to replace with 48 and 23, respectively.


http://wwwoa.ees.hokudai.ac.jp/people/vinu/share_pal/model/solution.gif


-
Thanks

[nonicad]

Do have code fortran prog fom the book :Zygmuny kowalik,

[vinu]

Sorry, I dont have the code you requested,

-Vinu

[cd264]

Dear vinu:
quite useful post for me!
long for your new lectrue.:thumbsup: