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://www.oceanographers.net/upload/vinu/layermodel/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://www.oceanographers.net/upload/vinu/layermodel/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://www.oceanographers.net/upload/vinu/layermodel/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://www.oceanographers.net/upload/vinu/layermodel/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://www.oceanographers.net/upload/vinu/layermodel/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:

[cd264]

hi! vinu
Is the model can be apply for the globe circulation simulation?how large can the model to simulate?what about the boundary conditions and how to prevent the unstable?

Also,happy newyear！

[vinu]

Hi Tom,

The model we were discussing in earliers posts is only a simplified
representation of thermocline without any consideration on
non-linearity and other surface ocean processes. Moreover, the
application of betaplane is limitted to only a few degrees because the
beta-value has variation with latitude (cosine function of latitude).

The above model doesnot represent the mixed layer, and thermodynamics.
Additing those two, along with non-linearity becomes a powerful tool for
tropical or midlatitude circulation. Numerous literature exists describing
the layered model with dynamics and thermodynamics and used for
several application. McCreary et al. (1994), various papers of
Jensen (1990,94,98,2003), Murtugudde et al. (1998), Qu et al. (1994,..)
etc are some examples.

The model described here, thus, is not a suitable choice for any global
circulation study. But indeed, this model could resolve some of the
fundamental structure of ocean dynamcis, such as wind driven
gyre circulation for a limitted extend of latitudilan width for which the
beta plane approximation is applicable with a given value of beta or
equatorial-mid latutude wave dynamics etc.. I have noted it in the
orignal post (second paragraph). http://www.oceanographers.net/forums/showthread.php?p=1254

Literature exists using layer model with only dynamics and used for particular
application such as equatorial dynamics. Paper by Potemra et al. (2001), JGR,
Vol,106, 2407-2422 used layer model to study the propagation of
equatorlial rossby waves in the pacific and its contribution to Indian Ocean
Rossby waves.

-Thanks

[Malte]

Dear Forum,
i have tried to translate the rgmodel fortran code into a matlab file, but even though i have choosen a real small timestep compared to delta_x, the model keeps 'blowing' up, if I use a small disturbance of h as initial condition.

I tried really long finding the error in my code, but haven't been successful yet. I attach it to this post, maybe someone sees quick what goes wrong, without spending too much time on it. I am sitting on it for a few days now and still no insight...
Thank you very much in advance!
Cheers,
Malte

[Malte]

Maybe I'll post the numerical schemes here, so one doesn't have to download the file.

The integration scheme (centered in time and space):


for j=2:jmax-1
for i=2:imax-1

% calculate the coriolis force at the same position, one step in the past (km):
coriolis = f(j)*v(i,j,km);

% calculate the horizontal pressure gradient
% (finite diff. centered in y space of the first derivative)::
horiz_pressure = -g_reduced*(h(i+1,j,k) - h(i-1,j,k))/d2x;

% Diffusion term: A_H * Laplace*U
% Laplace*U = (d^2 U / d x^2) + (d^2 U / d y^2)
% Centered finite difference of the second derivative
diffx = (u(i-1,j,k) - 2*u(i,j,k) + u(i+1,j,k))/dx2;
diffy = (u(i-1,j,k) - 2*u(i,j,k) + u(i+1,j,k))/dy2;
diffusion = A_H*(diffx + diffy);

% windstress
windstress = 0;

% calculate U:
u(i,j,kp)=u(i,j,km)+delta_t*(coriolis-horiz_pressure+diffusion+windstress);
end
end



% Calculate meridional transport (y-direction)
% Finite Difference Formulation of dV/dt:
for j=2:jmax-1
for i=2:imax-1

% calculate the coriolis force at the same position, one step in the past (km):
coriolis = -f(j)*u(i,j,km);

% calculate the horizontal pressure gradient
% (finite diff. centered in y space of the first derivative):
horiz_pressure = -g_reduced*(h(i,j+1,k) - h(i,j-1,k))/d2y;


% Diffusion term: A_H * Laplace*V
% Laplace*V = (d^2 V / d x^2) + (d^2 V / d y^2)
% Centered finite difference of the second derivative
diffx = (v(i-1,j,k) - 2*v(i,j,k) + v(i+1,j,k))/dx2;
diffy = (v(i-1,j,k) - 2*v(i,j,k) + v(i+1,j,k))/dx2;
diffusion = A_H*(diffx + diffy);

