Potential solver io4

From OpenGGCM

io4 solves the ionosphere potential equation nab*sig*nab*Phi = J_par*sin(I). See papers for details. Solution is on entire sphere but internally each hemisphere is solved separately. Phi=0 at the equator. Uses finite element (Galerkin) method. The elements are Hermite polyniminals in theta, non-equally spaced, and Fourier (sin,cos) in azimuth. Leads to a dense matrix which is solved by a direct method.

More specifically, it solves:

c
c   galerkin pseudospectral routine to solve the ionosphere
c   potential equation
c
c     d/dx( a(x,y) * df/dx + b(x,y) * df/dy )
c   + d/dy( c(x,y) * df/dx + d(x,y) * df/dy )  =  rhs
c
c   with boundary conditions:  f(0,.)=0, df(xx2,.)/dx=0, periodic y
c
c   basis functions are:
c
c   B(i,j,x,y) = H(mod(i-1,4),x) * F(j,y) , i=1,..,4*nx;  j=1,...,2*ny
c
c   with:   H(i,x) : Hermite polynominal with support t1,..,t2
c                    cubics with continous first derivative
c           F(j,y) : sin((((j-1)/2)+1)*y), j odd
c                    cos( ((j-1)/2   )*y), j even
c

   subroutine io4(ijob,nx,ny,ngord,npany,sigp,sigh,xjpa,pot,nphi,
  +               nthe,wk,nwk,t1,t2,t3)

   ijob=0:   initialize
   ijob=1:   set up matrix and factor it
   ijob=3:   solve equation with matrix factored
   ijob=4:   evaluate solution on reg grid xjpa

   nx:       number of panels in latitude for hermite FE
   ny:       number of modes in azimuth

   ngord:    order of gauss integration
   npany:    no of panels in y for integrations

   sigp:     array holding pedersen conductivity
   sigh:     array holding hall conductivity
   xjpa:     array holding parallel current, for ijob=4: potential

   nphi,nthe:  dimensions of sigp,sigh,xjpa
 
   wk:       workspace, must not be touched between calls.
 
   nwk:      size of workspace, will stop if not enough.

   t1,t2,t3:  nphi*nthe  work arrays