Dynlib: Difference between revisions

From gfi
(Added more function descriptions)
mNo edit summary
Line 65: Line 65:
=== ''lap2'': 2-D laplacian of a scalar ===
=== ''lap2'': 2-D laplacian of a scalar ===


(res)=lap2(dat,dx,dy)
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.
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.
Line 71: Line 71:
=== ''vor'': 2-D vorticity ===
=== ''vor'': 2-D vorticity ===


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


Calculates the z component of vorticity of (u,v), using centred differences.  
Calculates the z component of vorticity of (u,v), using centred differences.  
Line 77: Line 77:
=== ''div'': 2-D divergence===
=== ''div'': 2-D divergence===


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


Calculates the 2-D divergence of (u,v), using centred differences.  
Calculates the 2-D divergence of (u,v), using centred differences.  
Line 83: Line 83:
=== ''def_shear'': shear deformation ===
=== ''def_shear'': shear deformation ===


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


Calculates the shear (antisymmetric) deformation of (u,v), using centred differences.  
Calculates the shear (antisymmetric) deformation of (u,v), using centred differences.  
Line 89: Line 89:
=== ''def_stretch'': stretch deformation ===
=== ''def_stretch'': stretch deformation ===


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


Calculates the stretch (symmetric) deformation of (u,v), using centred differences.  
Calculates the stretch (symmetric) deformation of (u,v), using centred differences.  
Line 95: Line 95:
=== ''def_total'': total deformation ===
=== ''def_total'': total deformation ===


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


Calculates the total (rotation-independent) deformation of (u,v), using centred differences.  
Calculates the total (rotation-independent) deformation of (u,v), using centred differences.  
Line 101: Line 101:
=== ''def_angle'': deformation angle ===
=== ''def_angle'': deformation angle ===


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


Calculates the angle between the x-axis and the dilatation axis of the deformation of (u,v).  
Calculates the angle between the x-axis and the dilatation axis of the deformation of (u,v).  
Line 107: Line 107:
=== ''isopv_angle'': iso-PV line angle ===
=== ''isopv_angle'': iso-PV line angle ===


(res)=isopv_angle(pv,dx,dy)
res=isopv_angle(pv,dx,dy)


Calculates the angle between the x-axis and the iso-lines of PV.
Calculates the angle between the x-axis and the iso-lines of PV.
Line 113: Line 113:
=== ''beta'': angle between dilatation axis and iso-PV lines===
=== ''beta'': angle between dilatation axis and iso-PV lines===


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


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


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


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


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


=== ''geop_from_montgp'': geopotential ===
=== ''geop_from_montgp'': geopotential ===
Line 144: Line 152:
=== ''rev'': PV gradient reversal ===
=== ''rev'': PV gradient reversal ===


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


Gradient reversal: At each (i,j,k) grid point, finds the reversals of PV y-gradient.
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)  
                    and classes them as:
                      c (cyclonic)
                      a (anticyclonic)  


Arguments:          pv: Potential vorticity pv(k,j,i) on (time, lat, lon) grid.
Arguments:           
 
                    pv: Potential vorticity pv(k,j,i) on (time, lat, lon) grid.
             highenough: array of flags, highenough(k,j,i) = {0 or 1} (type int*1)
             highenough: array of flags, highenough(k,j,i) = {0 or 1} (type int*1)
                         denoting whether to test the point for reversal. This is
                         denoting whether to test the point for reversal. This is
Line 164: Line 170:
              
              


Returns: int*1 ::  revc,  reva  (reversal flag) (threshold test applied)
Returns:  
 
        int*1 ::  revc,  reva  (reversal flag) (threshold test applied)
           real ::  revci,  revai  (reversal absolute gradient)  
           real ::  revci,  revai  (reversal absolute gradient)  
                                   (threshold test applied)
                                   (threshold test applied)
Line 176: Line 184:
res = prepare_fft(thedata,dx,dy)
res = prepare_fft(thedata,dx,dy)


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


Line 187: Line 194:
Calculates sum along k dimension for k values which are flagged in kix vector (length nz)
Calculates sum along k dimension for k values which are flagged in kix vector (length nz)
   
   
returns:
returns:  
 
  res(ny,nx) - thedata summed over k where kix==1
  res(ny,nx) - thedata summed over k where kix==1
  nres      - sum(kix)
  nres      - sum(kix)
Line 198: Line 206:


Arguments:
Arguments:
  zdata(nz,ny,nx) : geopotentials of all gridpoints
  zdata(nz,ny,nx) : geopotentials of all gridpoints
  ztest(1,ny,nx)  : geopotential of topography
  ztest(1,ny,nx)  : geopotential of topography
Line 203: Line 212:


Returns:  
Returns:  
  res(nz,ny,nx)  : 3-D flag array set to:
  res(nz,ny,nx)  : 3-D flag array set to:
   1 if zdata(t,y,x) > (ztest(1,y,x) + zthres)
   1 if zdata(t,y,x) > (ztest(1,y,x) + zthres)