% windstress
windstress = 0;

% calculate V:
v(i,j,kp)=v(i,j,km)+delta_t*(coriolis+horiz_pressu re+diffusion+windstress);
end
end




% Continuity equation.
% Using forward finite difference instead of centered finite difference!
for j=2:jmax-1
for i=2:imax-1

% hdiv = dU/dx + dV/dy
% use centered finite difference:
% hdiv = (-H*(u(i+1,j,k) - u(i-1,j,k))/d2x) + ((v(i,j+1,k) - v(i,j-1,k))/d2y);
hdiv = (-H*((u(i+1,j,k)-u(i-1,j,k))/d2x+(v(i,j+1,k)-v(i,j-1,k))/d2y)+ A_H*((h(i+1,j,k)+h(i-1,j,k)-2.0*h(i,j,k))/dx2 + (h(i,j+1,k)+h(i,j-1,k)-2.0*h(i,j,k))/dy2));
% new elevation of the thermocline h:
h(i,j,kp) = h(i,j,km) + delta_t*hdiv;


end
end



with the boundary cond. applied every timestep:




for j=1:jmax

% West:
u(1,j,kp) = -u(3,j,kp); % no normal flow across the rigid boundary
u(2,j,kp) = 0.0; % no-slip at the boundary

v(1,j,kp) = -v(3,j,kp); % same as above but for v-velocity
v(2,j,kp) = 0.0; %

h(1,j,kp) = h(3,j,kp); % thickness restored to neighbouring value

% East
u(imax,j,kp) = -u(imax-2,j,kp);
u(imax-1,j,kp) = 0.0;

v(imax,j,kp) = -v(imax-2,j,kp);
v(imax-1,j,kp) = 0.0;

h(imax,j,kp) = h(imax-2,j,kp);
end

for i=1:imax

% North:
u(i,1,kp) = -u(i,3,kp);
u(i,2,kp) = 0.0;

v(i,1,kp) = -v(3,j,kp);
v(i,2,kp) = 0.0;

h(i,j,kp) = h(i,3,kp);


% South:
u(i,jmax,kp) = -u(i,jmax-2,kp);
u(i,jmax-1,kp) = 0.0;

v(i,jmax,kp) = -v(i,jmax-2,kp);
v(i,jmax-1,kp) = 0.0;

h(i,jmax,kp) = h(i,jmax-2,kp);

end

[vinu]

HI Malte,

I have seen one of your (other post?) asking about C.F.L. limits.
I think it is this post only.

I havent checked your code in detail, but a quick look found that

diffx = (u(i-1,j,k) - 2*u(i,j,k) + u(i+1,j,k))/dx2;
diffy = (u(i-1,j,k) - 2*u(i,j,k) + u(i+1,j,k))/dy2;

your diffx and diffy are doing same operation. The first should be derivative
in x and second should be in Y. that means. u(i,j-1,k) - 2*u + u(i,j+1,k)

If I could find any other fixes I will report,

Also please see the "Numerical ocean model" post for C.F.L. details.




-Vinu

[Malte]

Hi Vinu,

thank you very much! That was indeed the problem. Finding mistakes like that in code can be quite frustrating.
Thanks, I already found your post about the CFL shortly after made my post and edited it afterwards.
When I keep running the m file i have written in a simmilar domain like you used (central latitude:30, total lat. of 30, total lon. of 40)
the model shows results like I would expect (with a zonal windfield forcing). The velocity vectors u and v follow the windfield pattern and a clockwise gyre develops:
http://g.imagehost.org/t/0149/a1_9.jpg (http://g.imagehost.org/view/0149/a1_9)
http://g.imagehost.org/t/0274/a2.jpg (http://g.imagehost.org/view/0274/a2)


Then I tried to use a domain that would more or less represent the north atlantic (central latitude: 30, total lat. of 60, total lon. of 80).
But the result is only a flow in x-direction on the bottom of the domain following the windfield. At the top there is no flow following the Westerlies:
http://f.imagehost.org/t/0771/b1.jpg (http://f.imagehost.org/view/0771/b1)
http://g.imagehost.org/t/0976/b2.jpg (http://g.imagehost.org/view/0976/b2)

