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 a₀^k = a^k(0), which correspond to T₀(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