Line 215: Line 225:
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.
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.
   
   
Inputs:  
Arguments:  
 
     pv_in(nz,ny,nx) : isentropic pv. Should be on a regular lat-lon grid  
     pv_in(nz,ny,nx) : isentropic pv. Should be on a regular lat-lon grid  
                         and 180W must be the first longitude.  
                         and 180W must be the first longitude.  
Line 224: Line 235:
                 ncon : number of contours to test, normally 41 or 21
                 ncon : number of contours to test, normally 41 or 21
                 lev : potential temperature of the level
                 lev : potential temperature of the level
Outputs:
Returns:
 
  beta_a_out(nz,ny,nx) : flag array, =1 if anticyclonic wave breaking  
  beta_a_out(nz,ny,nx) : flag array, =1 if anticyclonic wave breaking  
  beta_c_out(nz,ny,nx) : flag array, =1 if cyclonic wave breaking  
  beta_c_out(nz,ny,nx) : flag array, =1 if cyclonic wave breaking  
Line 249: Line 261:


Arguments:
Arguments:
  u(nz,ny,nx)    : zonal velocity
  u(nz,ny,nx)    : zonal velocity
  v(nz,ny,nx)    : meridional velocity
  v(nz,ny,nx)    : meridional velocity
Line 261: Line 274:


Arguments:
Arguments:
  u(nz,ny,nx)    : zonal velocity
  u(nz,ny,nx)    : zonal velocity
  v(nz,ny,nx)    : meridional velocity
  v(nz,ny,nx)    : meridional velocity
Line 267: Line 281:


Returns:  
Returns:  
  respr(nz,ny,nx)  : Real part of positive eigenvlaue
  respr(nz,ny,nx)  : Real part of positive eigenvlaue
  respi(nz,ny,nx)  : Imaginary part of positive eigenvlaue
  respi(nz,ny,nx)  : Imaginary part of positive eigenvlaue
Line 279: Line 294:


Arguments:
Arguments:
  u(nz,ny,nx)    : zonal velocity
  u(nz,ny,nx)    : zonal velocity
  v(nz,ny,nx)    : meridional velocity
  v(nz,ny,nx)    : meridional velocity

Revision as of 10:47, 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 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 2-D slice of the 3-D data.

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: gradient of a scalar

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

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.

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)

where:

stretch

= fractional PV gradient stretching rate

= 1/|gradPV| * d/dt(|gradPV|)

= gamma, 'stretching rate' (Lapeyre Klein Hua)

= -1/|gradPV| * F_n (Keyser Reeder Reed) where 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| * F_s (Keyser Reeder Reed) where 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:

                   pv: Potential vorticity pv(k,j,i) on (time, lat, lon) grid.
           highenough: array of flags, highenough(k,j,i) = {0 or 1} (type int*1)
                       denoting whether to test the point for reversal. This is
                       typically the output of highenough() funtion, which returns
                       1 where the surface is sufficiently above ground level and
                       0 elsewhere.
            latitudes: vector of latitudes of the pv array
             ddythres: Cutoff y-gradient for pv. The magnitude of (negative) 
                       d(pv)/dy must be above ddythres for reversal to be detected;
                       applies to revc, reva, revci,revai.
            

Returns:

        int*1 ::  revc,   reva   (reversal flag) (threshold test applied)
         real ::  revci,  revai  (reversal absolute gradient) 
                                 (threshold test applied)		
         real ::  revciy, revaiy (reversal absolute y-gradient) 
                                 (no threshold test applied)
         int*1::  tested (flag to 1 all tested points: where highenough==1 
                         and not on edge of grid)

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.

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)

returns:

res(ny,nx) - thedata summed over k where kix==1
nres       - sum(kix)

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:

zdata(nz,ny,nx) : geopotentials of all gridpoints
ztest(1,ny,nx)  : geopotential of topography
zthres          : threshold geopotential height

Returns:

res(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:

    pv_in(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(nx) : vector of longitudes
      latvalues(ny) : vector of latitudes
               ncon : number of contours to test, normally 41 or 21
                lev : potential temperature of the level

Returns:

beta_a_out(nz,ny,nx) : flag array, =1 if anticyclonic wave breaking 
beta_c_out(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 [resx=u_g,resy=v_g] = v_g(mont,lat)

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

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

laccel: Lagrangian acceleration

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

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

Arguments:

u(nz,ny,nx)     : zonal velocity
v(nz,ny,nx)     : meridional velocity
mont(nz,ny,nx)  : Montgomery potential
lat(ny)         : latitude

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:

u(nz,ny,nx)     : zonal velocity
v(nz,ny,nx)     : meridional velocity
mont(nz,ny,nx)  : Montgomery potential
lat(ny)         : latitude

Returns:

respr(nz,ny,nx)   : Real part of positive eigenvlaue
respi(nz,ny,nx)   : Imaginary part of positive eigenvlaue
resmr(nz,ny,nx)   : Real part of negative eigenvlaue
resmi(nz,ny,nx)   : Imaginary part of negative eigenvlaue

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:

u(nz,ny,nx)     : zonal velocity
v(nz,ny,nx)     : meridional velocity
mont(nz,ny,nx)  : Montgomery potential
lat(ny)         : latitude