I have no clue why this happens. I think i got the coriolis force right for the domain and also the wind pattern.
I only changed the domain specifications.

With kind regards,
Malte

(Attached is the corrected code if anyone is interested in it)







, it keeps running normal

[vinu]

Hi Malte,


But when you attempt to stretch the code for a domain with meridional
width of 60 degrees, putting a central latitude at 30, you are violating
the beta-plane assumption. The beta-plane assumption to hold, the
latitudinal variation of Coriolis force should be linear on either side of
the central latitude. This satisfies (approximately), when you limits your
'wings' into a domain of 15 degrees on either sides of central latitude.
Beyond this arbitrary limits, the beta value itself is changing, so your
condition of f=f0+beta.dy will be nonlinear. So, we hardly expect a single
Gyre dynamics in this geometry.

If you are into application of this model for a particular research purpose,
(1) either you should choose a domain satisfying beta-plane assumption
(2) or you should convert the model Cartesian coordinates into spherical coordinate,
thereby, taking Coriolis force f=2.omega.sin(theta) directly. There you
get f-plane, beta-plane and changing beta value all in one geometry.

In doing so, you will come across with the problems of converging merdions
at high latitude because of variable zonal resolution as you move northward.
You should choose a smaller time-step for integration to satisfy CFL in
converging meridian, or improving your laplacian diffusion suitably by adjusting
a diffusion coefficient as a function of latitude.

The question you have raised (why is the present model unable to give a
successful gyre with a stretched geometry) has to be answered with
the backup of gyre dynamics. In the Rossby-adjustment problems in
spherical co-ordinate, the time-scale for wave propagations are related
to the local beta-value. In the equatorial area, the Rossby radius of deformation
is larger and time-scale is smaller, while in the higher latitude the
time scale is larger for a wave to propagate from east to west. In this case,
the adjustment of ocean in-response to the wind forcing will have
gyres particularly related to the domain geometry and forcing.
Interesting solutions in this case are studied in classical papers of
ocean spin-up problems by several oceanography pioneers
starting from Sverdrup, Munk and so on...

Still, if you keep on spinning up your ocean in the present geometry
you will find your circulations changing, because, I think the domain
will take a longer time to equilibrate.

Zheng Shaobin ( summerpca_at_yahoo.com ) also followed a similar
track of creating an rgmodel in Matlab, later moved to a model with
spherical goemtry and non-linearity in dynamics. I think he was into an
application of the reduced gravity model for a particular research purpose,
as of the latest information, he found a spherical coordinate, fully
non-linear, linux & fortran based reduced gravity ocean model
from Texas-Am university.


PS:

Hello Zheng, do you have any updates to share in this forum about your
recent advancement in using your reduced gravity model? We havent
communicated recently after our previous discussions.

-Vinu

[Malte]

Hi Vinu,

thank you for the explanation! I haven't thought about that.
I have found another application of the model at:
http://ic.ucsc.edu/~ammoore/ocea290g/MATLAB/ (http://ic.ucsc.edu/%7Eammoore/ocea290g/MATLAB/)
(the file is called 'shallow water template'). The PDE's are just slightly different from your model. They use a 60 degree latitude domain size, but centered at the equator. So the variation of the coriolis parameter is just changing by 30 degrees in each hemisphere. So i guess that is the reason the beta plane linearization is justified in this case?

Cheers,
Malte

[vinu]

Hi Malte,

That is quite right. The beta-plane models are suitable for tropical
oceans, same as the case shown just above. Choosing a beta of equator
and a +/- 30 wings of the tropics is a normal practice.
The f-plane approximation is quite invalid in equatorial ocean
as it becomes zero near the equator. But, the beta, which is gradient of f
in y attains maximum amplitude (because of cos(theta)) at the equator.

Although, the approximation still has error but only a small percentage.
Considering the sphericity of the earth, half of the earth lies at latitudes of less
than 30-degrees and a maximum percentage error in the beta-plane
approximation in that range of latitudes is only 14%. (Gill, 1982, Chapter-11.4,
Tropics).

