Problem II.5.2: Instructions and useful hints
==================================================================

1. Format of Data file
   -------------------

The accompanying file (btacho_1.dat) contains the steady-state
solution for the magnetized tachocline problem of section II.5.4,
for parameter values Rv=Rm=1000 and B0=1 Gauss.

The solution was obtained on a 201 x 41 grid in [r,mu], where
mu=cos(theta).

The data was generated by the following FORTRAN-77 statements:

      parameter (nm=41,nr=201)
      dimension omega(nm,nr),bphi(nm,nr),xmu(nm),r(nr)

      write(*,10) nr,nm
      write(*,11) (xmu(i),i=1,nm)
      write(*,11) (r(i)  ,i=1,nr)
      do 1 j=1,nr
	 write(*,11) (omega(i,j),i=1,nm)
	 write(*,11) ( bphi(i,j),i=1,nm)
    1 continue
   10 format(2i5)
   11 format(1p6e2.4)

If you can't decipher the above, then here's the long version:

The first line contains the two integers defining the size
of the spatial grid

The next series of lines contain 41 floating point numbers,
defining the mu-part of the 2-D grid

The next series of lines contain 201 floating point numbers,
defining the r-part of the 2-D grid

The next series of lines contain 41 floating point numbers,
corresponding to the values of Omega(r,mu) as a function of mu
for the first depth value r(1)

The next series of lines contain 41 floating point numbers,
corresponding to the values of B_phi(r,mu) as a function of mu
for the first depth value r(1)

The next series of lines contain 41 floating point numbers,
corresponding to the values of Omega(r,mu) as a function of mu
for the second depth value r(2)

...and so on, all the way to the last block of lines, which
contains the 41 floating point numbers corresponding to the
values of B_phi(r,mu) as a function of mu
for the last depth value r(201)

2. Finite difference estimates for partial derivatives
   ---------------------------------------------------

Let the index i=1,...,41 define the location of a node along
the mu-direction of the grid, and the index j=1,...,201 define
the location of the node along the r-direction of the grid.
We will write as Omega(i,j) the value of Omega at this node,
with the same notation for B_phi. 

Second-order centered finite difference formulae for partial
derivatives of Omega with respect to mu or r at node (i,j)
are as follows:

d Omega         Omega(i+1,j)-Omega(i-1,j)
-------(i,j) =  -------------------------
 d mu                mu(i+1)-mu(i-1)

d Omega         Omega(i,j+1)-Omega(i,j-1)
-------(i,j) =  -------------------------
 d r                 r(j+1)-r(j-1)

d2 Omega            Omega(i+1,j)-2*Omega(i,j)+Omega(i-1,j)
--------(i,j) =  4* --------------------------------------
 d mu2                      (mu(i+1)-mu(i-1))**2

d2 Omega            Omega(i,j+1)-2*Omega(i,j)+Omega(i,j-1)
--------(i,j) =  4* --------------------------------------
 d r2                        (r(j+1)-r(j-1))**2



Note: the above formulae are strictly valid only for equally spaced
spatial grid. The mu part of your grid is equally spaced, but not
the r; however, the grid spacing in the r-direction varies slowly
and smoothly, so the use of the above formulae remains warranted
Note however that such centered expression cannot be used along
the boundaries.