diagnostics_antenna.fpp Source File

!> A module to calculate the properties of the
!> antenna used to drive turbulence
module diagnostics_antenna
  implicit none

  private

  public :: init_diagnostics_antenna, finish_diagnostics_antenna, write_jext
  public :: write_lorentzian

  real, dimension(:,:,:), allocatable ::  j_ext_hist

contains
  
  !> FIXME : Add documentation  
  subroutine init_diagnostics_antenna(gnostics)
    use diagnostics_config, only: diagnostics_type
    use kt_grids, only: ntheta0, naky
    use mp, only: proc0
    implicit none
    type(diagnostics_type), intent(in) :: gnostics
    if (gnostics%write_jext) then
       allocate (j_ext_hist(ntheta0, naky,0:gnostics%navg-1))
       j_ext_hist = 0.
    end if

    if (proc0.and.gnostics%ascii_files%write_to_jext) &
      write (gnostics%ascii_files%jext, *) &
      "Warning: the contents of this file differ &
      & to the contents generated by the old diagnostics module. In &
      & particular, the factor of 0.5 is no longer applied to then ky=0 &
      & mode."
  end subroutine init_diagnostics_antenna

  !> FIXME : Add documentation  
  subroutine finish_diagnostics_antenna(gnostics)
    use diagnostics_config, only: diagnostics_type
    implicit none
    type(diagnostics_type), intent(in) :: gnostics
    if (gnostics%write_jext.and.allocated(j_ext_hist)) deallocate(j_ext_hist)
  end subroutine finish_diagnostics_antenna

  !> FIXME : Add documentation  
  subroutine write_jext(gnostics)
    use kt_grids, only: ntheta0, naky
    use diagnostics_config, only: diagnostics_type
    use mp, only: proc0
    use warning_helpers, only: is_not_zero
#ifdef NETCDF
    use neasyf, only: neasyf_write
#endif
    use gs2_io, only: starts, kx_dim, ky_dim, time_dim
    implicit none
    type(diagnostics_type), intent(in) :: gnostics
    real:: t
    integer:: ik, it

    !GGH J_external
    real, dimension(:,:), allocatable ::  j_ext

    allocate (j_ext(ntheta0, naky)); j_ext=0.

    !Get z-centered j_ext at current time
    call calc_jext(gnostics,j_ext)

    !Write to netcdf
#ifdef NETCDF
    if (proc0) call neasyf_write(gnostics%file_id, "antenna_j_ext", j_ext, &
         dim_names=[kx_dim, ky_dim, time_dim], start=starts(3, gnostics%nout), &
         long_name="Time averaged external current in the antenna, real(kperp^2 A_antenna)", units="radians")
#endif
    !Write to ascii
    if (proc0 .and. gnostics%ascii_files%write_to_jext .and. (.not. gnostics%create)) then
       t = gnostics%user_time
       do ik = 1, naky
          do it = 1, ntheta0
             if (is_not_zero(j_ext(it, ik))) then
                write (unit=gnostics%ascii_files%jext, fmt="(es12.4,i4,i4,es12.4)")  &
                     t,it,ik,j_ext(it,ik)
             end if
          end do
       end do
    end if

    deallocate(j_ext)
  end subroutine write_jext

  !> FIXME : Add documentation  
  subroutine write_lorentzian(gnostics)
    use antenna, only: antenna_w
    use diagnostics_config, only: diagnostics_type
    use gs2_io, only: starts, netcdf_write_complex, complex_dim, time_dim
    use mp, only: proc0
    implicit none
    type(diagnostics_type), intent(in) :: gnostics
    if (.not. proc0) return
    call netcdf_write_complex(gnostics%file_id, "antenna_w", antenna_w(), &
         dim_names=[complex_dim, time_dim], start=starts(2, gnostics%nout))
  end subroutine write_lorentzian

  !> A subroutine to calculate the time-averaged antenna current
  !> j_ext = kperp^2 A_antenna.
  !>
  !> @note This looks like it may be reproducing the [[get_jext]] routine
  !> from [[antenna]]
  subroutine calc_jext (gnostics, j_ext)
    use mp, only: proc0
    use antenna, only: antenna_apar
    use volume_averages, only: average_theta
    use theta_grid, only: ntgrid
    use kt_grids, only: kperp2, ntheta0, naky
    use diagnostics_config, only: diagnostics_type
    implicit none
    !Passed
    type(diagnostics_type), intent(in) :: gnostics
    integer:: istep
    real, dimension(:,:), intent(out) ::  j_ext
    complex, dimension(:,:,:), allocatable :: j_extz
    !Local 
    integer :: i

    istep = gnostics%istep

    !call get_jext(j_ext)    
    allocate (j_extz(-ntgrid:ntgrid, ntheta0, naky)) ; j_extz = 0.
    call antenna_apar (kperp2, j_extz)       
    j_extz(ntgrid-1,:,:) = 0.
    call average_theta(real(j_extz), j_ext, gnostics%distributed)
    deallocate (j_extz)

    !Do averages with a history variable
    if (proc0) then
       !Save variable to history
       if (gnostics%navg > 1) then
          ! This looks a little odd as we ignore the first step
          if (istep > 1) j_ext_hist(:,:,mod(istep,gnostics%navg)) = j_ext(:,:)
          
          !Use average of history
          if (istep >= gnostics%navg) then
             j_ext=0.
             do i=0,gnostics%navg-1
                j_ext(:,:) = j_ext(:,:) + j_ext_hist(:,:,i) / real(gnostics%navg)
             end do
          end if
       end if
    end if
  end subroutine calc_jext
end module diagnostics_antenna