But putting a wing of +/-30 degrees from the mid-latitude ocean
(i.e. a central latitude = 30) the error in beta-plane will be larger because
variation of beta with latitude is larger as we move poleward.

Note:
The primitive equations in spherical coordinate will have geomtric terms
in additions to the terms explained in the beginning of this post. There, the
derivates are found with respect to phi and theta (i.e. lon and lat angle).

The details of primitive equations of reduced gravity models in spherical coordinates
are given in Jensen (1991), JGR. Vol. 96, NO.C12, 22151-22167.

-Vinu

[Malte]

Dear Forum,

I tried to force a spatial NPZ model (based on franks (1986)) with an advection-diffusion tracer equation, using centered space and centered time (leapfrog) finite differences. The numerical routine is based in great parts on the one published here for the rgmodel. Somehow I only get reasonable results if I start with the same initial conditions (concentrations for N,P and Z) at every grid point. As soon as I use different initial concentrations, the routine fails in my program. I am just starting to learn Fortran, therefore the problem might be based on my lack of understanding of some basic concepts.

Hope someone might be able to help me with the problem.
Thanks!

The main program npz.f90:



PROGRAM NPZ


USE solver

IMPLICIT NONE

!--------------------------------------------!
! DECLARATIONS
!--------------------------------------------!

REAL :: N,P,Z,Mass,Masstotal,u,v,t,d2x,d2y,dx2,dy2,day2sec ,d2t
INTEGER :: i,j,tsteps,loop,l,ltemp,lminus,lplus

!--------------------------------------------!
! GLOBAL variables and par.
!--------------------------------------------!

INTEGER, PARAMETER :: imax = 20
INTEGER, PARAMETER :: jmax = 20
COMMON /var/ N(imax,jmax,3), P(imax,jmax,3), Z(imax,jmax,3), Mass(imax,jmax), u(imax,jmax), v(imax,jmax),t
COMMON /par/ d2x,d2y,dx2,dy2,day2sec,d2t,l,ltemp,lminus,lplus

!--------------------------------------------!
! INPUT NAMELIST
!--------------------------------------------!

REAL :: alpha, vmaxNP, vmaxPZ, Ks, mP, mZ, Phalf, Pzero, dt, Kx, Ky
INTEGER :: tstart,tend,dx,dy,check

NAMELIST /NPZNAME/ alpha, vmaxNP, vmaxPZ, Ks, mP, mZ, Phalf, Pzero, tstart, tend, dt, dx, dy, Kx, Ky, check

OPEN(UNIT=41,FILE='npz.nml',STATUS='OLD')
READ(41,NML=NPZNAME)
CLOSE(41)



!--------------------------------------------!
! HELP VARIABLES
!--------------------------------------------!

d2x = 2.0*dx
d2y = 2.0*dy
dx2 = dx*dx
dy2 = dy*dy
day2sec = 86400.0 ! seconds per day




!--------------------------------------------!
! OUTPUT
!--------------------------------------------!

OPEN(UNIT=11,FILE='ijtnpz.out',STATUS='UNKNOWN')
OPEN(UNIT=12,FILE='tmass.out',STATUS='UNKNOWN')

!--------------------------------------------!
! INITIAL CONDITIONS
!--------------------------------------------!

!N = 7.32 ! (mmoles m^-3)
!P = 1.01 ! (mmoles m^-3)
!Z = 0.23 ! (mmoles m^-3)

! Allocate space for the field variables

! Initial Cond. at inner points for state variables

DO i = 2,imax-1
DO j = 2,jmax-1
N(i,j,1) = 0.0
P(i,j,1) = 0.0
Z(i,j,1) = 0.0
END DO
END DO


DO i = 5,8
DO j = 5,8
N(i,j,1) = 7.32
P(i,j,1) = 1.01
Z(i,j,1) = 0.23
END DO
END DO



! Initial Cond. at outer points (Boundary)
DO i=1,imax
N(i,1,1) = 0.0 ! at the west border
P(i,1,1) = 0.0
Z(i,1,1) = 0.0

N(i,jmax,1) = 0.0 ! at the east border
P(i,jmax,1) = 0.0
Z(i,jmax,1) = 0.0
END DO


