get_source_term Subroutine

  subroutine get_source_term &
       (phi, apar, bpar, phinew, aparnew, bparnew, include_explicit, &
        ab_coefficients, isgn, iglo,ik,it,il,ie,is, sourcefac, source)
    use dist_fn_arrays, only: aj0, aj1, vperp2, vpac, g, wstar
    use dist_fn_arrays, only: gexp_1, gexp_2, gexp_3
    use theta_grid, only: ntgrid
    use kt_grids, only: aky
    use le_grids, only: lmax, grid_has_trapped_particles, is_ttp, can_be_ttp
    use species, only: spec, nonmaxw_corr
    use run_parameters, only: fphi, fapar, fbpar
    use gs2_time, only: code_dt, wunits
    use hyper, only: D_res
    use constants, only: zi
    implicit none
    complex, dimension (-ntgrid:,:,:), intent (in) :: phi,    apar,    bpar
    complex, dimension (-ntgrid:,:,:), intent (in) :: phinew, aparnew, bparnew
    logical, intent (in) :: include_explicit
    real, dimension(:), allocatable, intent(in) :: ab_coefficients
    integer, intent (in) :: isgn, iglo, ik, it, il, ie, is
    complex, intent (in) :: sourcefac
    complex, dimension (-ntgrid:), intent (out) :: source

    integer :: ig
    complex, dimension (-ntgrid:ntgrid) :: phigavg, apargavg

!CMR, 4/8/2011
! apargavg and phigavg combine to give the GK EM potential chi. 
!          chi = phigavg - apargavg*vpa(:,isgn,iglo)*spec(is)%stm
! phigavg  = phi J0 + 2 T_s/q_s . vperp^2 bpar/bmag J1/Z
! apargavg = apar J0 
! Both quantities are decentred in time and evaluated on || grid points
!
    phigavg  = (fexp(is)*phi(:,it,ik)   + (1.0-fexp(is))*phinew(:,it,ik)) &
                *aj0(:,iglo)*fphi &
             + (fexp(is)*bpar(:,it,ik) + (1.0-fexp(is))*bparnew(:,it,ik))&
                *aj1(:,iglo)*fbpar*2.0*vperp2(:,iglo)*spec(is)%tz
    apargavg = (fexp(is)*apar(:,it,ik)  + (1.0-fexp(is))*aparnew(:,it,ik)) &
                *aj0(:,iglo)*fapar

! source term in finite difference equations
    select case (source_option_switch)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Default choice: solve self-consistent equations
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    case (source_option_full)
       if (il <= lmax) then
          call set_source
       else
          source = 0.0
       end if       

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Solve self-consistent terms + include external i omega_d * F_0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    case(source_option_phiext_full)
       if (il <= lmax) then
          call set_source
          if (aky(ik) < epsilon(0.0)) then
             source = source &
                  - zi * wdrift(:ntgrid-1, isgn, iglo) * nonmaxw_corr(ie, is) &
                  * 2.0 * phi_ext * sourcefac * aj0(:ntgrid-1, iglo)
          end if
       else
          source = 0.0
       end if
    case(source_option_homogeneous)
       source = 0.0
    end select

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Do matrix multiplications
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    if (il <= lmax) then
       if (isgn == 1) then
          do ig = -ntgrid, ntgrid-1
             source(ig) = source(ig) &
                  + b(ig,1,iglo)*g(ig,1,iglo) + a(ig,1,iglo)*g(ig+1,1,iglo)
          end do
       else
          do ig = -ntgrid, ntgrid-1
             source(ig) = source(ig) &
                  + a(ig,2,iglo)*g(ig,2,iglo) + b(ig,2,iglo)*g(ig+1,2,iglo)
          end do
       end if
    end if

