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

### fourier.f

• This file, fourier.f, contains implementation of two major subroutines used by main. These are subroutine get_diffusion_matrix, which evaluates matrix F and subroutine expand_temp_profile, which evaluates coefficients a0k = ak(0), which correspond to T0(x, y).
• There are also two internal (PRIVATE) subroutines defined on that file, subroutine factor_nodes and function minpower2. These are used internally by get_diffusion_matrix and by expand_temp_profile, and do not have to be visible to external programs.

Here's the listing:

MODULE fourier

USE parameters
USE diffusion
USE pessl_hpf
USE hpf_library
IMPLICIT NONE

PRIVATE

PUBLIC expand_temp_profile, get_diffusion_matrix

PRIVATE minpower2, factor_nodes

CONTAINS

SUBROUTINE get_diffusion_matrix(f)
REAL(long), INTENT(out) :: f(:,:)
!HPF\$ DISTRIBUTE *(*,BLOCK) :: f
COMPLEX(long), ALLOCATABLE :: df(:,:)
!HPF\$ DISTRIBUTE (*,BLOCK) :: df
INTEGER, ALLOCATABLE :: ixi(:), iyi(:)
REAL(long) :: scale_f
INTEGER :: nx, ny, ix, iy, ixp, iyp, istat
INTEGER :: ixdiff, iydiff, num_errors
INTEGER :: naux1, naux2, i, j, factor1, factor2, idum
INTEGER :: pn_fac(5)   ! prime factors of number of nodes

pn_fac = 0
CALL factor_nodes(pn_fac)
factor1 = 3**pn_fac(2) * 5**pn_fac(3) * 7**pn_fac(4) * &
11**pn_fac(5)
factor2 = (4*(dif_nx+1) + factor1 -1)/factor1
nx = minpower2( factor2,idum) * factor1
factor2 = (4*(dif_ny+1) + factor1 -1)/factor1
ny = minpower2( factor2,idum) * factor1
scale_f = 1.d0/ REAL(nx*ny, long)
ALLOCATE(df(nx,ny))
DO j = 1, ny
DO i =1,  nx
df(i,j) = diff_coef((twopi*(i-1))/nx,(twopi*(j-1))/ny)
ENDDO
ENDDO
CALL fft(df,scale=scale_f,transpose='N')
num_errors=0
ALLOCATE(ixi(dif_nx*dif_ny), stat=istat)
IF( istat .NE. 0 ) num_errors = num_errors + 1
ALLOCATE(iyi(dif_nx*dif_ny), stat=istat)
IF( istat .NE. 0 ) num_errors = num_errors + 1
DO ix = 1, dif_nx
DO iy = 1, dif_ny
i = (iy-1)* dif_nx + ix
ixi(i) = ix
iyi(i) = iy
ENDDO
ENDDO
DO j = 1, dif_nx*dif_ny
iyp = iyi(j)
ixp = ixi(j)
DO i = 1, dif_nx*dif_ny
iy = iyi(i)
ix = ixi(i)
ixdiff = iabs(ix-ixp) + 1
iydiff = iabs(iy-iyp) + 1
f(i,j) = ( ( ix*ixp + iy*iyp*dif_ly_ratio ) *   &
REAL(df(ixdiff, iydiff)    - df(ix+ixp+1,iy+iyp+1)) &
+     ( ix*ixp - iy*iyp*dif_ly_ratio ) * &
REAL(df(ix+ixp+1,iydiff )  - df(ixdiff,iy+iyp+1)))
ENDDO
ENDDO
DEALLOCATE(df)
DEALLOCATE(ixi)
DEALLOCATE(iyi)
RETURN
END SUBROUTINE get_diffusion_matrix

SUBROUTINE expand_temp_profile(a)
REAL(long), INTENT(out) :: a(:)
!HPF\$ DISTRIBUTE *(BLOCK) :: a
COMPLEX(long), ALLOCATABLE :: atmp(:,:)
!HPF\$ DISTRIBUTE (*,BLOCK) :: atmp
INTEGER :: i,j, nx, ny, istat, naux1, naux2
INTEGER :: idum, jb, jx, jy
REAL(long) :: x, y, scale_f

nx = minpower2( 2*(dif_nx+1), idum)
ny = minpower2( 2*(dif_ny+1), idum)
scale_f = -twopi / REAL( nx*ny,long)
ALLOCATE(atmp(nx,ny), stat=istat )
IF( istat .NE. 0 )  THEN
WRITE(*,*) 'allocation failure in expand_temp_profile'
STOP
ENDIF
DO i = 1, nx
DO j = 1, ny
atmp(i,j) = init_temp((twopi*(i-1))/nx, (twopi*(j-1))/ny)
ENDDO
ENDDO
CALL fft(atmp,scale=scale_f,transpose='N')
DO j =  1, dif_nx * dif_ny
jx = MOD(j-1, dif_nx) + 2
jy = (j-1) / dif_nx + 2
a(j) = REAL(atmp(jx,jy),long)
ENDDO
DEALLOCATE(atmp)
RETURN
END SUBROUTINE expand_temp_profile

SUBROUTINE factor_nodes(pn_fac)
INTEGER    :: pn_fac(:)
INTEGER n1, n2, l2

n2 = NUMBER_OF_PROCESSORS()
n1 = n2/11
IF( n1*11 .EQ. n2) THEN
pn_fac(5) = 1
n2 = n1
ENDIF
n1 = n2/7
IF( n1*7 .EQ. n2) THEN
pn_fac(4) = 1
n2 = n1
ENDIF
n1 = n2/5
IF( n1*5 .EQ. n2) THEN
pn_fac(3) = 1
n2 = n1
ENDIF
n1 = n2/3
IF( n1*3 .EQ. n2) THEN
IF ( (n1/3)*3 .EQ. n1 ) THEN
pn_fac(2) = 2
n2 = n1/3
ELSE
pn_fac(2) = 1
n2 = n1
ENDIF
ENDIF
n1 = minpower2(n2,l2)
pn_fac(1) = l2
IF( n1 .NE. n2) THEN
WRITE(*,*) 'Invalid number of nodes, it must have the form:'
WRITE(*,*) '2**n * 3**m * 5**i * 7**j * 11**k, where '
WRITE(*,*) ' n >= 0, 0<=m<=2 and 0<= i,j,k <=1 '
WRITE(*,*) ' choose a power of 2 for best performance'
STOP
ENDIF
END SUBROUTINE factor_nodes

FUNCTION minpower2( n, log2n )
INTEGER n, minpower2, log2n
INTEGER m, i
m=n

IF( n < 0 ) WRITE(*,*) 'n cannot be negative'
powerloop: DO i= 1, BIT_SIZE(n)
m = m / 2
IF( m .EQ. 0 ) EXIT powerloop
ENDDO powerloop
IF ( 2**(i-1) .NE. n ) THEN
IF( 2**i < 0) WRITE(*,*) 'n too large'
log2n = i
minpower2 = 2**i
ELSE
log2n = i-1
minpower2 = n
ENDIF
RETURN
END FUNCTION minpower2

END MODULE fourier

Zdzislaw Meglicki
2001-02-26