DO j=1,jmax
N(1,j,1) = 0.0 ! at the north border
P(1,j,1) = 0.0
Z(1,j,1) = 0.0

N(imax,j,1) = 0.0 ! at the south border
P(imax,j,1) = 0.0
Z(imax,j,1) = 0.0
END DO

! Set all future states zero

DO i = 1,imax
DO j = 1,jmax
N(i,j,2) = 0.0
P(i,j,2) = 0.0
Z(i,j,2) = 0.0
N(i,j,3) = 0.0
P(i,j,3) = 0.0
Z(i,j,3) = 0.0
END DO
END DO





t = 0 ! Initial time


DO i=2,imax-1
DO j=2,jmax-1
! Write initial conditions in file:
WRITE(11,'(2I5,4E12.4)')i,j,t,N(i,j,1),P(i,j,1),Z( i,j,1)
END DO
END DO



! Initialize time index
d2t = dt ! at the begining forward intergration for first time step
lplus = 3 ! future time index
l = 2 ! present time index
lminus = 1 ! past time index


!--------------------------------------------!
! MAIN PROGRAM
!--------------------------------------------!


! Time Loop
DO loop = 1,2000
t = loop*dt/day2sec ! actual time (days)



!--------------------------------------------!
! Call Solver !
!--------------------------------------------!

CALL leapfrog(loop)




END DO



ENDFILE(UNIT=11) ! end output file
CLOSE(UNIT=11) ! close output fie
ENDFILE(UNIT=12) ! end output file
CLOSE(UNIT=12) ! close output fie


1000 END PROGRAM NPZ


integration routine:



MODULE solver
CONTAINS

SUBROUTINE leapfrog(loop)


!--------------------------------------------!
! GLOBAL variables and par.
!--------------------------------------------!

INTEGER, PARAMETER :: imax = 20
INTEGER, PARAMETER :: jmax = 20

COMMON /var/ N(imax,jmax,3), P(imax,jmax,3), Z(imax,jmax,3), Mass(imax,jmax), u(imax,jmax), v(imax,jmax), t
COMMON /par/ d2x,d2y,dx2,dy2,day2sec,d2t,l,ltemp,lminus,lplus


!--------------------------------------------!
! LOCAL DECLARATIONS
!--------------------------------------------!

REAL :: NPuptake(imax,jmax), PNmort(imax,jmax), ZNmort(imax,jmax)
REAL :: PZgraz(imax,jmax), dNdt(imax,jmax), dPdt(imax,jmax), dZdt(imax,jmax)


!--------------------------------------------!
! GLOBAL DECLARATIONS
!--------------------------------------------!

REAL :: N,P,Z,Mass,Masstotal,u,v,t,d2x,d2y,dx2,dy2,day2sec
INTEGER :: i,j,tsteps,loop

!--------------------------------------------!
! READ NAMELIST
!--------------------------------------------!


REAL :: alpha, vmaxNP, vmaxPZ, Ks, mP, mZ, Phalf, Pzero, dt, Kx, Ky
INTEGER :: tstart,tend,dx,dy,check

NAMELIST /NPZNAME/ alpha, vmaxNP, vmaxPZ, Ks, mP, mZ, Phalf, Pzero, tstart, tend, dt, dx, dy, Kx, Ky, check
OPEN(UNIT=42,FILE='npz.nml',STATUS='OLD')
READ(42,NML=NPZNAME)
CLOSE(42)



!--------------------------------------------!
! READ velocity vectors (u,v) !
!--------------------------------------------!

! assign constant u and v:
DO i = 1,imax
DO j = 1,jmax
u(i,j) = 0.01 ! 0.0 0.1
v(i,j) = 0.01 ! 0.0 0.1
END DO
END DO



!--------------------------------------------!
! Apply Boundary Cond. at each timestep !
!--------------------------------------------!

! Cyclic Boundary conditions

DO i=1,imax
N(i,1,l) = N(i,jmax-1,l) ! at the west border
P(i,1,l) = P(i,jmax-1,l)
Z(i,1,l) = Z(i,jmax-1,l)

N(i,jmax,l) = N(i,2,l) ! at the east border
P(i,jmax,l) = P(i,2,l)
Z(i,jmax,l) = Z(i,2,l)
END DO


