Discussion of diffusion.f

Right at the beginning of the file there is a keyword

PRIVATE

which means that everything in this module is PRIVATE, i.e., invisible to programs using this module, unless specified otherwise. The otherwise specifications can then be seen below:

  PUBLIC init_diffusion, init_temp, diff_coef
  REAL, PUBLIC :: dif_ly_ratio
  INTEGER, PUBLIC :: dif_nx, dif_ny, dif_npts, dif_ntemps
  INTEGER, PUBLIC :: dif_numx, dif_numy
  REAL(long), PUBLIC :: dif_delta_t

This means that functions init_diffusion, init_temp, and diff_coef are visible to calling programs, as are variables dif_ly_ratio (aka l_r), dif_nx (aka n_x, i.e., the number of harmonics in the x direction), dif_ny (n_y, the number of harmonics in the y direction), dif_npts (n_p, the total number of grid points), dif_ntemps (the total number of time snapshots), dif_numx (number of grid points in the x direction), dif_numy (number of grid points in the y direction), and dif_delta_t (the $\Delta t$ that separates time snapshots).

However, the following two variables:

  INTEGER :: init_f=1, diff_f=1

are private to module diffusion. They are used only by functions init_temp, and diff_coef and don't have to be accessible to any callers. Observe that although init_f and diff_f are private to this module, they are also global within this module, i.e., all functions and subroutines of the module have access to these two variables.

init_diffusion

There is one new element in this subroutine that we haven't covered yet. It is the NAMELIST and its use. This is how it is declared and then used:

NAMELIST /input/ ly_ratio, delta_t, numx, numy, nx, ny, numt, init_f, diff_f
...	 
INQUIRE ( file='diffus.naml', exist=ex)
IF( ex ) THEN
   OPEN( 10, file='diffus.naml', action='read')
   READ( 10, input)
   CLOSE(10)
ENDIF
...

The declaration of NAMELIST specifies certain input format   that is quite unique to Fortran. It is also very convenient, although not always portable.

The way you would have to input your data on file diffus.naml is as follows:

&input 
ly_ratio = 1.0d0,
delta_t = 0.05,
numx = 5,
numy = 6,
nx = 14,
ny = 14,
numt = 40,
init_f = 3,
diff_f = 3
/

You can always check easily, if your data has been read in correctly. To do so simply insert:

WRITE( 6, input )

immediately after the READ( 10, input ) statement. What you will get, when the program runs, should look as follows:

 &INPUT
 LY_RATIO=1.00000000000000000, DELTA_T=0.500000000000000028E-01, 
 NUMX=5, NUMY=6, NX=14, NY=14, NUMT=40, INIT_F=3, DIFF_F=3
 /

The subroutine works as follows. All variables to be read from the namelist are pre-initialised to default values in the declaration part:

INTEGER :: numx=5, numy=5, nx=7, ny=7, numt=20
REAL(long) :: ly_ratio=1.d0, delta_t=0.1

Then the subroutine checks if the input file diffus.naml exists. The answer to that question is returned on a logical variable:

LOGICAL :: ex

If the file exists then values of ly_ratio, delta_t, numx, numy, nx, ny, numt, init_f, and diff_f are read from the namelist. If the file does not exist then we stay with the defaults.   All the variables listed above are PRIVATE. Before the subroutine exits, the data is transferred from   PRIVATE to PUBLIC variables:

    dif_ly_ratio = ly_ratio
    dif_npts = numx*numy
    dif_delta_t = delta_t
    dif_ntemps  = numt
    dif_nx = nx
    dif_ny = ny
    dif_numx = numx
    dif_numy = numy

init_temp

This function returns T₀(x, y) for any x and y. The function is really defined on $[0, \pi]\times[0, \pi]$ and then extended to $[0, 2\pi]\times[0, 2\pi]$ so as to be odd about $\pi$. This ensures that the Fourier decomposition of T₀(x,y) is expressed in terms of the same coefficients a_l(0) both for $[0, \pi]\times[0, \pi]$ and for $[0, 2\pi]\times[0, 2\pi]$, because the $\sin{x}$ function itself is odd about $\pi$.

