Add QGPV Dipole solution as initial condition

[JeffreyEarly]

Note to self: check the FPSimulator code, look at the addDipole function. This appears to set PV = -beta*y inside a circle, and zero outside. It then does PV inversion to get the stream function.

[JeffreyEarly]

Nope, here's the relevant code,

        double c = 1;
    double a = 2/sqrt(2);
    double q = a * sqrt( 1/c + 1);
    double k = 4.0732;

    double D1 = -c * a / bessk(1, q);
    double B1 = c * (1/c + 1) * a * a * a / ( k * k * bessj(1, k) );

    double kovera = k/a;
    double linearCoeff = ( 1.0 + ( kovera * kovera + 1 ) * c )/(kovera*kovera);

    double xdomain = [self xmaxKM] - [self xminKM];
    double ydomain = [self ymaxKM] - [self yminKM];

    int i,j;
    for ( i=0; i< [self xgrid]; i++ ) {
        double x = xdomain * ( (double) i) / ((double) [self xgrid]) + [self xminKM];
        for ( j=0; j<[self ygrid]; j++ ) {
            double y = ( ydomain * ( (double) j) / ((double) [self ygrid]) + [self yminKM]);
            double r2 = fmod( fabs(x - self.gaussianXcenter), xdomain)*fmod( fabs(x - self.gaussianXcenter), xdomain)+fmod( fabs(y - self.gaussianYcenter), ydomain)*fmod( fabs(y - self.gaussianYcenter), ydomain);
            double r = sqrt( r2 )/(self.lengthScale / 1000);
            double theta = atan2( y - self.gaussianYcenter,x - self.gaussianXcenter);
            if ( r <= a ) {
                psiMatrix[i][j] += (B1 * bessj(1, kovera*r) - linearCoeff * r)*sin(theta) / self.jacobianCoefficient;
            } else {
                psiMatrix[i][j] += (D1*bessk(1, sqrt(1/c+1) * r))*sin(theta) / self.jacobianCoefficient;
            }
        }
    }

At least the interior part is equation 5.3 in Flierl 1980,