DO j=1,jmax
N(1,j,l) = N(imax-1,j,l) ! at the north border
P(1,j,l) = P(imax-1,j,l)
Z(1,j,l) = Z(imax-1,j,l)

N(imax,j,l) = N(2,j,l) ! at the south border
P(imax,j,l) = P(2,j,l)
Z(imax,j,l) = Z(2,j,l)

END DO


!--------------------------------------------!
! Convert from daily rates to rates per second
!--------------------------------------------!


vmaxNP = vmaxNP/day2sec
vmaxPZ = vmaxPZ/day2sec
mP = mP/day2sec
mZ = mZ/day2sec




! Mass blance at each time step
Mass = 0.0
Masstotal = 0.0


!--------------------------------------------!
! Loop over the grid !
!--------------------------------------------!



DO i = 2,imax-1
iplus = i+1
iminus = i-1
DO j = 2,jmax-1
jplus = j+1
jminus = j-1



!--------------------------------------------!
! LOCAL BIOLOGY PROCESSES (at current time)
!--------------------------------------------!

NPuptake(i,j) = vmaxNP*(N(i,j,l)/(Ks + N(i,j,l)))*P(i,j,l)

! linear mortality:
PNmort(i,j) = mP*P(i,j,l)
ZNmort(i,j) = mZ*Z(i,j,l)

! Ivlev Saturation:
PZgraz(i,j) = vmaxPZ*(1 - exp(-1.0*P(i,j,l)))*Z(i,j,l)




!--------------------------------------------!
! LOCAL BIOLOGY SOURCE TERMS (S)
!--------------------------------------------!


dNdt(i,j) = -NPuptake(i,j) +PNmort(i,j) +ZNmort(i,j) +(1-alpha)*PZgraz(i,j)
dPdt(i,j) = NPuptake(i,j) -PNmort(i,j) -PZgraz(i,j)
dZdt(i,j) = alpha*PZgraz(i,j) -ZNmort(i,j)

!--------------------------------------------!
! Centered Time and Space
!--------------------------------------------!


N(i,j,lplus) = N(i,j,lminus) + ( &
-u(i,j)*((N(iplus,j,l) - N(iminus,j,l))/d2x) & ! Advection in x-direction
-v(i,j)*((N(i,jplus,l) - N(i,jminus,l))/d2y) & ! Advection in y-direction
+Kx*((N(iplus,j,l) + N(iminus,j,l) - 2.0*N(i,j,l))/dx2) & ! Diffusion in x-direction
+Ky*((N(i,jplus,l) + N(i,jminus,l) - 2.0*N(i,j,l))/dy2) & ! Diffusion in y-direction
+dNdt(i,j) & ! Local change due to Biology
)*d2t


P(i,j,lplus) = P(i,j,lminus) + ( &
-u(i,j)*((P(iplus,j,l) - P(iminus,j,l))/d2x) & ! Advection in x-direction
-v(i,j)*((P(i,jplus,l) - P(i,jminus,l))/d2y) & ! Advection in y-direction
+Kx*((P(iplus,j,l) + P(iminus,j,l) - 2.0*P(i,j,l))/dx2) & ! Diffusion in x-direction
+Ky*((P(i,jplus,l) + P(i,jminus,l) - 2.0*P(i,j,l))/dy2) & ! Diffusion in y-direction
+dPdt(i,j) & ! Local change due to Biology
)*d2t


Z(i,j,lplus) = Z(i,j,lminus) + ( &
-u(i,j)*((Z(iplus,j,l) - Z(iminus,j,l))/d2x) & ! Advection in x-direction
-v(i,j)*((Z(i,jplus,l) - Z(i,jminus,l))/d2y) & ! Advection in y-direction
+Kx*((Z(iplus,j,l) + Z(iminus,j,l) - 2.0*Z(i,j,l))/dx2) & ! Diffusion in x-direction
+Ky*((Z(i,jplus,l) + Z(i,jminus,l) - 2.0*Z(i,j,l))/dy2) & ! Diffusion in y-direction
+dZdt(i,j) & ! Local change due to Biology
)*d2t


