! ! Compile with ! f95 -O4 -r8 -o burg burg.f90 ! program burg ! ! Reduce typographical errors with implicit none: ! implicit none ! ! For convenience, define i and pi ! complex, parameter :: zi = ( 0.0 , 1.0 ) real, parameter :: pi = 3.1415926535897931 ! ! Defer specification of problem size until runtime ! complex, dimension(:), allocatable :: a, astar, aold, nl real, dimension(:), allocatable :: gam ! ! Declare a few control variables ! real :: dt, nu, gamma, pert, time, fac integer :: nsteps, n, kmax, nout, kinj integer :: k ! ! Read inputs, initialize arrays ! call init ! ! Main loop ! do n = 1, nsteps call step ! ! Every nout steps, do some output to make time-dependent graphs later. ! if (mod(n,nout) == 0) call output if (mod(n,nout) == 0) call final_output end do ! ! Write final info to a file ! call final_output ! ! This is the end of the main program block. ! contains ! ! Subroutines that are "contained" have access to all ! variables declared in the containing procedure, so ! common blocks are not needed. ! subroutine init ! ! Declare some local variables. ! real :: x, y ! ! Define the input namelist ! namelist /input/ dt, nu, kinj, gamma, kmax, nsteps, pert, nout, fac ! ! Open the input file, read the namelist, and close it back up. ! open (unit=10, file='burg.in', status='old') read (10, input) close (10) ! ! Allocate storage for the calculation. Note that we do not need to keep ! information for k < 0 because the reality condition determines them from ! the relation a(-k) = conjg(a(k)). The total amount of information is ! the same as in the real-space representation, since the main arrays are ! complex. ! allocate (a (0:kmax)) allocate (aold (0:kmax)) allocate (astar(0:kmax)) allocate (nl (0:kmax)) allocate (gam (0:kmax)) ! ! Set allocated variables to zero initially. ! a = 0. aold = 0. astar = 0. nl = 0. gam = 0. ! ! Put in initial condition. ! if (pert > 0.) then do k = 1, kmax call random_number(x) call random_number(y) aold(k) = cmplx(x,y) * pert end do else aold(1) = -pert end if ! ! Set the growth rate and k for the injected energy. ! gam(kinj) = gamma ! ! Set time = 0. ! time = 0. ! ! Open unit 11 for later writes. ! open (unit=11, file='burg.out') end subroutine init subroutine step ! ! This is a second-order Runge-Kutta scheme in time. ! ! ! Take half-step first with aold on the RHS: ! ! Linear terms: ! do k = 0, kmax astar(k) = aold(k) * (1. - 0.5*nu*k**2*dt + 0.5*gam(k)*dt) end do ! ! Nonlinear term: ! call nonlinear (aold, nl) do k = 0, kmax astar(k) = astar(k) - 0.5*zi*dt*nl(k)*fac end do ! ! Then take full step from aold, using astar on the RHS: ! do k = 0, kmax a(k) = aold(k) - nu*k**2*dt*astar(k) + gam(k)*dt*astar(k) end do call nonlinear (astar, nl) do k = 0, kmax a(k) = a(k) - zi*dt*nl(k)*fac end do ! ! Update the "old" value ! aold = a ! ! Update the time for diagnostics. ! time = time + dt end subroutine step subroutine nonlinear (a, nl) ! ! Note that the variable "a" is locally defined in this subroutine. ! complex, dimension(0:) :: a, nl complex :: aptmp, aqtmp integer :: p, q ! ! Calculate the nonlinear term from a given array of "a" values ! nl = 0. do p = -kmax, kmax do q = -kmax, kmax if (abs(p+q) > kmax) cycle if (p+q < 0) cycle if (p < 0) then aptmp = conjg(a(-p)) else aptmp = a(p) end if if (q < 0) then aqtmp = conjg(a(-q)) else aqtmp = a(q) end if nl(p+q) = nl(p+q) + p * aptmp * aqtmp end do end do end subroutine nonlinear subroutine output ! ! Put quantities that you want time histories of here: ! write(11,*) time, real(a(1)), aimag(a(1)), real(sum(a*conjg(a))), & real(a(1)*conjg(a(1))) end subroutine output subroutine final_output ! ! These are output statements called at the end of the program. ! complex, allocatable, dimension(:) :: u integer :: j logical :: first = .true. open (unit=12, file='burg.out2') do k = 0, kmax write(12,*) k, real(a(k)), aimag(a(k)), real(a(k)*conjg(a(k))) end do if (first) then open (unit=13, file='burg.x') write(13,*) '# x u(x)' first = .false. else write(13,*) end if ! ! Construct real-space representation of u from spectral coefficients: ! allocate (u(-kmax:kmax)) u = 0. do j = -kmax, kmax ! ! These loops could be combined: ! do k = 1, kmax u(j) = u(j) + a(k)*exp(zi*k*pi*j/kmax) end do do k = 1, kmax u(j) = u(j) + conjg(a(k))*exp(-zi*k*pi*j/kmax) end do end do do j = -kmax, kmax write(13,*) pi*j/kmax, real(u(j))/2. end do deallocate (u) end subroutine final_output end program burg