DIVGEOM is the jOceans module of jLab.

 DIVGEOM  Geometric decomposition of eddy vorticity flux divergence.
 
    [F1,F2,F3,F4]=DIVGEOM(DX,DY,K,L,THETA) returns the geometric decomposition
    of eddy vorticity flux divergence associated with variance ellipses 
    having kinetic energy K, anisotropy L, and orientation THETA.
 
    K, L, and THETA are matrices of the same size.  These are defined on an
    X-Y grid with x oriented in *columns* and y oriented in *rows*.  DX and
    DY are the sampling intervals in the X and Y directions, respectively.
 
    Note that K and L are related to ellipse parameters KAPPA and LAMBDA
    used elsewhere in JLAB by K=KAPPA^2 and L=LAMBDA*KAPPA^2.
 
    F1, F2, F3, and F4 are four different contributions to the eddy 
    vorticity flux divergence, as follows:
  
        F1     Quadratic variations in the linear energy L
        F2     Product of linear variations in THETA and L
        F3     Quadratic variations in the orientation THETA
        F4     Product of linear variations in orientation THETA
 
    For details, see Waterman and Lilly (2015), Geometric decomposition of
    eddy-mean flow feedbacks in barotropic systems, J. Phys. Oceanogr. 
 
    [F1,F2,F3,F4,F]=DIVGEOM(...) also returns the total eddy flux 
    divergence F, calculated directly, with F1+F2+F3+F4 = F apart from 
    numerical error.
 
    By default, DIVGEOM calculates derivates with repeated applications of
    a first central difference.  DIVGEOM(...,'arakawa') alternately uses a 
    modified first central difference appropriate for models that employ an
    Awakawa advection scheme. 
    __________________________________________________________________
 
    Usage: [f1,f2,f3,f4]=divgeom(dx,dy,K,L,theta);
    __________________________________________________________________
    This is part of JLAB --- type 'help jlab' for more information
    (C) 2013--2015 J.M. Lilly --- type 'help jlab_license' for details