gca_update_payload Subroutine

private subroutine gca_update_payload(igrid, xpart, upart, qpart, mpart, mypayload, mynpayload, particle_time)

Uses

  • proc~~gca_update_payload~~UsesGraph proc~gca_update_payload gca_update_payload module~mod_global_parameters mod_global_parameters proc~gca_update_payload->module~mod_global_parameters module~mod_basic_types mod_basic_types module~mod_global_parameters->module~mod_basic_types module~mod_connectivity mod_connectivity module~mod_global_parameters->module~mod_connectivity module~mod_constants mod_constants module~mod_global_parameters->module~mod_constants module~mod_physicaldata mod_physicaldata module~mod_global_parameters->module~mod_physicaldata module~mod_variables mod_variables module~mod_global_parameters->module~mod_variables mpi mpi module~mod_global_parameters->mpi module~mod_variables->module~mod_basic_types

Update payload subroutine

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: igrid
double precision, intent(in) :: xpart(1:ndir)
double precision, intent(in) :: upart(1:ndir)
double precision, intent(in) :: qpart
double precision, intent(in) :: mpart
double precision, intent(out) :: mypayload(mynpayload)
integer, intent(in) :: mynpayload
double precision, intent(in) :: particle_time

Calls

proc~~gca_update_payload~~CallsGraph proc~gca_update_payload gca_update_payload proc~cross cross proc~gca_update_payload->proc~cross proc~get_bfield get_bfield proc~gca_update_payload->proc~get_bfield proc~get_efield get_efield proc~gca_update_payload->proc~get_efield proc~get_vec get_vec proc~gca_update_payload->proc~get_vec proc~mpistop mpistop proc~cross->proc~mpistop proc~get_bfield->proc~get_vec proc~interpolate_var interpolate_var proc~get_bfield->proc~interpolate_var proc~get_efield->proc~get_bfield proc~get_efield->proc~get_vec proc~get_efield->proc~interpolate_var proc~get_vec->proc~interpolate_var proc~interpolate_var->proc~mpistop mpi_abort mpi_abort proc~mpistop->mpi_abort

Called by

proc~~gca_update_payload~~CalledByGraph proc~gca_update_payload gca_update_payload proc~gca_create_particles gca_create_particles proc~gca_create_particles->proc~gca_update_payload proc~gca_integrate_particles gca_integrate_particles proc~gca_integrate_particles->proc~gca_update_payload proc~particles_create particles_create proc~particles_create->proc~gca_create_particles proc~main main proc~main->proc~particles_create program~agile agile program~agile->proc~main

Source Code

  subroutine gca_update_payload(igrid,xpart,upart,qpart,mpart,mypayload,&
     mynpayload,particle_time)
    use mod_global_parameters
    integer, intent(in)           :: igrid,mynpayload
    double precision, intent(in)  :: xpart(1:ndir),upart(1:ndir),qpart,mpart,&
       particle_time
    double precision, intent(out) :: mypayload(mynpayload)
    double precision, dimension(1:ndir) :: vE, e, b, bhat, vfluid, current
    double precision, dimension(1:ndir) :: drift1, drift2
    double precision, dimension(1:ndir) :: drift3, drift4, drift5, drift6,&
        drift7
    double precision, dimension(1:ndir) :: bdotgradb, vEdotgradb, gradkappaB
    double precision, dimension(1:ndir) :: bdotgradvE, vEdotgradvE
    double precision, dimension(1:ndir) :: gradBdrift, reldrift,&
        bdotgradbdrift
    double precision, dimension(1:ndir) :: vEdotgradbdrift, bdotgradvEdrift
    double precision, dimension(1:ndir) :: vEdotgradvEdrift
    double precision                    :: kappa, upar, absb, gamma, vpar,&
        vEabs
    double precision                    :: gradBdrift_abs, reldrift_abs, epar
    double precision                    :: bdotgradbdrift_abs,&
        vEdotgradbdrift_abs
    double precision                    :: bdotgradvEdrift_abs,&
        vEdotgradvEdrift_abs
    double precision                    :: momentumpar1, momentumpar2,&
        momentumpar3, momentumpar4

    call get_bfield(igrid, xpart(1:ndir), particle_time, b)
    call get_efield(igrid, xpart(1:ndir), particle_time, e)
