antenna_amplitudes Subroutine

  subroutine antenna_amplitudes (apar)
    use mp, only: broadcast, proc0
    use gs2_time, only: code_dt, code_time
    use kt_grids, only: naky, ntheta0, enforce_reality
    use theta_grid, only: theta, ntgrid, gradpar
    use ran, only: ranf
    use constants, only: zi
    
    complex, dimension (-ntgrid:,:,:), intent(out) :: apar
    complex :: force_a, force_b
    real :: dt, sigma, time
    integer :: i, it, ik

    apar=0.

    if (no_driver) return

    if (.not. allocated(apar_new)) then
       allocate (apar_new(-ntgrid:ntgrid,ntheta0,naky)) ; apar_new = 0.
       allocate (apar_old(-ntgrid:ntgrid,ntheta0,naky)) 
    end if

    apar_old = apar_new

    dt = code_dt
    time = code_time

    if (time > t0) then       
       wtmp = w_antenna + w_dot*(time-t0)
    else
       wtmp = w_antenna
    end if

! GGH fixed a bug with the frequency sweeping here.  11.17.05
    do i = 1, nk_stir
       w_stir(i) = gradpar(0)*abs(kz_stir(i))*sqrt(1./beta_s) * wtmp
    end do

    apar = 0.

    do i=1, nk_stir
       
       it = 1 + mod(kx_stir(i)+ntheta0,ntheta0)
       ik = 1 + mod(ky_stir(i)+naky,naky)

       if (proc0) then
          !GGH Decorrelation
          sigma = sqrt(12.*abs(aimag(w_stir(i)))/dt)*amplitude
          
          force_a = cmplx(ranf() - 0.5, ranf() - 0.5) * sigma
          if (trav(i)) then
             force_b = cmplx(ranf() - 0.5, ranf() - 0.5) * sigma
          else
             force_b = force_a
          end if
          
          if (trav(i) .and. abs(aimag(w_stir(i))) < epsilon(0.0)) then
             a_ant(i) = a_ant(i)*exp(-zi*w_stir(i)*dt)+force_a*dt
             b_ant(i) = 0.
          else
             a_ant(i) = a_ant(i)*exp(-zi*w_stir(i)*dt)+force_a*dt
             b_ant(i) = b_ant(i)*exp(zi*conjg(w_stir(i))*dt)+force_b*dt
          end if
       end if

! potential existed for a bug here, if different processors 
! got different random numbers; antennas driving different parts of
! velocity space at the same k could be driven differently.  BD 6.7.05

       call broadcast (a_ant)
       call broadcast (b_ant)

!       if (time < t0) then
!          apar(:,it,ik) = apar(:,it,ik) &
!               + (a_ant(i)+b_ant(i))/sqrt(2.)*exp(zi*kz_stir(i)*theta) &
!               * (0.5-0.5*cos(time*pi/t0)
!       else
! NOTE: Is the apar(:,it,ik) unnecessary on the RHS here (=0)?
! ANSWER: No, we are inside the nk_stir loop.  In general case, apar /= 0 here.
          apar(:,it,ik) = apar(:,it,ik) &
               + (a_ant(i)+b_ant(i))/sqrt(2.)*exp(zi*kz_stir(i)*theta) 
!       end if

       if (write_antenna) then
         if (proc0) write(out_unit, fmt='(8(1x,e13.6),1x,i2)') &
              & time, real((a_ant(i)+b_ant(i))/sqrt(2.)), &
              &     aimag((a_ant(i)+b_ant(i))/sqrt(2.)),real(a_ant(i)), aimag(a_ant(i)), &
              &     real(b_ant(i)), aimag(b_ant(i)), real(wtmp), i
       end if
    end do

    if (enforce_reality) then
       apar(:,1,1) = 0.0
       
       do it = 1, ntheta0/2
          apar(:,it+(ntheta0+1)/2,1) = conjg(apar(:,(ntheta0+1)/2+1-it,1))
       enddo
    end if

    !GGH NOTE: Here apar_new is equal to apar+ apar_ext
    apar_new = apar

  end subroutine antenna_amplitudes