! Mass Balance
Mass(i,j) = N(i,j,lplus) + P(i,j,lplus) + Z(i,j,lplus)
Masstotal = Masstotal + Mass(i,j)


! Write output at each timestep for each spatial point
WRITE(11,'(2I5,4E12.4)')i,j,t,N(i,j,lplus),P(i,j,l plus),Z(i,j,lplus)





END DO ! End i loop
END DO ! End j loop



WRITE(*,*)t,Masstotal
WRITE(12,'(2E12.4)')t,Masstotal




IF(check==1) THEN

DO i = 1,imax
DO j = 1,jmax

!--------------------------------------------!
! CHECK if all state variables are non-negative
!--------------------------------------------!

IF (N(i,j,lplus) < 0) THEN
PRINT *, 'WARNING, negative state variable N!'
END IF

IF (P(i,j,lplus) < 0) THEN
PRINT *, 'WARNING, negative state variable P!'
END IF

IF (Z(i,j,lplus) < 0) THEN
PRINT *, 'WARNING, negative state variable Z!'
END IF



END DO ! End i loop
END DO ! End j loop

END IF



! apply asselin-robert filter

DO i=2,imax-1
DO j=2,jmax-1
N(i,j,l) = N(i,j,l) + 0.1*( &
N(i,j,lplus) -2.0*N(i,j,l) + N(i,j,lminus))

P(i,j,l) = P(i,j,l) + 0.1*( &
P(i,j,lplus) -2.0*P(i,j,l) + P(i,j,lminus))

Z(i,j,l) = Z(i,j,l) + 0.1*( &
Z(i,j,lplus) -2.0*Z(i,j,l) + Z(i,j,lminus))
END DO
END DO



! use leap-frog scheme from second time step and rotate the time index
d2t = 2.0*dt
ltemp = lminus
lminus = l
l = lplus
lplus = ltemp

WRITE(*,*)lminus,l,lplus



RETURN
END SUBROUTINE leapfrog

END MODULE solver



Namelist:



&NPZNAME
check = 1, ! if check==1, checks if all state variables are non-negative
dt = 500, ! (s)
dx = 2000, ! (m)
dy = 2000, ! (m)
Kx = 0.01, ! 0.0 5.0e3
Ky = 0.01, ! 0.0 5.0e3
alpha = 0.7, ! ();
vmaxNP = 2., ! (day^-1)
vmaxPZ = 1.5, ! (day^-1)
Ks = 1., ! (mmoles m^-3)
mP = 0.1, ! (day^-1)
mZ = 0.2, ! (day^-1)
Phalf = 0.5, ! (mmoles m^-3)
Pzero = 0.3 ! (mmoles m^-3)
/

[Malte]

These are the dynamics that should appear (result when initial cond. are the same on all grid points):

[vinu]

!--------------------------------------------!
! READ velocity vectors (u,v) !
!--------------------------------------------!

! assign constant u and v:
DO i = 1,imax
DO j = 1,jmax
u(i,j) = 0.01 ! 0.0 0.1
v(i,j) = 0.01 ! 0.0 0.1
END DO
END DO


Hi Malte,


I am just copying above set of codes, where you assumes an offline
currents (prescribed flow) a priori, and then using this currents and
trying to advect your N.P,Z tracers. Also the flow is two dimentional.

But the question is, what boundary condition does this flow satisfy?
If you assume u=0.1 and v=0.1, for a squar domain, the flow is
diogonal. But at boundaries, you prescribed cyclic condition. Let us
imagin, the domain is cyclic zonally, that means you can wrap the
domain by connecting east and west edges, so you end up in in
a cylinder shape. But what happens to the north and south boundaries.
They wont join. That means, applying cyclic on four boundaries wont
make the flow really cyclic unless you respect the geometry as a sphere.

Another trouble I could see is the incorrect representation of western
and eastern edges. In your domain u(i,j,1), which is the zonal dimension?.
I assume 'i' is the zonal dimension. In that case, your western boundary
is u(1,j,1) and eastern boundary is u(imax,j,1). Along the western boundary,
the 'j' varies from 1 to jmax.

