getfieldeq1 Subroutine

private subroutine getfieldeq1(phi, apar, bpar, antot, antota, antotp, fieldeq, fieldeqa, fieldeqp, local_only)

FIXME : Add documentation

Arguments

Type IntentOptional Attributes Name
complex, intent(in), dimension (-ntgrid:,:,:) :: phi
complex, intent(in), dimension (-ntgrid:,:,:) :: apar
complex, intent(in), dimension (-ntgrid:,:,:) :: bpar
complex, intent(in), dimension (-ntgrid:,:,:) :: antot
complex, intent(in), dimension (-ntgrid:,:,:) :: antota
complex, intent(in), dimension (-ntgrid:,:,:) :: antotp
complex, intent(inout), dimension (-ntgrid:,:,:) :: fieldeq
complex, intent(inout), dimension (-ntgrid:,:,:) :: fieldeqa
complex, intent(inout), dimension (-ntgrid:,:,:) :: fieldeqp
logical, intent(in), optional :: local_only

Contents

Source Code


Source Code

  subroutine getfieldeq1 (phi, apar, bpar, antot, antota, antotp, &
       fieldeq, fieldeqa, fieldeqp, local_only)
    use theta_grid, only: ntgrid, bmag
    use kt_grids, only: naky, ntheta0, kperp2
    use run_parameters, only: has_phi, has_apar, has_bpar, beta
    use species, only: spec, has_electron_species
    use optionals, only: get_option_with_default
    use gs2_layouts, only: g_lo
    implicit none
    complex, dimension (-ntgrid:,:,:), intent (in) :: phi, apar, bpar
    complex, dimension (-ntgrid:,:,:), intent (in) :: antot, antota, antotp
    complex, dimension (-ntgrid:,:,:), intent (in out) ::fieldeq,fieldeqa,fieldeqp
    logical, intent(in), optional :: local_only
    logical :: local_local_only
    integer :: ik, it
    integer :: it_llim, it_ulim, ik_llim, ik_ulim

    local_local_only = get_option_with_default(local_only, .false.)

    if (local_local_only) then
       it_llim = g_lo%it_min ; it_ulim = g_lo%it_max
       ik_llim = g_lo%ik_min ; ik_ulim = g_lo%ik_max
    else
       it_llim = 1 ; it_ulim = ntheta0
       ik_llim = 1 ; ik_ulim = naky
    end if

    if (.not. allocated(fl_avg)) allocate (fl_avg(ntheta0, naky))
    fl_avg = 0.

    if (.not. has_electron_species(spec)) then
       if (adiabatic_option_switch == adiabatic_option_fieldlineavg) then
          call calculate_flux_surface_average(fl_avg,antot)
       end if
    end if

    if (has_phi) then
       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(ik, it) &
       !$OMP SHARED(ik_llim, ik_ulim, it_llim, it_ulim, fieldeq, antot, &
       !$OMP bpar, gamtot1, gamtot, gridfac1, phi, fl_avg) &
       !$OMP COLLAPSE(2) &
       !$OMP SCHEDULE(static)
       do ik = ik_llim, ik_ulim
          do it = it_llim, it_ulim
             fieldeq(:,it,ik) = antot(:,it,ik) + &
                  bpar(:,it,ik)*gamtot1(:,it,ik) - &
                  gamtot(:,it,ik)*gridfac1(:,it,ik)*phi(:,it,ik) + &
                  fl_avg(it,ik)
          end do
       end do
       !$OMP END PARALLEL DO
    end if

    if (has_apar) then
       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(ik, it) &
       !$OMP SHARED(ik_llim, ik_ulim, it_llim, it_ulim, fieldeqa, antota, &
       !$OMP kperp2, gridfac1, apar) &
       !$OMP COLLAPSE(2) &
       !$OMP SCHEDULE(static)
       do ik = ik_llim, ik_ulim
          do it = it_llim, it_ulim
             fieldeqa(:, it, ik) = antota(:, it, ik) - &
                  kperp2(:, it, ik)*gridfac1(:, it, ik)*apar(:, it, ik)
          end do
       end do
       !$OMP END PARALLEL DO
    end if

    ! bpar == delta B_parallel / B_0(theta) b/c of the factor of 1/bmag(theta)**2
    ! in the following
    if (has_bpar) then
       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(ik, it) &
       !$OMP SHARED(ik_llim, ik_ulim, it_llim, it_ulim, fieldeqp, antotp, bpar, &
       !$OMP gamtot2, phi, gamtot1, beta, apfac, bmag, gridfac1)  &
       !$OMP COLLAPSE(2) &
       !$OMP SCHEDULE(static)
       do ik = ik_llim, ik_ulim
          do it = it_llim, it_ulim
             fieldeqp(:,it,ik) = (antotp(:, it, ik) + &
                  bpar(:, it, ik)*gamtot2(:, it, ik) + &
                  0.5*phi(:, it, ik)*gamtot1(:, it, ik))*(beta*apfac/bmag**2) + &
                  bpar(:, it, ik)*gridfac1(:, it, ik)
          end do
       end do
       !$OMP END PARALLEL DO
    end if
  end subroutine getfieldeq1