next up previous index
Next: Discussion of main.f Up: The Code Previous: The Code

main.f

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



Zdzislaw Meglicki
2001-02-26