Difference between revisions of "Dynlib"
m (→prepare_fft: make data periodic in y for FFT) 
m (→okuboweiss: OkuboWeiss criterion) 

Line 377:  Line 377:  
<code>res = okuboweiss(u,v,dx,dy)</code>  <code>res = okuboweiss(u,v,dx,dy)</code>  
−  Calculates OkuboWeiss criterion lambda_0=1/4 (sigma^2omega^2)= 1/4 W  +  Calculates OkuboWeiss criterion lambda_0=1/4 * (sigma^2omega^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)  +  This is the square of the eigenvalues in Okubo's paper (assumes divergence is negligible). 
=== ''laccel'': Lagrangian acceleration===  === ''laccel'': Lagrangian acceleration=== 
Revision as of 15:36, 21 January 2013
Contents
 1 Documentation
 2 Obtaining dynlib
 3 Quick start to developing with dynlib
 4 Dynlib functions
 4.1 ddx: partial x derivative
 4.2 ddy: partial y derivative
 4.3 grad: gradient of a scalar
 4.4 lap2: 2D laplacian of a scalar
 4.5 vor: 2D vorticity
 4.6 div: 2D divergence
 4.7 def_shear: shear deformation
 4.8 def_stretch: stretch deformation
 4.9 def_total: total deformation
 4.10 def_angle: deformation angle
 4.11 isopv_angle: isoPV line angle
 4.12 beta: angle between dilatation axis and isoPV lines
 4.13 stretch_stir: fractional stretching rate and angular rotation rate of grad(PV)
 4.14 geop_from_montgp: geopotential
 4.15 rev: PV gradient reversal
 4.16 prepare_fft: make data periodic in y for FFT
 4.17 sum_kix: sum along k for flagged kvalues
 4.18 high_enough: flags points which are sufficiently above ground
 4.19 contour_rwb: detects RWB events, Riviere algorithm
 4.20 v_g: geostrophic velocity
 4.21 okuboweiss: OkuboWeiss criterion
 4.22 laccel: Lagrangian acceleration
 4.23 accgrad_eigs: Lagrangian acceleration gradient tensor eigenvalues
 4.24 dphidt: Lagrangian derivative of compression axis angle
Documentation
The steps necessary to obtain dynlib are described below. A more thorough documentation is compiled in the main documentation page.
Obtaining dynlib
 Copying the source code repository
git clone /Data/gfi/users/tsp065/lib/dynlib.git
 Change into the dynlib folder
cd dynlib
 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 diffcomparisons 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 workinprogress.
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 operate on real arrays with dimension (nz,ny,nx) where nz is number of times or levels, and ny and nx are number of latitudes and longitudes, respectively. Typically, the results for each level or time are computed individually as a 2D slice of the 3D data.
ddx: partial x derivative
res=ddx(dat,dx,dy)
Calculates the partial x derivative of dat, using centred differences. For a nonEWcyclic 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 nonEWcyclic grid, 0 is returned on all edges of the x,y domain. For an EWcyclic grid, 0 is returned on the first and last latitudes.
grad: gradient of a scalar
(resx,resy)=grad(dat,dx,dy)
Calculates the gradient of dat, using centred differences. For a nonEWcyclic grid, 0 is returned on all edges of the x,y domain. For an EWcyclic grid, 0 is returned on the first and last latitudes.
lap2: 2D laplacian of a scalar
res=lap2(dat,dx,dy)
Calculates the 2D laplacian of dat, using centred differences. For a nonEWcyclic grid, 0 is returned on all edges of the x,y domain. For an EWcyclic grid, 0 is returned on the first and last latitudes.
vor: 2D vorticity
res=vor(u,v,dx,dy)
Calculates the z component of vorticity of (u,v), using centred differences.
div: 2D divergence
res=div(u,v,dx,dy)
Calculates the 2D 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 (rotationindependent) deformation of (u,v), using centred differences.
def_angle: deformation angle
res=def_angle(u,v,dx,dy)
Calculates the angle between the xaxis and the dilatation axis of the deformation of (u,v).
isopv_angle: isoPV line angle
res=isopv_angle(pv,dx,dy)
Calculates the angle between the xaxis and the isolines of PV.
beta: angle between dilatation axis and isoPV lines
res=beta(u,v,pv,dx,dy)
Calculates the angle between the dilatation axis and the isolines 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(DE*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 ygradient and classes them as c (cyclonic) or a (anticyclonic)
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 ygradient 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 ygradient 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 4E12.
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 ygradient (no threshold test applied) 
revaiy  real  (nz,ny,nx)  Anticyclonic reversal absolute PV ygradient (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.
Type  Dim  Description  

thedata  real  (nz,ny,nx)  input data 
Type  Dim  Description  

res  real  (nz,2*ny2,nx)  output data 
sum_kix: sum along k for flagged kvalues
(res,nres) = sum_kix(thedata,kix,dx,dy)
Calculates sum along k dimension for k values which are flagged in kix vector (length nz).
Type  Dim  Description  

thedata  real  (nz,ny,nx)  input data 
kix  int  (nz)  index flag for summation 
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)
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 
Type  Dim  Description  

res  int*1  (nz,ny,nx) 

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 wavebreaking events from a PV field on isentropic coordinates.
Reference: RiviÃ¨re (2009, hereafter R09): Effect of latitudinal variations in lowlevel baroclinicity on eddy life cycles and uppertropospheric wavebreaking processes. J. Atmos. Sci., 66, 1569â€“1592. See the appendix C.
Type  Dim  Description  

pv_in  real  (nz,ny,nx)  isentropic pv. Should be on a regular latlon 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 
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: OkuboWeiss criterion
res = okuboweiss(u,v,dx,dy)
Calculates OkuboWeiss criterion lambda_0=1/4 * (sigma^2omega^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.
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 
ncon  int  0  number of contours to test, normally 41 or 21 
lev  real  0  potential temperature of the level 
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.
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 
ncon  int  0  number of contours to test, normally 41 or 21 
lev  real  0  potential temperature of the level 
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.
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 
ncon  int  0  number of contours to test, normally 41 or 21 
lev  real  0  potential temperature of the level 