Dynlib: Difference between revisions

From gfi
mNo edit summary
Line 104: Line 104:
Calculates the partial y derivative of dat, using centred differences. For a non-EW-cyclic grid, 0 is returned on all edges of the x,y domain. For an EW-cyclic grid, 0 is returned on the first and last latitudes.
Calculates the partial y derivative of dat, using centred differences. For a non-EW-cyclic grid, 0 is returned on all edges of the x,y domain. For an EW-cyclic grid, 0 is returned on the first and last latitudes.


=== ''grad'' : gradient of a scalar ===
=== ''grad'' : 2-D gradient of a scalar ===


<code>(resx,resy)=grad(dat,dx,dy)</code>
<code>(resx,resy)=grad(dat,dx,dy)</code>


Calculates the gradient of dat, using centred differences. For a non-EW-cyclic grid, 0 is returned on all edges of the x,y domain. For an EW-cyclic grid, 0 is returned on the first and last latitudes.
Calculates the 2-D gradient of dat, using centred differences in x and y. For a non-EW-cyclic grid, 0 is returned on all edges of the x,y domain. For an EW-cyclic grid, 0 is returned on the first and last latitudes.


=== ''lap2'' : 2-D laplacian of a scalar ===
=== ''lap2'' : 2-D laplacian of a scalar ===

Revision as of 13:55, 21 January 2013

Documentation

The steps necessary to obtain dynlib are described below. A more thorough documentation is compiled in the main documentation page.

Obtaining dynlib

  1. Copying the source code repository
    git clone /Data/gfi/users/tsp065/lib/dynlib.git
  2. Change into the dynlib folder
    cd dynlib
  3. Compile the library
    ./compile

Quick start to developing with dynlib

Editing the Fortran code

The fortran code lives in the main source code directory. At the moment there are six source code files

$ ls *.f95
dynlib_config.f95 dynlib_const.f95 dynlib_conv.f95 dynlib_diag.f95 dynlib_kind.f95 dynlib_stat.f95

The most important are dynlib_diag.f95 which contains subroutines that calculate various diagnostics, and dynlib_stat.f95 which contains statistical functions. Changed Fortran sources need to be recompiled, again using

./compile

Version control

The changes you made to the source code files can be listed by

git status

or viewed in detailed diff-comparisons by

git diff

or for one file only

git diff [filename]

Commit your changes from time to time and give a sensible and brief description of your changes in the editor that is opened (automatically)

git commit -a

The commit is then stored in your copy of the source code repository, but not yet available for others, which allows you to also commit work-in-progress.

A more thorough introduction to the version control system is given here or on the official documentation.

Using the Fortran functions

An example python script which calculates deformation using the Fortran function is provided with deformation.py.

Dynlib functions

The functions generally operate on real arrays with dimension (nz,ny,nx) where nz is number of times or levels, and ny and nx are the number of latitudes and longitudes, respectively. The function descriptions below contain detailed descriptions of arguments and returns where there is any deviation from this pattern; otherwise they may be assumed to be of the form:

Arguments:
Type Dim Description
u real (nz,ny,nx) Zonal velocity
v real (nz,ny,nx) Meridional velocity
pv real (nz,ny,nx) Potential vorticity
Returns:
Type Dim Description
res real (nz,ny,nx) output data

The ubiquitous inputs dx and dy are all of the form

Type Dim Description
dx real (ny,nx) dx(j,i)=x(j,i+1)-x(j,i-1) (in metres)
dy real (ny,nx) dy(j,i)=y(j+1,i)-y(j-1,i) (in metres)

Typically, the results for each level or time are computed individually in 2-D fashion, though they are returned as a 3-D array or the same size as the input.

ddx : partial x derivative

res=ddx(dat,dx,dy)

Calculates the partial x derivative of dat, using centred differences. For a non-EW-cyclic grid, 0 is returned on the edges of the x domain.

ddy : partial y derivative

res=ddy(dat,dx,dy)

Calculates the partial y derivative of dat, using centred differences. For a non-EW-cyclic grid, 0 is returned on all edges of the x,y domain. For an EW-cyclic grid, 0 is returned on the first and last latitudes.

grad : 2-D gradient of a scalar

(resx,resy)=grad(dat,dx,dy)

Calculates the 2-D gradient of dat, using centred differences in x and y. For a non-EW-cyclic grid, 0 is returned on all edges of the x,y domain. For an EW-cyclic grid, 0 is returned on the first and last latitudes.

lap2 : 2-D laplacian of a scalar

res=lap2(dat,dx,dy)

