PROGRAM main
USE parameters; USE diffusion; USE fourier; USE pessl_hpf
IMPLICIT NONE
INTEGER :: n, ix, jx, iy, jy, k, i, j, stat, it, ij
INTEGER :: num_errors
REAL(long) :: t1, t2, delx, dely
INTEGER :: j1, j2
REAL(long), ALLOCATABLE :: dif_x(:), dif_y(:)
REAL(long), ALLOCATABLE :: lambda(:)
REAL(long), ALLOCATABLE :: gap(:)
REAL(long), ALLOCATABLE :: f(:,:), b(:,:), a(:), a_rep(:)
REAL(long), ALLOCATABLE :: x(:,:), ab(:), abt(:,:)
REAL(long), ALLOCATABLE :: xsines(:,:), ysines(:,:), temp(:,:)
INTEGER :: biga_index, num_eigvalues, num_vectors
!HPF$ DISTRIBUTE (*,BLOCK) :: f, b, x, abt, temp
!HPF$ DISTRIBUTE (BLOCK) :: a, ab
CALL init_diffusion
ALLOCATE( dif_x(dif_numx*dif_numy) ); ALLOCATE( dif_y(dif_numy*dif_numx) )
delx = PI/ ( dif_numx + 1.d0); dely = PI/ ( dif_numy + 1.d0)
DO i = 1, dif_numx
DO j = 1, dif_numy
ij = dif_numx*(j-1) + i
dif_x(ij) = delx* i; dif_y(ij) = dely * j
ENDDO
ENDDO
num_errors = 0; n = dif_nx * dif_ny
ALLOCATE( lambda(n), stat=stat); IF( stat .NE. 0) num_errors = num_errors + 1
ALLOCATE( a(n), stat=stat); IF( stat .NE. 0) num_errors = num_errors + 1
ALLOCATE( ab(n), stat=stat); IF( stat .NE. 0) num_errors = num_errors + 1
ALLOCATE( a_rep(n), stat=stat); IF( stat .NE. 0) num_errors = num_errors + 1
ALLOCATE( xsines(dif_npts, dif_nx) , stat=stat ); IF( stat .NE. 0) &
num_errors = num_errors + 1
ALLOCATE( ysines(dif_npts, dif_ny) , stat=stat ); IF( stat .NE. 0) &
num_errors = num_errors + 1
ALLOCATE( temp(dif_npts, dif_ntemps) , stat=stat ); IF( stat .NE. 0) &
num_errors = num_errors + 1
ALLOCATE( f(n,n) , stat=stat ); IF( stat .NE. 0) num_errors = num_errors + 1
ALLOCATE( b(n,n) , stat=stat ); IF( stat .NE. 0) num_errors = num_errors + 1
ALLOCATE( abt(n,dif_ntemps) , stat=stat ); IF( stat .NE. 0) &
num_errors = num_errors + 1
ALLOCATE( x(n,dif_ntemps) , stat=stat ); IF( stat .NE. 0) &
num_errors = num_errors + 1
IF( num_errors .GT. 0 ) THEN
WRITE(*,*) 'Error in allocating arrays in main'; STOP
ENDIF
CALL expand_temp_profile(a)
a_rep = a
CALL get_diffusion_matrix(f)
DO i = 1, dif_nx
DO j = 1, dif_npts
t1 = i*dif_x(j); xsines(j,i) = SQRT(2.d0/pi) * dsin(t1)
ENDDO
ENDDO
DO i = 1, dif_ny
DO j = 1, dif_npts
t2 = i*dif_y(j); ysines(j,i) = SQRT(2.d0/pi) * dsin( t2)
ENDDO
ENDDO
IF( num_errors .GT. 0 ) THEN
WRITE(*,*) 'Error in allocating arrays for syevx in main'; STOP
ENDIF
CALL syevx(f,lambda,'u', z=b)
ab(i) = 0.d0
DO i = 1, n
DO j = 1, n
ab(i) = b(j,i)*a_rep(j) + ab(i)
ENDDO
ENDDO
DO it = 1, dif_ntemps
i = 0
DO i = 1, n
abt(i,it) = EXP( -lambda(i) *(it-1)*dif_delta_t) * ab(i)
ENDDO
ENDDO
CALL gemm(1.d0,b,abt,0.d0,x,'N','N')
DO it = 1, dif_ntemps
DO k = 1, dif_npts
temp(k, it) = 0.d0
ENDDO
ENDDO
DO it = 1, dif_ntemps
DO i = 1, n
iy = (i-1)/dif_nx+1; ix = MOD(i-1,dif_nx)+1
DO k = 1, dif_npts
temp(k,it) = temp(k,it) + xsines(k,ix) * ysines(k,iy) * x(i,it)
ENDDO
ENDDO
ENDDO
WRITE(*,*) ' point # X Y'
DO i = 1, dif_npts
WRITE(*,'(2x, i6, 2x, 2f11.4)') i, dif_x(i), dif_y(i)
ENDDO
WRITE(*,*)
DO k = 1, dif_npts, 6
WRITE(*,*); WRITE(*,'(30X,''Points'')')
WRITE(*,'('' TIME '',6(5x,''#'', i4))') (i, i=k, k+5)
DO i = 1, dif_ntemps
WRITE(*,'(7f10.5)') i*dif_delta_t, (temp(j,i),j=k,MIN(k+5,dif_npts))
ENDDO
ENDDO
DEALLOCATE(xsines); DEALLOCATE(ysines); DEALLOCATE(lambda)
DEALLOCATE(temp); DEALLOCATE(a); DEALLOCATE(abt); DEALLOCATE(a)
DEALLOCATE(f); DEALLOCATE(dif_x); DEALLOCATE(dif_y)
STOP
END PROGRAM main