!CMR, 21/7/2014: removed redundant line here:    source(ntgrid)=source(-ntgrid) 
!                as source(ntgrid) should never be used.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! special source term for totally trapped particles
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!CMR, 13/10/2014: 
! Upper limit of following loops setting source changed from ntgrid to ntgrid-1
! Source is allocated as: source(-ntgrid:ntgrid-1), so ntgrid is out of bounds.

    if (source_option_switch == source_option_full .or. &
        source_option_switch == source_option_phiext_full) then
       if (grid_has_trapped_particles() .and. isgn == 2 .and. can_be_ttp(il)) then
          do ig = -ntgrid, ntgrid-1
             if (.not. is_ttp(ig, il)) cycle
             source(ig) &
                  ! Note we could replace the explicit "2" appearing below with isgn for consistency.
                  = g(ig,2,iglo)*a(ig,2,iglo) &
                  - zi*(wdriftttp(ig, 2, iglo))*nonmaxw_corr(ie,is)*phigavg(ig) &
                  + zi*wstar(ik,ie,is)*phigavg(ig)
          end do

          if (source_option_switch == source_option_phiext_full .and. &
               aky(ik) < epsilon(0.0)) then
             do ig = -ntgrid, ntgrid-1
                if (.not. is_ttp(ig, il)) cycle
                source(ig) = source(ig) - zi* &
                     wdriftttp(ig, isgn, iglo)*nonmaxw_corr(ie,is)*2.0*phi_ext*sourcefac*aj0(ig,iglo)
             end do
          endif

          if (include_explicit) then
             select case (size(ab_coefficients))
             case (1)
                do ig = -ntgrid, ntgrid-1
                   if (.not. is_ttp(ig, il)) cycle
                   source(ig) = source(ig) + 0.5*code_dt*(&
                        ab_coefficients(1)*gexp_1(ig,isgn,iglo))
                end do
             case (2)
                do ig = -ntgrid, ntgrid-1
                   if (.not. is_ttp(ig, il)) cycle
                   source(ig) = source(ig) + 0.5*code_dt*( &
                        ab_coefficients(1)*gexp_1(ig,isgn,iglo) + &
                        ab_coefficients(2)*gexp_2(ig,isgn,iglo))
                end do
             case (3)
                do ig = -ntgrid, ntgrid-1
                   if (.not. is_ttp(ig, il)) cycle
                   source(ig) = source(ig) + 0.5*code_dt*( &
                        ab_coefficients(1)*gexp_1(ig,isgn,iglo) + &
                        ab_coefficients(2)*gexp_2(ig,isgn,iglo) + &
                        ab_coefficients(3)*gexp_3(ig,isgn,iglo))
                end do
             end select
          end if
       end if
    end if

  contains
    !> FIXME : Add documentation
    subroutine set_source

      use species, only: spec, nonmaxw_corr
      use theta_grid, only: itor_over_B
      implicit none
      complex :: apar_p, apar_m, phi_p, phi_m
      real :: bd, bdfac_p, bdfac_m
      bd = bkdiff(is)
      bdfac_p = 1.+bd*(3.-2.*real(isgn))
      bdfac_m = 1.-bd*(3.-2.*real(isgn))