Calculates the 2-D laplacian of dat, using centred differences. For a non-EW-cyclic grid, 0 is returned on all edges of the x,y domain. For an EW-cyclic grid, 0 is returned on the first and last latitudes.

vor : 2-D vorticity

res=vor(u,v,dx,dy)

Calculates the z component of vorticity of (u,v), using centred differences.

div : 2-D divergence

res=div(u,v,dx,dy)

Calculates the 2-D divergence of (u,v), using centred differences.

def_shear : shear deformation

res=def_shear(u,v,dx,dy)

Calculates the shear (antisymmetric) deformation of (u,v), using centred differences.

def_stretch : stretch deformation

res=def_stretch(u,v,dx,dy)

Calculates the stretch (symmetric) deformation of (u,v), using centred differences.

def_total : total deformation

res=def_total(u,v,dx,dy)

Calculates the total (rotation-independent) deformation of (u,v), using centred differences.

def_angle : deformation angle

res=def_angle(u,v,dx,dy)

Calculates the angle between the x-axis and the dilatation axis of the deformation of (u,v).

isopv_angle : iso-PV line angle

res=isopv_angle(pv,dx,dy)

Calculates the angle between the x-axis and the iso-lines of PV.

beta : angle between dilatation axis and iso-PV lines

res=beta(u,v,pv,dx,dy)

Calculates the angle between the dilatation axis and the iso-lines of PV.

stretch_stir : fractional stretching rate and angular rotation rate of grad(PV)

(stretch,stir)=stretch_stir(u,v,pv,dx,dy)

Returns real arrays, dim (nz,ny,nx):

stretch 
= fractional PV gradient stretching rate 
= 1/|gradPV| * d/dt(|gradPV|)
= gamma, 'stretching rate' (Lapeyre Klein Hua)
= -1/|gradPV| * Fn (Keyser Reeder Reed)
Fn = 0.5*|gradPV|(D-E*cos(2*beta))
   =  1/|gradPV| * F (Markowski Richardson)
stir 
= angular rotation rate of grad(PV) (aka stirring rate)
= d(theta)/dt (Lapeyre Klein Hua)
= 1/|gradPV| * Fs  (Keyser Reeder Reed)
Fs = 0.5*|gradPV|(vort+E*sin(2*beta))

geop_from_montgp : geopotential

res = geop_from_montgp(m,theta,p,dx,dy)

Calculates geopotential (res) from montgomery potential (m), potential temperature (theta) and pressure (p)

rev : PV gradient reversal

(resa,resc,resai,resci,resaiy,resciy,tested) = rev(pv,highenough,latitudes,ddythres,dx,dy)

Gradient reversal: At each (i,j,k) grid point, finds the reversals of PV y-gradient and classes them as c (cyclonic) or a (anticyclonic)

Arguments:
Type Dim Description
pv real (nz,ny,nx) potential vorticity
highenough int*1 (nz,ny,nx) array of flags denoting whether to test the point for reversal
latitudes real (ny) vector of latitudes
ddythres real 0 Cutoff y-gradient for pv
latitudes real (ny) vector of latitudes

Highenough is typically the output of highenough funtion, which returns 1 where the surface is sufficiently above ground level and 0 elsewhere.

ddythres is the cutoff y-gradient for pv. The magnitude of (negative) d(pv)/dy must be above ddythres for reversal to be detected; this applies to revc, reva, revci,revai. Typical value 4E-12.

Returns:
Type Dim Description
revc int*1 (nz,ny,nx) Cyclonic reversal flag (threshold test applied)
reva int*1 (nz,ny,nx) Anticyclonic reversal flag (threshold test applied)
revci real (nz,ny,nx) Cyclonic reversal absolute PV gradient (threshold test applied)
revai real (nz,ny,nx) Anticyclonic reversal absolute PV gradient (threshold test applied)
revciy real (nz,ny,nx) Cyclonic reversal absolute PV y-gradient (no threshold test applied)
revaiy real (nz,ny,nx) Anticyclonic reversal absolute PV y-gradient (no threshold test applied)
tested int*1 (nz,ny,nx) flag to 1 all tested points: where highenough==1 and point not on grid edge

prepare_fft : make data periodic in y for FFT

res = prepare_fft(thedata,dx,dy)

Returns the data extended along complementary meridians (for fft). For each lon, the reflected (lon+180) is attached below so that data is periodic in x and y. NOTE: Input data must be lats -90 to 90, and nx must be even.

Arguments:
Type Dim Description
thedata real (nz,ny,nx) input data
Returns:
Type Dim Description
res real (nz,2*ny-2,nx) output data

sum_kix : sum along k for flagged k-values

