Hint 9 above is actually somewhat misleading. At some stage I got confused myself by those sizing manipulations, which looked quite puzzling to me (and for a good reason, as it turned out: they have lose ends).

First, the Fourier tranform of function *D*(*x*,*y*)is going to be real, because *D*(*x*,*y*)is *even* about .
You can see it as follows:

If

`expand_temp_profile`

and in `get_diffusion_matrix`

, subroutine `fft`

is called
in the same way, i.e.,CALL fft(df,scale=scale_f,transpose='N') ... CALL fft(atmp,scale=scale_f,transpose='N')This form of calling

`fft`

and the output, which will be written `fft`

and
then extract the real part from the result.
The reasons for oversizing the transform have to do with negative wave numbers.

Still the first part of the hint, referring you to the discussion of
subroutine `fft`

in the PESSL manual was fine, and if you read
that discussion carefully, you would have noticed that both
invocations to `fft`

in our program
are done in the same, complex-to-complex, way. Consequently the
sizing differences associated with special calling conventions
for a complex-to-real case wouldn't apply here anyway.

It is instructive to have a look at the Fourier transform of the diffusion coefficient for the following input:

&input ly_ratio = 1.0d0, delta_t = 0.2, numx = 5, numy = 5, nx = 7, ny = 7, numt = 10, init_f = 8, diff_f = 2, /Diffusion coefficient model 2 is

The values of matrix

`df`

before calling `fft`

look as follows:df( 1, 1) = (0.165E+02, 0.000E+00) df( 1, 2) = (0.175E+02, 0.000E+00) df( 1, 3) = (0.184E+02, 0.000E+00) df( 1, 4) = (0.194E+02, 0.000E+00) df( 1, 5) = (0.204E+02, 0.000E+00) df( 1, 6) = (0.214E+02, 0.000E+00) df( 1, 7) = (0.223E+02, 0.000E+00) df( 1, 8) = (0.233E+02, 0.000E+00) ... df(31,25) = (0.223E+02, 0.000E+00) df(31,26) = (0.214E+02, 0.000E+00) df(31,27) = (0.204E+02, 0.000E+00) df(31,28) = (0.194E+02, 0.000E+00) df(31,29) = (0.184E+02, 0.000E+00) df(31,30) = (0.175E+02, 0.000E+00) df(31,31) = (0.165E+02, 0.000E+00)On exit from

`fft`

the non-zero coefficients look as follows:df( 2, 1) = (0.339E+00, 0.158E-16) df( 2, 3) = (0.386E-01, -.640E-16) df( 2, 5) = (0.146E-01, 0.166E-17) df( 2, 7) = (0.809E-02, -.433E-19) df( 2, 9) = (0.545E-02, -.206E-16) df( 2,11) = (0.418E-02, 0.673E-16) df( 2,13) = (0.355E-02, -.405E-16) df( 2,15) = (0.329E-02, -.555E-16) df( 2,17) = (0.329E-02, -.665E-16) df( 2,19) = (0.355E-02, -.736E-16) df( 2,21) = (0.418E-02, 0.347E-16) df( 2,23) = (0.545E-02, -.394E-16) df( 2,25) = (0.809E-02, 0.717E-17) df( 2,27) = (0.146E-01, 0.396E-16) df( 2,29) = (0.386E-01, -.342E-16) df( 2,31) = (0.339E+00, 0.231E-16) df( 4, 1) = (0.880E-01, -.301E-16) df( 4, 3) = (0.100E-01, 0.105E-15) df( 4, 5) = (0.381E-02, -.190E-16) df( 4, 7) = (0.210E-02, 0.445E-16) df( 4, 9) = (0.142E-02, -.327E-16) df( 4,11) = (0.109E-02, 0.409E-16) df( 4,13) = (0.924E-03, 0.615E-16) df( 4,15) = (0.854E-03, -.151E-16) df( 4,17) = (0.854E-03, -.167E-16) df( 4,19) = (0.924E-03, 0.642E-16) df( 4,21) = (0.109E-02, 0.461E-16) df( 4,23) = (0.142E-02, -.322E-16) df( 4,25) = (0.210E-02, 0.436E-16) df( 4,27) = (0.381E-02, -.181E-16) df( 4,29) = (0.100E-01, 0.104E-15) df( 4,31) = (0.880E-01, -.316E-16) ...There are quite a few of them, and they are all indeed real (with ``zero'' evaluating to something like ). The input function here is not a simple trigonometric function, so the Fourier spectrum gets scattered all over the place. But we can make things cleaner by adding the following model for the diffusion coefficient:

CASE (4) diff_coef = COS(x1) * COS(y1)and changing the input namelist to

&input ly_ratio = 1.0d0, delta_t = 0.2, numx = 5, numy = 5, nx = 7, ny = 7, numt = 10, init_f = 8, diff_f = 4, /This is a very strange kind of a diffusion coefficient, because it is negative in places. But the Fourier transform of the diffusion coefficient now has terms different from zero only for:

df( 1, 1) = (0.250E+00, -.573E-17) df( 1,31) = (0.250E+00, 0.187E-17) df(31, 1) = (0.250E+00, -.632E-19) df(31,31) = (0.250E+00, 0.526E-17)The scale factor here is different than it was for the initial temperature. There we had:

scale_f = -twopi / REAL( nx*ny,long)whereas here it is

scale_f = 1.d0/ REAL(nx*ny, long)which means that the result should be smaller, i.e.,

Last, you may still wonder, why the Fourier coefficients obtained
both for the initial temperature distribution and for the
diffusion coefficient are positive in this case, even though
the scale factor for the initial temperature distribution is
negative. The reason for this is that *T*_{0}(*x*,*y*) is odd about
,
whereas *D*(*x*,*y*) is even about .
This means that
the integral

will select the imaginary component this time. But, the same will happen in the