!    call get_vec(bp, igrid,xpart(1:ndir),particle_time,b)
!    call get_vec(vp, igrid,xpart(1:ndir),particle_time,vfluid)
!    call get_vec(jp, igrid,xpart(1:ndir),particle_time,current)
!    e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) + particles_eta*current(1)
!    e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) + particles_eta*current(2)
!    e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) + particles_eta*current(3)

    absb         = sqrt(sum(b(:)**2))
    bhat(1:ndir) = b(1:ndir) / absb
    epar         = sum(e(:)*bhat(:))
    call cross(e,bhat,vE)
    vE(1:ndir)   = vE(1:ndir) / absb
    vEabs = sqrt(sum(vE(:)**2))
    if (relativistic) then
      kappa = 1.d0/sqrt(1.0d0 - sum(vE(:)**2)/c_norm**2)
    else
      kappa = 1.d0
    end if
    vpar = upart(1)/upart(3)
    upar = upart(1)

    call get_vec(b_dot_grad_b, igrid,xpart(1:ndir),particle_time,bdotgradb)
    call get_vec(vE_dot_grad_b, igrid,xpart(1:ndir),particle_time,vEdotgradb)
    call get_vec(grad_kappa_B, igrid,xpart(1:ndir),particle_time,gradkappaB)
    call get_vec(b_dot_grad_vE, igrid,xpart(1:ndir),particle_time,bdotgradvE)
    call get_vec(vE_dot_grad_vE, igrid,xpart(1:ndir),particle_time,&
       vEdotgradvE)

    drift1(1:ndir) = bhat(1:ndir)/(absb/kappa**2)
    drift2(1:ndir) = upart(2)/(upart(3)*qpart)*gradkappaB(1:ndir)

    call cross(drift1,drift2,gradBdrift)
    gradBdrift_abs = sqrt(sum(gradBdrift(:)**2))

    drift3(1:ndir) = upar*epar/upart(3)*vE(1:ndir)
    call cross(drift1,drift3,reldrift)
    reldrift_abs = sqrt(sum(reldrift(:)**2))

    drift4(1:ndir) = mpart/qpart* ( upar**2/upart(3)*bdotgradb(1:ndir))
    call cross(drift1,drift4,bdotgradbdrift)
    bdotgradbdrift_abs = sqrt(sum(bdotgradbdrift(:)**2))

    drift5(1:ndir) = mpart/qpart* ( upar*vEdotgradb(1:ndir))
    call cross(drift1,drift5,vEdotgradbdrift)
    vEdotgradbdrift_abs = sqrt(sum(vEdotgradbdrift(:)**2))

    drift6(1:ndir) = mpart/qpart* ( upar*bdotgradvE(1:ndir))
    call cross(drift1,drift6,bdotgradvEdrift)
    bdotgradvEdrift_abs = sqrt(sum(bdotgradvEdrift(:)**2))

    drift7(1:ndir) = mpart/qpart* (upart(3)*vEdotgradvE(1:ndir))
    call cross(drift1,drift7,vEdotgradvEdrift)
    vEdotgradvEdrift_abs = sqrt(sum(vEdotgradvEdrift(:)**2))

    momentumpar1 = qpart/mpart*epar
    momentumpar2 = -(upart(2)/mpart/upart(3))*sum(bhat(:)*gradkappaB(:))
    momentumpar3 = upar*sum(vE(:)*bdotgradb(:))
    momentumpar4 = upart(3)*sum(vE(:)*vEdotgradb(:))

    ! Payload update
    if (mynpayload > 0) then
      ! current gyroradius
      mypayload(1) = sqrt(2.0d0*mpart*upart(2)*absb)/abs(qpart*absb)
    end if
    if (mynpayload > 1) then
      ! pitch angle
      mypayload(2) = atan2(sqrt((2.0d0*upart(2)*absb)/(mpart*upart(3)**2)),&
         vpar)
    end if
    if (mynpayload > 2) then
      ! particle v_perp
      mypayload(3) = sqrt((2.0d0*upart(2)*absb)/(mpart*upart(3)**2))
    end if
    if (mynpayload > 3) then
      ! particle parallel momentum term 1
      mypayload(4) = momentumpar1
    end if
    if (mynpayload > 4) then
      ! particle parallel momentum term 2
      mypayload(5) = momentumpar2
    end if
    if (mynpayload > 5) then
      ! particle parallel momentum term 3
      mypayload(6) = momentumpar3
    end if
    if (mynpayload > 6) then
      ! particle parallel momentum term 4
      mypayload(7) = momentumpar4
    end if
    if (mynpayload > 7) then
      ! particle ExB drift
      mypayload(8) = vEabs
    end if
    if (mynpayload > 8) then
      ! relativistic drift
      mypayload(9) = gradBdrift_abs
    end if
    if (mynpayload > 9) then
      ! gradB drift
      mypayload(10) = reldrift_abs
    end if
    if (mynpayload > 10) then
      ! bdotgradb drift
      mypayload(11) = bdotgradbdrift_abs
    end if
    if (mynpayload > 11) then
      ! vEdotgradb drift
      mypayload(12) = vEdotgradbdrift_abs
    end if
    if (mynpayload > 12) then
      ! bdotgradvE drift
      mypayload(13) = bdotgradvEdrift_abs
    end if
    if (mynpayload > 13) then
      ! vEdotgradvE drift
      mypayload(14) = vEdotgradvEdrift_abs
    end if

  end subroutine gca_update_payload