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., l_r, number of physical grid points in x and y directions, number of harmonics n_x and n_y in x and y directions, requested number of time snapshots and time separation, $\Delta t$ of those snapshots. Requested models for T₀(x, y) and for D(x, y) are read too. The models are selected by choosing an integer in $[1\ldots7]$ for T₀(x, y) and an integer in [1, 2, 3] for D(x, y). Any other integer yields a default model for T₀(x,y) and for D(x,y), so, all in all, we have 8 models for T₀(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
       OPEN( 10, file='diffus.naml', action='read')
       READ( 10, input)
       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