Here is how this is accomplished:

    isign = 1
    x1 = x
    IF( x .GT. pi ) THEN
       isign = -isign
       x1 = twopi - x
    ENDIF
    y1 = y
    IF( y .GT. pi ) THEN
       isign = -isign
       y1 = twopi - y
    ENDIF

Now we have the definitions of eight models of T₀(x, y) (in truth only 7, because model 5 is the same as model 1) implemented using the SELECT CASE statement (which is much the same as SWITCH in C):

1.

$T_0(x, y) = x(\pi - x)y(\pi - y)$

    CASE (1)
       init_temp = isign*(x1*(pi-x1))*y1*(pi-y1)

2.

$T_0(x, y) = x(\pi - x)y(\pi - y) y$

    CASE (2)
       init_temp = isign*(x1*(pi-x1))*y1*(pi-y1)*y1

3.

$T_0(x, y) = x(\pi - x)y(\pi - y)x$

    CASE (3)
       init_temp = isign*(x1*(pi-x1))*y1*(pi-y1)*x1

4.

$T_0(x, y) = x(\pi - x)y(\pi - y)x y$

    CASE (4)
       init_temp = isign*(x1*(pi-x1))*y1*(pi-y1)*x1*y1

5.

$T_0(x, y) = x(\pi - x)y(\pi - y)$

    CASE (5)
       init_temp = isign*(x1*(pi-x1))*y1*(pi-y1)

This model is really the same as model 1.

6.

$T_0(x, y) = \left(x(\pi - x)\right)^2 y (\pi - y)$

    CASE (6)
       init_temp = isign*(x1*(pi-x1))**2 *y1*(pi-y1)

7.

$T_0(x, y) = x(\pi - x) \left(y(\pi - y)\right)^2$

    CASE (7)
       init_temp = isign*(x1*(pi-x1))*(y1*(pi-y1))**2

8.

If init_f is not in $[1\ldots7]$ then the default distribution is assumed, which is $T_0(x, y) = \sin{x}\, \sin{y}$

    CASE default
       init_temp = isign*SIN(x1)*SIN(y1)

diff_coef

Function diff_coef defines D(x,y). The diffusion coefficient is defined to be even about $\pi$. The reason for this is that we are really considering the following integral:

$$\begin{eqnarray*}&&\frac{4}{\pi^2}\int_0^\pi\int_0^\pi\big( kk'D \cos(kx)\cos(k'... ...(jy)\cos(j'y)\sin(kx)\sin(k'x)\big) \,\textrm{d}x\,\textrm{d}y \end{eqnarray*}$$

and then extending it to $[0, 2\pi]$ by absorbing 4 into $\int_0^{2\pi}\int_0^{2\pi}$. Now, this is OK, but only if D(x,y) is even about $\pi$, because $\sin(kx)\sin(k'x)$ and $\sin(jy)\sin(j'y)$ are both even about $\pi$.

So, here is how D(x,y) is made even about $\pi$:

    x1 = x
    IF( x .GT. pi )  x1 = twopi - x
    y1 = y
    IF( y .GT. pi ) y1 = twopi - y

The three models and a default model for D(x,y) are as follows:

1.

$D(x,y) = \frac{1}{2} + \frac{x + y}{2\pi}$

    CASE (1)
       diff_coef = .5d0 + (x1 + y1) / (2 * twopi)

2.

$D(x,y) = \frac{(1+x)(\pi - x + 1)(y + \pi)}{3\pi}$

    CASE (2)
       diff_coef = ((1.d0 + x1)*(pi - x1 + 1.d0)*(y1 + pi))/ 3*pi

3.

$D(x,y) = \frac{(y+\pi)\pi}{(\pi + x)(2\pi - x)}$

    CASE (3)
       diff_coef = (y1 + pi) * pi/((pi + x1) * (2* pi - x1))

4.

The default is simply D(x,y) = 1:

    CASE default
       diff_coef = 1.d0