vspace_derivatives Subroutine

  subroutine vspace_derivatives (g, gold, g1, phi, bpar, phinew, bparnew, diagnostics, gtoc, ctog)
    use array_utils, only: copy, zero_array
    use redistribute, only: gather, scatter
    use dist_fn_arrays, only: g_adjust, to_g_gs2, from_g_gs2
    use gs2_layouts, only: g_lo, le_lo
    use gs2_time, only: code_dt
    use theta_grid, only: ntgrid
    use le_grids, only: integrate_moment, negrid, g2le
    use collisions, only: solfp1, colls, adjust, heating, hyper_colls, use_le_layout, nxi_lim, c_rate
    use job_manage, only: time_message
    use mp, only: get_mp_times
    use optionals, only: get_option_with_default
    implicit none

    complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent (in out) :: g, gold, g1
    complex, dimension (-ntgrid:,:,:), intent (in) :: phi, bpar, phinew, bparnew
    integer, optional, intent (in) :: diagnostics

    !CMR, 12/9/2013:
    !CMR   New logical optional input parameters gtoc, ctog used to set
    !CMR   flags (g_to_c and c_to_g) to control whether redistributes required
    !CMR   to map g_lo to collision_lo, and collision_lo to g_lo.
    !CMR   All redistributes are performed by default.
    !CMR

    !CMR, 3/10/2013:
    !CMR   Just realised we need to be careful with g_adjust if we avoid
    !CMR   mapping g_lo <=> le_lo    Mmmm, a little more to think about here ;-(
    !CMR
    logical, intent(in), optional :: gtoc, ctog
    logical :: g_to_c, c_to_g
    complex, dimension (:,:,:), allocatable :: gle
    complex, dimension (:,:,:), allocatable :: gc1, gc2, gc3
    logical :: heating_flag
    real :: mp_total, mp_total_after

    call time_message(.false., time_vspace_derivatives, "Vspace derivative")
    call get_mp_times(total_time = mp_total)
    g_to_c = get_option_with_default(gtoc, .true.)
    c_to_g = get_option_with_default(ctog, .true.)

    heating_flag = heating .and. present(diagnostics)

    if(use_le_layout) then

       if (colls) then
          if (adjust) then
             call g_adjust (g, phinew, bparnew, direction = from_g_gs2)
             if (heating_flag) call g_adjust (gold, phi, bpar, direction = from_g_gs2)
          end if
          if (heating_flag) then
             allocate (gc3(-ntgrid:ntgrid,2,g_lo%llim_proc:g_lo%ulim_alloc))
             call copy(g, gc3)
          end if

          allocate (gle(nxi_lim,negrid+1,le_lo%llim_proc:le_lo%ulim_alloc))
          call zero_array(gle)

          ! map data from g_layout to le_layout
          if (g_to_c) call gather (g2le, g, gle)

          if (colls) then
             ! update distribution function to take into account collisions
             call solfp1 (gle, diagnostics)
          end if

          ! remap from le_layout to g_layout
          if (c_to_g) call scatter (g2le, gle, g)
          deallocate (gle)
          if (heating_flag) then
             ! form (h_i+1 + h_i)/2 * C(h_i+1) and integrate.
             gc3 = 0.5*conjg(g+gold)*(g-gc3)/code_dt
             ! We only want the real part of this result so we could
             ! either provide a real argument overload or reduce to
             ! a complex temporary before storing just the real part.
             ! This would allow c_rate to be declared real, halving the
             ! memory requirements.
             call integrate_moment (gc3, c_rate(:,:,:,:,3))
             deallocate (gc3)
          endif

          if (adjust) then
             call g_adjust (g, phinew, bparnew,  direction = to_g_gs2)
             if (heating_flag) call g_adjust (gold, phi, bpar, direction = to_g_gs2)
          end if
       end if
    else
       if (colls) then
          if (heating_flag) then
             allocate (gc1(-ntgrid:ntgrid,2,g_lo%llim_proc:g_lo%ulim_alloc))
             if (hyper_colls) then
                allocate (gc2(-ntgrid:ntgrid,2,g_lo%llim_proc:g_lo%ulim_alloc))
             else
                allocate (gc2(1,1,1))
             end if
             allocate (gc3(-ntgrid:ntgrid,2,g_lo%llim_proc:g_lo%ulim_alloc))
          else
             allocate (gc1(1,1,1))
             allocate (gc2(1,1,1))
             allocate (gc3(1,1,1))
          end if
          call zero_array(gc1)
          call zero_array(gc2)
          call zero_array(gc3)

          if (adjust) then
             call g_adjust (g, phinew, bparnew, direction = from_g_gs2)
             if (heating_flag) call g_adjust (gold, phi, bpar, direction = from_g_gs2)
          end if

          if (heating_flag) call copy(g, gc3)

          ! update distribution function to take into account collisions
          call solfp1 (g, g1, gc1, gc2, diagnostics)

          if (heating_flag) then
             call integrate_moment (gc1, c_rate(:,:,:,:,1))
             deallocate (gc1)

             ! We only want the real part of this result so we could
             ! either provide a real argument overload or reduce to
             ! a complex temporary before storing just the real part.
             ! This would allow c_rate to be declared real, halving the
             ! memory requirements.
             if (hyper_colls) call integrate_moment (gc2, c_rate(:,:,:,:,2))
             deallocate (gc2)

             ! form (h_i+1 + h_i)/2 * C(h_i+1) and integrate.
             gc3 = 0.5*conjg(g+gold)*(g-gc3)/code_dt

             call integrate_moment (gc3, c_rate(:,:,:,:,3))

             deallocate (gc3)
          end if

          if (adjust) then
             call g_adjust (g, phinew, bparnew, direction = to_g_gs2)
             if (heating_flag) call g_adjust (gold, phi, bpar, direction = to_g_gs2)
          end if
       end if

    end if

    call time_message(.false., time_vspace_derivatives, "Vspace derivative")
    call get_mp_times(total_time = mp_total_after)
    time_vspace_derivatives_mpi = time_vspace_derivatives_mpi + (mp_total_after - mp_total)

  end subroutine vspace_derivatives