!CMR, 4/8/2011:
! Some concerns, may be red herrings !
! (1) no bakdif factors in phi_m, apar_p, apar_m, vpar !!! 
!                        (RN also spotted this for apar_p)
! (2) can interpolations of products be improved? 
!
!  Attempt at variable documentation:
! phigavg  = phi J0 + 2 T_s/q_s . vperp^2 bpar/bmag J1/Z
! apargavg = apar J0                        (decentered in t) 
! NB apargavg and phigavg combine to give the GK EM potential chi
! chi = phigavg - apargavg*vpa(:,isgn,iglo)*spec(is)%stm
! phi_p = 2 phigavg                      .... (roughly!)
! phi_m = d/dtheta (phigavg)*DTHETA 
! apar_p = 2 apargavg  
! apar_m = 2 d/dt (apar)*DELT  (gets multiplied later by J0 and vpa when included in source)
! => phi_p - apar_p*vpa(:,isgn,iglo)*spec(is)%stm = 2 chi  .... (roughly!)  
! vparterm = -2.0*vpar (IN ABSENCE OF LOWFLOW TERMS)
! wdfac = wdrift + wcoriolis/spec(is)%stm (IN ABSENCE OF LOWFLOW TERMS)
! wstarfac = wstar  (IN ABSENCE OF LOWFLOW TERMS)
! Line below is not true for alphas. vpar = q_s/Tstar * v_|| * the rest. EGH/GW
! vpar = q_s/sqrt{T_s m_s} (v_||^GS2). \gradpar(theta)/DTHETA . DELT (centred)
! wdrift =    q_s/T_s  v_d.\grad_perp . DELT 
! wcoriolis = q_s/T_s  v_C.\grad_perp . DELT 
!
! Definition of source:= 2*code_dt*RHS of GKE
! source     appears to contain following physical terms
!   -2q_s/T_s v||.grad(J0 phi + 2 vperp^2 bpar/bmag J1/Z T_s/q_s).delt 
!   -2d/dt(q v|| J0 apar / T).delt
!   +hyperviscosity
!   -2 v_d.\grad_perp (q J0 phi/T + 2 vperp^2 bpar/bmag J1/Z).delt 
!   -coriolis terms
!   2{\chi,f_{0s}}.delt  (allowing for sheared flow)
!CMRend
!GJW (2018): 
! Ensuring apar respects bakdif. This fixes a numerical instability that prevents 
! nonlinear electromagnetic runs.

      do ig = -ntgrid, ntgrid-1
         phi_p = bdfac_p*phigavg(ig+1)+bdfac_m*phigavg(ig)
         phi_m = phigavg(ig+1)-phigavg(ig)
         ! RN> bdfac factors seem missing for apar_p
         apar_p = bdfac_p*apargavg(ig+1)+bdfac_m*apargavg(ig)
         apar_m = bdfac_p * (aparnew(ig+1,it,ik) - apar(ig+1,it,ik)) * fapar + &
              bdfac_m * (aparnew(ig,it,ik) - apar(ig,it,ik)) * fapar

!MAB, 6/5/2009:
! added the omprimfac source term arising with equilibrium flow shear  
!CMR, 26/11/2010:
! Useful to understand that all source terms propto aky are specified here 
! using Tref=mref vtref^2. See 
! [note by BD and MK on "Microinstabilities in Axisymmetric Configurations"].
! This is converted to  the standard internal gs2 normalisation, 
! Tref=(1/2) mref vtref^2, by wunits, which contains a crucial factor 1/2.
! (Would be less confusing if always used same Tref!)
!
         source(ig) = (-2.0*vpar(ig,isgn,iglo)*nonmaxw_corr(ie,is)*phi_m &
              -spec(is)%zstm*vpac(ig,isgn,iglo)*nonmaxw_corr(ie,is) &
              *((aj0(ig+1,iglo) + aj0(ig,iglo))*0.5*apar_m  &
              + D_res(it,ik)*apar_p) &
              -zi*wdrift(ig,isgn,iglo)*nonmaxw_corr(ie,is)*phi_p) &
              + zi*(wstar(ik,ie,is) &
              + vpac(ig,isgn,iglo)*code_dt*wunits(ik)*ufac(ie,is) &
              -2.0*omprimfac*vpac(ig,isgn,iglo)*nonmaxw_corr(ie,is)*code_dt*wunits(ik)*g_exb*itor_over_B(ig)/spec(is)%stm) &
              *(phi_p - apar_p*spec(is)%stm*vpac(ig,isgn,iglo))
      end do

      if (include_explicit) then
         select case (size(ab_coefficients))
         case (1)
            do ig = -ntgrid, ntgrid-1
               source(ig) = source(ig) + 0.5*code_dt*(&
                    ab_coefficients(1)*gexp_1(ig,isgn,iglo))
            end do
         case (2)
            do ig = -ntgrid, ntgrid-1
               source(ig) = source(ig) + 0.5*code_dt*( &
                    ab_coefficients(1)*gexp_1(ig,isgn,iglo) + &
                    ab_coefficients(2)*gexp_2(ig,isgn,iglo))
            end do
         case (3)
            do ig = -ntgrid, ntgrid-1
               source(ig) = source(ig) + 0.5*code_dt*( &
                    ab_coefficients(1)*gexp_1(ig,isgn,iglo) + &
                    ab_coefficients(2)*gexp_2(ig,isgn,iglo) + &
                    ab_coefficients(3)*gexp_3(ig,isgn,iglo))
            end do
         end select
      end if
    end subroutine set_source

  end subroutine get_source_term