(res,nres) = sum_kix(thedata,kix,dx,dy)

Calculates sum along k dimension for k values which are flagged in kix vector (length nz).

Arguments:
Type Dim Description
thedata real (nz,ny,nx) input data
kix int (nz) index flag for summation
Returns:
Type Dim Description
res real (ny,nx) (summed) output data
nres int 0 Number of data summed = sum(kix)

Sum_kix is typically used for calculating seasonal means. To do this, kix is set to 1 for times in the relevant season and 0 elsewhere. After summing res and nres over all years, res/nres gives the mean for the season for all years.

high_enough : flags points which are sufficiently above ground

res = high_enough(zdata,ztest,zthres,dx,dy)

Arguments:
Type Dim Description
zdata real (nz,ny,nx) geopotential of gridpoints
ztest real (1,ny,nx) geopotential of topography
zthres real 0 threshold geopotential height difference


Returns:
Type Dim Description
res int*1 (nz,ny,nx)
3-D flag array set to:
1 if zdata(t,y,x) > (ztest(1,y,x) + zthres)
0 otherwise

contour_rwb : detects RWB events, Riviere algorithm

(beta_a_out,beta_c_out) = contour_rwb(pv_in,lonvalues,latvalues,ncon,lev,dx,dy)

Detects the occurrence of anticyclonic and cyclonic wave-breaking events from a PV field on isentropic coordinates.

Reference: Rivière (2009, hereafter R09): Effect of latitudinal variations in low-level baroclinicity on eddy life cycles and upper-tropospheric wave-breaking processes. J. Atmos. Sci., 66, 1569–1592. See the appendix C.

Arguments:
Type Dim Description
pv_in real (nz,ny,nx) isentropic pv. Should be on a regular lat-lon grid and 180W must be the first longitude. (If 180W is not the first longitude, the outputs will have 180W as the first, so must be rearranged)
lonvalues real (nx) vector of longitudes
latvalues real (ny) vector of latitudes
ncon int 0 number of contours to test, normally 41 or 21
lev real 0 potential temperature of the level
Returns:
Type Dim Description
beta_a_out int (nz,ny,nx) flag array, =1 if anticyclonic wave breaking
beta_c_out int (nz,ny,nx) flag array, =1 if cyclonic wave breaking

v_g : geostrophic velocity

(resx,resy) = v_g(mont,lat,dx,dy)

Calculates geostrophic velocity. Returns zero on equator.

okuboweiss : Okubo-Weiss criterion

res = okuboweiss(u,v,dx,dy)

Calculates Okubo-Weiss criterion lambda_0=1/4 * (sigma^2-omega^2)= 1/4 W, where sigma is total deformation and omega is vorticity.

This is the square of the eigenvalues in Okubo's paper (assumes divergence is negligible).

laccel : Lagrangian acceleration

(resx,resy) = laccel(u,v,mont,lat,dx,dy)

Calculates Lagrangian acceleration on the isentropic surface, based on Montgomery potential.

Arguments:
Type Dim Description
u real (nz,ny,nx) zonal velocity
v real (nz,ny,nx) meridional velocity
mont real (nz,ny,nx) Montgomery potential
lat real (ny) vector of latitudes

accgrad_eigs : Lagrangian acceleration gradient tensor eigenvalues

(respr,respi,resmr,resmi) = accgrad_eigs(u,v,mont,lat,dx,dy)

Calculates eigenvalues of the lagrangian acceleration gradient tensor.

Arguments:
Type Dim Description
u real (nz,ny,nx) zonal velocity
v real (nz,ny,nx) meridional velocity
mont real (nz,ny,nx) Montgomery potential
lat real (ny) vector of latitudes
Returns:
Type Dim Description
respr real (nz,ny,nx) Real part of positive eigenvalue
respi real (nz,ny,nx) Imaginary part of positive eigenvalue
resmr real (nz,ny,nx) Real part of negative eigenvalue
resmi real (nz,ny,nx) Imaginary part of negative eigenvalue
ncon int 0 number of contours to test, normally 41 or 21
lev real 0 potential temperature of the level

dphidt : Lagrangian derivative of compression axis angle

res = dphidt(u,v,mont,lat,dx,dy)

Calculates Lagrangian time derivative of compression axis angle: d(phi)/dt (ref Lapeyre et. al 1999), from deformation and Lagrangian acceleration tensor.

Arguments:
Type Dim Description
u real (nz,ny,nx) zonal velocity
v real (nz,ny,nx) meridional velocity
mont real (nz,ny,nx) Montgomery potential
lat real (ny) vector of latitudes