If you still want to keep u, v as 0.1 and constant througout the
domain, my recommendation would be to use a radiating boundary
rather than a cyclic boundary. That means, let the flow leave the
boundary as it reaches the edge and let new flow to come in from
the entry points. In this case, your tracers wont be conserved in
the domain. To achieve radiation, simply way is to restore the
edges with its immediate neighbours.

Hope some of this might help you,

-Vinu

[Malte]

Hi Vinu,

I chose the cyclic boundary conditions mainly as a 'debugging' help, to see if I introduced some unintentional sources or sinks due to faulty programming in my scheme. Figuring out the appropriate boundary conditions caused me a lot of headaches. My final goal is to develop a numerical scheme that is appropriate to describe the dynamics of a big lake or nearly closed estuary.
The 2d code was more of a training exercise of trying to code different numerical schemes. In the end, I want to use a 3d scheme which is able to describe upwelling and downwelling at the closed boundaries. However, I am still not sure how to treat the bc for the transport equation in that case.
For example if the domain (x and y direction) looks like this (x = boundary, o = inner grid points):

x x x x x x x x
x o o o o o o x
x o o o o o o x
x i i+1 o o o o x
x o o o o o o x
x o o o o o o x
x x x x x x x x

When I calculate the tracer concentration at the point i based on centered differences (i-1 and i+1) and the concentration at x(i-1) is just set to zero, this will certainly yield to a wrong concentration estimate for the inner point.



Another trouble I could see is the incorrect representation of western
and eastern edges. In your domain u(i,j,1), which is the zonal dimension?.
I assume 'i' is the zonal dimension. In that case, your western boundary
is u(1,j,1) and eastern boundary is u(imax,j,1). Along the western boundary,
the 'j' varies from 1 to jmax.


I wanted to use the 'normal' convention first. The zonal index i for west-east direction and the meridional index j for north-south direction, which would would give as you said:

western boundary: (1,j)
eastern boundary: (imax,j)

However, doesn't treat Fortran the matrix indices the opposite way?:

(1,1) (1,2) (1,3)
(2,1) (2,2) (2,3)
(3,1) (3,2) (3,3)

Thank you very much (as always) for your help!

[vinu]

I wanted to use the 'normal' convention first. The zonal index i for west-east direction and the meridional index j for north-south direction, which would would give as you said:

western boundary: (1,j)
eastern boundary: (imax,j)

However, doesn't treat Fortran the matrix indices the opposite way?:

(1,1) (1,2) (1,3)
(2,1) (2,2) (2,3)
(3,1) (3,2) (3,3)

Thank you very much (as always) for your help!

Hi Malte,

The matrix appear in the above way in Fortran only when you print it. Because,
in Fortran first index vary fastest. To understand this, better to
try following simple code. In the following code, we are sure, our longitude are
10, and latitudes are 4. I set a 'sin' cycle from west to east. Now, apply
some boundary condition and see the values appear in screen. This
excersice may help you to know the structure of fortran arrays.

cccccccccccccccccccccccccccccccccccccccccccccccccc c
c test program for boundary conditions.

parameter (lon=10, lat=4)
real u(lon, lat)

pi = 3.14
omega = 2.0*pi/(lon*1.0)

do i = 1,lon ! eas-west
do j = 1,lat ! north-south

u(i,j) = sin(i*omega)

enddo
enddo

c
c apply boundary conditions
c

do j=1,lat
u(1,j) = -99.999 ! western boundary
enddo
do j=1,lat
u(lon,j) = -88.888 ! eastern boundary
enddo

c
c print the result
c

do i=1,lon
write(*,10) (u(i,j),j=1,lat)
enddo
10 format (4(1x,f7.3))
end

cccccccccccccccccccccccccccccccccccccc

[Malte]

thanks, I wasn't aware of that convention.

I used the following code to print it on screen in the accustomed way:


uu = transpose(u)

!
! print the result
!

do j=1,lat
write(*,20) (uu(j,i),i=1,lon)
enddo
20 format (10(1x,f7.3))


Do you have any suggestions regarding the boundary conditions problem in the tracer equation, to assure mass conservative behavior in a closed basin?

Thanks.

[Malte]

Just wanted to let you know that I used reflective boundary conditions and the model works now.

e.g. northern boundary:

do i=1,lon
T(i,lat) = T(i,lat-1)
enddo