Next: Discussion of diffusion.f Up: The Code Previous: param.f

### diffusion.f

This file, which lives on

```/afs/ovpit.indiana.edu/common/www/htdocs/gustav/diffusion/diffusion.f
```
contains one simple subroutine and two simple functions.

The subroutine we have already encountered in the `main`. It is `init_diffusion`, and its task is to read the input that defines basic parameters for the run, i.e., lr, number of physical grid points in x and y directions, number of harmonics nx and ny in x and y directions, requested number of time snapshots and time separation, of those snapshots. Requested models for T0(x, y) and for D(x, y) are read too. The models are selected by choosing an integer in for T0(x, y) and an integer in [1, 2, 3] for D(x, y). Any other integer yields a default model for T0(x,y) and for D(x,y), so, all in all, we have 8 models for T0(x,y) and 4 for D(x,y).

The two functions `init_temp` and `diff_coef` implement those models.

Here's the listing of `diffusion.f`

```MODULE diffusion

USE parameters
IMPLICIT NONE

PRIVATE

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
INTEGER :: init_f=1, diff_f=1

CONTAINS

SUBROUTINE init_diffusion

NAMELIST /input/ ly_ratio, delta_t, numx, numy, nx, ny, numt, &
init_f, diff_f
INTEGER :: numx=5, numy=5, nx=7, ny=7, numt=20
REAL(long) :: ly_ratio=1.d0, delta_t=0.1
LOGICAL :: ex

INQUIRE ( file='diffus.naml', exist=ex)
IF( ex ) THEN
CLOSE(10)
ENDIF
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
RETURN

END SUBROUTINE init_diffusion

FUNCTION init_temp(x, y)
REAL(long), INTENT(in) :: x, y
REAL(long) :: init_temp

INTEGER :: isign
REAL(long) :: x1, y1

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
SELECT CASE (init_f)
CASE (1)
init_temp = isign*(x1*(pi-x1))*y1*(pi-y1)
CASE (2)
init_temp = isign*(x1*(pi-x1))*y1*(pi-y1)*y1
CASE (3)
init_temp = isign*(x1*(pi-x1))*y1*(pi-y1)*x1
CASE (4)
init_temp = isign*(x1*(pi-x1))*y1*(pi-y1)*x1*y1
CASE (5)
init_temp = isign*(x1*(pi-x1))*y1*(pi-y1)
CASE (6)
init_temp = isign*(x1*(pi-x1))**2 *y1*(pi-y1)
CASE (7)
init_temp = isign*(x1*(pi-x1))*(y1*(pi-y1))**2
CASE default
init_temp = isign*SIN(x1)*SIN(y1)
END SELECT
RETURN

END FUNCTION init_temp

FUNCTION diff_coef(x, y)
REAL(long), INTENT(in) :: x, y
REAL(long) :: diff_coef

REAL(long) :: x1, y1

x1 = x
IF( x .GT. pi )  x1 = twopi - x
y1 = y
IF( y .GT. pi ) y1 = twopi - y
SELECT CASE (diff_f)
CASE (1)
diff_coef = .5d0 + (x1 + y1) / (2 * twopi)
CASE (2)
diff_coef = ((1.d0 + x1)*(pi - x1 + 1.d0)*(y1 + pi))/ 3*pi
CASE (3)
diff_coef = (y1 + pi) * pi/((pi + x1) * (2* pi - x1))
CASE default
diff_coef = 1.d0
END SELECT
RETURN
END FUNCTION diff_coef

END MODULE diffusion
```

Zdzislaw Meglicki
2001-02-26