! $Id$
!
!  This module integrates drag forces between particles and gas.
!
!  Reference:
!    Yang, C.-C., & Johansen, A. 2016, ApJS, in press (arXiv:1603.08523)
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
!
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: lparticles_drag = .true.
!
!***************************************************************
module Particles_drag
!
  use Cdata
  use Cparam
  use Messages
  use Particles_cdata
  use Particles_map
!
  implicit none
!
  include 'particles_drag.h'
!
!  Initialization parameters
!
  logical :: ldrag_equilibrium_global = .false.  ! Set global/local drag equilibrium when lset_drag_equilibrium = .true.
  logical :: ldrag_on_gas = .false.              ! Turn on/off drag on gas
  logical :: ldrag_on_par = .false.              ! Turn on/off drag on particles
  logical :: lset_drag_equilibrium = .false.     ! Initialize the velocities of the particles and the gas under equilibrium
  real :: gx_gas = 0.0        ! Background acceleration of gas in x direction
  real :: gz_par_coeff = 0.0  ! Coefficient for background acceleration of a particle in z direction
  real :: taus = 0.0          ! Dimensionless stopping time: Omega * tdrag
  real :: tdrag = 0.0         ! Drag timescale
!
  namelist /particles_drag_init_pars/ &
      ldrag_equilibrium_global, ldrag_on_gas, ldrag_on_par, lset_drag_equilibrium, gx_gas, gz_par_coeff, taus, tdrag
!
!  Runtime parameters
!
  logical :: ldrag_pm_back_reaction = .true.  ! Couple the gas by particle-meshing back reaction from particles
!
  namelist /particles_drag_run_pars/ ldrag_on_gas, ldrag_on_par, ldrag_pm_back_reaction, gx_gas, gz_par_coeff, taus, tdrag
!
!  Module variables
!
  real :: dv_gas = 0.0
  real :: epicycle_freq = 0.0
  real :: epicycle_ratio = 0.0
  real :: gzpc = 0.0
  real :: oneplustaus2inv = 0.0
  real :: taus2 = 0.0
  real :: twominusqtaus = 0.0
  real :: twotaus = 0.0
  real :: twoomega1 = 0.0
!
  contains
!***********************************************************************
    subroutine register_particles_drag()
!
!  Register this module.
!
!  14-dec-14/ccyang: coded.
!
      if (lroot) call svn_id("$Id$")
!
    endsubroutine register_particles_drag
!***********************************************************************
    subroutine initialize_particles_drag()
!
!  Perform any post-parameter-read initialization, i.e., calculate
!  derived parameters.
!
!  11-dec-15/ccyang: coded.
!
!  Sanity checks.
!
      if (coord_system /= 'cartesian' .and. Omega /= 0.0) &
          call fatal_error('initialize_particles_drag', 'non-inertial curvilinear system is not implemented. ')
      if (ldrag_on_gas .and. iuu == 0) call fatal_error('initialize_particles_drag', 'drag on gas requires uu. ')
!
!  Check the drag timescale.
!
      drag: if (ldrag_on_par) then
        settdrag: if (tdrag == 0.0) then
          if (taus == 0.0) call fatal_error('initialize_particles_drag', 'tdrag = 0')
          if (Omega == 0.0) call fatal_error('initialize_particles_drag', 'Omega = 0')
          tdrag = taus / Omega
        elseif (taus == 0.0) then settdrag
          taus = Omega * tdrag
        endif settdrag
      endif drag
      taus2 = 2.0 * (2.0 - qshear) * taus**2
      twotaus = 2.0 * taus
      twominusqtaus = (2.0 - qshear) * taus
      oneplustaus2inv = 1.0 / (1.0 + taus2)
      if (lroot) print *, 'initialize_particles_drag: tdrag = ', tdrag
!
!  Find the frequency and aspect ratio of the epicycles.
!
      epicycle_freq = sqrt(2.0 * (2.0 - qshear)) * Omega
      epicycle_ratio = sqrt(2.0 / (2.0 - qshear))
!
!  Get the velocity reduction of the gas.
!
      dv: if (gx_gas /= 0.0) then
        if (Omega == 0.0) call fatal_error('initialize_particles_drag', 'Omega = 0')
        if (.not. ldrag_on_gas .and. .not. ldrag_on_par) call fatal_error('initialize_particles_drag', 'gx_gas has no effect. ')
        twoomega1 = 0.5 / Omega
        dv_gas = twoomega1 * gx_gas
      endif dv
!
!  Initialization for particle_zaccel().
!
      gzpc = gz_par_coeff**2
!
!  Calculate and overwrite mp_swarm and rhop_swarm.
!
      mprhop: if (eps_dtog > 0.0) then
        mp_swarm = find_mp_swarm(eps_dtog)
        rhop_swarm = mp_swarm / product((/dx,dy,dz/), lactive_dimension)
        info: if (lroot) then
          print *, 'initialize_particles_drag: reset mp_swarm = ', mp_swarm
          print *, 'initialize_particles_drag: reset rhop_swarm = ', rhop_swarm
        endif info
      endif mprhop
!
    endsubroutine initialize_particles_drag
!***********************************************************************
    subroutine init_particles_drag(f, fp)
!
!  Set some initial conditions for the gas and the particles, if any.
!
!  08-may-16/ccyang: coded.
!
      real, dimension(mx,my,mz,mfarray), intent(inout) :: f
      real, dimension(mpar_loc,mparray), intent(inout) :: fp
!
      type(pic), dimension(nx,ny,nz) :: cell
      type(particle), dimension(:), pointer :: ghost
      integer, dimension(:), pointer :: sendlist
      integer, dimension(0:nproc_comm) :: ngp_send, ngp_recv
      integer :: npsend
!
      real :: ux, uy, vx, vy
!
!  Initialize the gas and the particle velocities with drag equilibrium, if requested.
!
      init: if (lset_drag_equilibrium) then
        global: if (ldrag_equilibrium_global) then
          if (lroot) print *, 'init_particles_drag: Set global drag equilibrium. '
          call get_nsh_solution(dv_gas, eps_dtog, ux, uy, vx, vy)
          f(l1:l2,m1:m2,n1:n2,iux) = f(l1:l2,m1:m2,n1:n2,iux) + ux
          f(l1:l2,m1:m2,n1:n2,iuy) = f(l1:l2,m1:m2,n1:n2,iuy) + uy
          fp(1:npar_loc,ivpx) = fp(1:npar_loc,ivpx) + vx
          fp(1:npar_loc,ivpy) = fp(1:npar_loc,ivpy) + vy
        else global
          if (lroot) print *, 'init_particles_drag: Set local drag equilibrium. '
          if (lparticles_blocks) call fatal_error('init_particles_drag', 'particles_blocks version is not implemented yet. ')
          nullify(ghost)
          call distribute_particles(f, fp, npsend, sendlist, ghost, cell, ngp_send, ngp_recv)
          call set_drag_equilibrium(cell)
          call collect_particles(f, fp, npsend, sendlist, ghost, cell, ngp_send, ngp_recv, &
                                 ldrag_on_par, ldrag_on_gas .or. Omega /= 0.0, .false.)
        endif global
      endif init
!
    endsubroutine init_particles_drag
!***********************************************************************
    subroutine read_particles_drag_init_pars(iostat)
!
!  Read initialization parameters from namelist particles_drag_init_pars.
!
!  27-aug-15/ccyang: coded.
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=particles_drag_init_pars, IOSTAT=iostat)
!
    endsubroutine read_particles_drag_init_pars
!***********************************************************************
    subroutine write_particles_drag_init_pars(unit)
!
!  Write initialization parameters to namelist particles_drag_init_pars.
!
!  14-dec-15/ccyang: coded.
!
      integer, intent(in) :: unit
!
      integer :: stat
!
      write(unit, NML=particles_drag_init_pars, IOSTAT=stat)
      if (stat /= 0) call fatal_error('write_particles_drag_init_pars', 'cannot write particles_drag_init_pars. ', force=.true.)
!
    endsubroutine write_particles_drag_init_pars
!***********************************************************************
    subroutine read_particles_drag_run_pars(iostat)
!
!  Read runtime parameters from namelist particles_drag_run_pars.
!
!  27-aug-15/ccyang: coded.
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=particles_drag_run_pars, IOSTAT=iostat)
!
    endsubroutine read_particles_drag_run_pars
!***********************************************************************
    subroutine write_particles_drag_run_pars(unit)
!
!  Write runtime parameters to namelist particles_drag_run_pars.
!
!  14-dec-15/ccyang: coded.
!
      integer, intent(in) :: unit
!
      integer :: stat
!
      write(unit, NML=particles_drag_run_pars, IOSTAT=stat)
      if (stat /= 0) call fatal_error('write_particles_drag_run_pars', 'cannot write particles_drag_run_pars. ', force=.true.)
!
    endsubroutine write_particles_drag_run_pars
!***********************************************************************
    subroutine integrate_drag(f, fp, dt)
!
!  Wrapper for the integration of the drag force between particles and
!  gas.
!
!  08-may-16/ccyang: coded.
!
      real, dimension(mx,my,mz,mfarray), intent(inout) :: f
      real, dimension(mpar_loc,mparray), intent(inout) :: fp
      real, intent(in) :: dt
!
      type(pic), dimension(nx,ny,nz) :: cell
      type(particle), dimension(:), pointer :: ghost
      integer, dimension(:), pointer :: sendlist
      integer, dimension(0:nproc_comm) :: ngp_send, ngp_recv
      integer :: npsend
!
      if (lparticles_blocks) call fatal_error('integrate_drag', 'particles_blocks version is not implemented yet. ')
!
!  Distribute particles to the cells.
!
      nullify(ghost)
      call distribute_particles(f, fp, npsend, sendlist, ghost, cell, ngp_send, ngp_recv)
!
!  Locally integrate the drag force as well as other source terms.
!
      if (ldrag_on_par .and. ldrag_on_gas) then
        call drag_on_both(cell, dt)
      elseif (ldrag_on_par) then
        call drag_on_particles(cell, dt)
      elseif (ldrag_on_gas) then
        call fatal_error('integrate_drag', 'drag on gas only is not implemented. ')
      endif
!
!  Assign the properties of the particles onto the grid.
!
      call map_particles(f, cell)
!
!  Collect the change in velocities.
!
      call collect_particles(f, fp, npsend, sendlist, ghost, cell, ngp_send, ngp_recv, &
                             ldrag_on_par, ldrag_on_gas .or. Omega /= 0.0, ldrag_on_gas .and. ldrag_pm_back_reaction)
!
    endsubroutine integrate_drag
!***********************************************************************
!***********************************************************************
! LOCAL SUBROUTINES START BELOW HERE.
!***********************************************************************
!***********************************************************************
    elemental subroutine drag_on_particles(cell, dt)
!
!  Finds the change in velocity for the particles under the drag.
!
!  11-dec-15/ccyang: coded.
!
      type(pic), intent(inout) :: cell
      real, intent(in) :: dt
!
!  Find the velocity change.
!
      rotation: if (Omega /= 0.0) then
        if (cell%np > 0) call drag_particle_omega(cell%p%v(1), cell%p%v(2), cell%u(1), cell%u(2), dt, cell%p%dv(1), cell%p%dv(2))
        call gas_epicycle(cell%u(1), cell%u(2), dt, cell%du(1), cell%du(2))
      else rotation
        call drag_particle(cell%p%v(1), cell%u(1), dt, cell%p%dv(1))
        call drag_particle(cell%p%v(2), cell%u(2), dt, cell%p%dv(2))
      endif rotation
!
      if (gz_par_coeff /= 0.0) then
        call drag_particle(cell%p%v(3), cell%u(3), dt, cell%p%dv(3), particle_zaccel(cell%p%x(3)))
      else
        call drag_particle(cell%p%v(3), cell%u(3), dt, cell%p%dv(3))
      endif
!
    endsubroutine drag_on_particles
!***********************************************************************
    elemental subroutine drag_on_both(cell, dt)
!
!  Finds the change in velocity for the gas and the particles under
!  mutual drags as well as other source terms.
!
!  13-sep-15/ccyang: coded.
!
      type(pic), intent(inout) :: cell
      real, intent(in) :: dt
!
      rotation: if (Omega /= 0.0) then
        call drag_mutual_omega(cell%p%eps, cell%p%v(1), cell%p%v(2), cell%u(1), cell%u(2), dt, &
                               cell%p%dv(1), cell%p%dv(2), cell%du(1), cell%du(2))
      else rotation
        if (gx_gas /= 0.0) then
          call drag_mutual(cell%p%eps, cell%p%v(1), cell%u(1), dt, cell%p%dv(1), cell%du(1), agas=gx_gas)
        else
          call drag_mutual(cell%p%eps, cell%p%v(1), cell%u(1), dt, cell%p%dv(1), cell%du(1))
        endif
        call drag_mutual(cell%p%eps, cell%p%v(2), cell%u(2), dt, cell%p%dv(2), cell%du(2))
      endif rotation
!
      if (gz_par_coeff /= 0.0) then
        call drag_mutual(cell%p%eps, cell%p%v(3), cell%u(3), dt, cell%p%dv(3), cell%du(3), apar=particle_zaccel(cell%p%x(3)))
      else
        call drag_mutual(cell%p%eps, cell%p%v(3), cell%u(3), dt, cell%p%dv(3), cell%du(3))
      endif
!
    endsubroutine drag_on_both
!***********************************************************************
    elemental subroutine drag_particle(v, u, dt, dv, ap)
!
!  Calculates the velocity change of each particle under the gas drag.
!
!  18-jun-15/ccyang: coded.
!
!  Input Arguments:
!      v        Initial particle velocity.
!      u        Gas velocity.
!      dt       Time step.
!  Output Arguments:
!      dv       Change in particle velocity in time dt.
!  Optional Argument:
!      ap       Background acceleration on the particle.
!
      use Sub, only: one_minus_exp
!
      real, intent(in) :: v, u, dt
      real, intent(out) :: dv
      real, intent(in), optional :: ap
!
      real :: t
!
      t = dt / tdrag
      getdv: if (present(ap)) then
        if (t * t > epsilon(1.0)) then
          dv = (u - v + ap * tdrag) * (1.0 - exp(-t))
        else
          dv = ((u - v) / tdrag + ap) * (1.0 - 0.5 * t) * dt
        endif
      else getdv
        dv = (u - v) * one_minus_exp(t)
      endif getdv
!
    endsubroutine drag_particle
!***********************************************************************
    pure subroutine drag_particle_omega(vx, vy, ux, uy, dt, dvx, dvy)
!
!  Calculates the velocity change of each particle under the gas drag
!  and rotation.
!
!  05-may-18/ccyang: coded.
!
!  Input Arguments:
!      vx       x-component of the initial velocities of the particles.
!      vy       y-component of the initial velocities of the particles.
!      ux       x-component of the gas velocity.
!      uy       y-component of the gas velocity.
!      dt       Time step.
!  Output Arguments:
!      dvx      Changes in x-component of particle velocities in time dt.
!      dvy      Changes in y-component of particle velocities in time dt.
!
      real, dimension(:), intent(in) :: vx, vy
      real, intent(in) :: ux, uy, dt
      real, dimension(:), intent(out) :: dvx, dvy
!
      real :: cosot, sinot1, sinot2
      real :: x, y, vxe, vye, vx0, vy0
!
!  Find the epicycle.
!
      x = epicycle_freq * dt
      cosot = cos(x)
      y = sin(x)
      sinot1 = epicycle_ratio * y
      sinot2 = y / epicycle_ratio
!
!  Find the decay factor.
!
      y = exp(-dt / tdrag)
      x = 1.0 - y
!
!  Find the terminal velocity.
!
      vye = -oneplustaus2inv * dv_gas
      vxe = twotaus * vye
!
      vx0 = vxe + (x * ux - y * vxe) * cosot + (x * (uy + dv_gas) - y * vye) * sinot1
      vy0 = vye + (x * (uy + dv_gas) - y * vye) * cosot - (x * ux - y * vxe) * sinot2
!
!  Find the velocity changes of the particles.
!
      x = 1.0 - y * cosot
      sinot1 = y * sinot1
      sinot2 = y * sinot2
!
      dvx = vx0 - (x * vx - vy * sinot1)
      dvy = vy0 - (x * vy + vx * sinot2)
!
    endsubroutine drag_particle_omega
!***********************************************************************
    pure subroutine drag_mutual(eps, v, u, dt, dv, du, apar, agas)
!
!  Calculates the velocity change of the gas and the particles under
!  mutual drag interaction.
!
!  13-sep-15/ccyang: coded.
!
!  Input Arguments:
!      eps      Particle-to-gas density ratio of each particle.
!      v        Initial particle velocity.
!      u        Initial gas velocity.
!      dt       Time step.
!  Output Arguments:
!      dv       Change in particle velocity in time dt.
!      du       Change in gas velocity in time dt, if
!          ldrag_pm_back_reaction = .false.; (1 + sum(eps)) times the
!          change in center-of-mass velocity, otherwise.
!  Optional Arguments:
!      apar     External acceleration of each of the particles.
!      agas     External acceleartion of the gas.
!
      use Sub, only: one_minus_exp
!
      real, dimension(:), intent(in) :: eps
      real, dimension(size(eps)), intent(in) :: v
      real, intent(in) :: u, dt
      real, dimension(size(eps)), intent(out) :: dv
      real, intent(out) :: du
      real, dimension(size(eps)), intent(in), optional :: apar
      real, intent(in), optional :: agas
!
      real :: epstot, norm, ucm, du0, acm, dag
      real :: t, x, y, z
!
!  Find weights.
!
      epstot = sum(eps)
      norm = 1.0 / (1.0 + epstot)
!
!  Find exponential decays.
!
      t = dt / tdrag
      x = one_minus_exp(t)
      if (.not. ldrag_pm_back_reaction) y = one_minus_exp((1.0 + epstot) * t)
      z = exp(-t) * one_minus_exp(epstot * t) / epstot
!
!  Find center-of-mass velocity.
!
      ucm = norm * (u + sum(eps * v))
      du0 = ucm - u
!
!  Find the velocity changes due to mutual drag.
!
      dv = (ucm - v) * x - du0 * z
      if (ldrag_pm_back_reaction) then
        if (.not. present(apar) .and. .not. present(agas)) du = 0.0
      else
        du = du0 * y
      endif
!
!  Add the velocity changes due to external accelerations.
!
      accel: if (present(apar) .or. present(agas)) then
!       Find center-of-mass acceleration.
        if (present(apar)) then
          acm = sum(eps * apar)
        else
          acm = 0.0
        endif
        if (present(agas)) acm = acm + agas
        if (ldrag_pm_back_reaction) du = acm * dt
        acm = norm * tdrag * acm
!       Find the velocity changes.
        if (present(apar)) dv = dv + (apar * tdrag - acm) * x
        if (present(agas)) then
          dag = norm * (agas * tdrag - acm)
        else
          dag = -norm * acm
        endif
        du0 = acm * t
        dv = dv + (du0 + dag * (x - z))
        if (.not. ldrag_pm_back_reaction) du = du + (du0 + dag * y)
      endif accel
!
    endsubroutine drag_mutual
!***********************************************************************
    pure subroutine drag_mutual_omega(eps, vx, vy, ux, uy, dt, dvx, dvy, dux, duy, agx)
!
!  Calculates the change in the horizontal components of velocity for
!  the gas and the particles under mutual drag interaction, shear,
!  Coriolis force, and background gas pressure gradient.
!
!  12-jun-15/ccyang: coded
!
!  Input Arguments:
!      eps      Particle-to-gas density ratio of each particle.
!      vx       x-component of the initial velocity of each particle.
!      vy       y-component of the initial velocity of each particle.
!      ux       x-component of the initial gas velocity.
!      uy       y-component of the initial gas velocity.
!      dt       Time step.
!  Output Arguments:
!      dvx      Change in x-component of particle velocity in time dt.
!      dvy      Change in y-component of particle velocity in time dt.
!      dux      Change in x-component of gas velocity in time dt, if
!          ldrag_pm_back_reaction = .false.; (1 + sum(eps)) times the
!          change in center-of-mass velocity, otherwise.
!      duy      Change in y-component of gas velocity in time dt, if
!          ldrag_pm_back_reaction = .false.; (1 + sum(eps)) times the
!          change in center-of-mass velocity, otherwise.
!  Optional Argument:
!      agx      x-component of the acceleartion of the gas.
!
      use Sub, only: one_minus_exp
!
      real, dimension(:), intent(in) :: eps
      real, dimension(size(eps)), intent(in) :: vx, vy
      real, intent(in) :: ux, uy, dt
      real, dimension(size(eps)), intent(out) :: dvx, dvy
      real, intent(out) :: dux, duy
      real, intent(in), optional :: agx
!
      real, dimension(size(eps)) :: vx0, vy0
      real :: epstot, dvg
      real :: ux0, uy0, uxe, uye, vxe, vye, vxcm, vycm
      real :: cosot, sinot1, sinot2
      real :: ot, t, ts, a0, a1, a2, a3, a4
!
!  Get the total solid-to-gas ratio.
!
      epstot = sum(eps)
!
!  Get the velocity reduction of the gas.
!
      if (present(agx)) then
        dvg = twoomega1 * (gx_gas + agx)
      else
        dvg = dv_gas
      endif
!
!  Get the deviation from the NSH solution.
!
      call get_nsh_solution(dvg, epstot, uxe, uye, vxe, vye)
      ux0 = ux - uxe
      uy0 = uy - uye
      vx0 = vx - vxe
      vy0 = vy - vye
!
!  Find the center-of-mass velocity of the particles.
!
      vcm: if (epstot > 0.0) then
        vxcm = sum(eps * vx0) / epstot
        vycm = sum(eps * vy0) / epstot
      else vcm
        vxcm = 0.0
        vycm = 0.0
      endif vcm
!
!  Get the epicycle motions.
!
      ot = epicycle_freq * dt
      cosot = cos(ot)
      a0 = sin(ot)
      sinot1 = a0 * epicycle_ratio
      sinot2 = a0 / epicycle_ratio
!
      uxe = ux0 * cosot + uy0 * sinot1
      uye = uy0 * cosot - ux0 * sinot2
!
      vxe = vxcm * cosot + vycm * sinot1
      vye = vycm * cosot - vxcm * sinot2
!
!  Find the velocity change for the particles.
!
      t = dt / tdrag
      a0 = exp(-t)
      a3 = 1.0 + epstot
      ts = a3 * t
      a4 = exp(-ts)
      if (abs(ts**4) > epsilon(1.0)) then
        a1 = (epstot + a4) / a3 - a0
      else
        a1 = 0.5 * epstot * t**2 * (1.0 - (t + ts) / 3.0)
      endif
      a2 = one_minus_exp(ts) / a3
!
      dvx = a1 * vxe + a2 * uxe + (a0 * (vx0 * cosot + vy0 * sinot1) - vx0)
      dvy = a1 * vye + a2 * uye + (a0 * (vy0 * cosot - vx0 * sinot2) - vy0)
!
!  Find the velocity change for the gas.
!
      gas: if (ldrag_pm_back_reaction) then
        uxe = ux0 + epstot * vxcm
        uye = uy0 + epstot * vycm
        a0 = ot * ot
        if (abs(a0) > epsilon(1.0)) then
          a0 = 1.0 - cosot
        else
          a0 = 0.5 * a0 * (1.0 - a0 / 12.0)
        endif
        dux = -a0 * uxe + sinot1 * uye
        duy = -a0 * uye - sinot2 * uxe
      else gas
        if (abs(ts * ts) > epsilon(1.0)) then
          a1 = (1.0 + epstot * a4) / a3
        else
          a1 = 1.0 - epstot * t
        endif
        a2 = epstot * a2
        dux = a1 * uxe + a2 * vxe - ux0
        duy = a1 * uye + a2 * vye - uy0
      endif gas
!
    endsubroutine drag_mutual_omega
!***********************************************************************
    pure subroutine gas_epicycle(ux, uy, dt, dux, duy)
!
!  Find the change in gas velocity due to rotation.
!
!  14-dec-15/ccyang: coded
!
!  Input Arguments
!      ux      x-component of the initial velocity of the gas.
!      uy      y-component of the initial velocity of the gas.
!      dt      Time step.
!  Output Arguments
!      dux     Change in x-component of the gas velocity in time dt.
!      duy     Change in y-component of the gas velocity in time dt.
!
      real, intent(in) :: ux, uy, dt
      real, intent(out) :: dux, duy
!
      real :: x, y, z
!
!  Find the epicycle.
!
      z = epicycle_freq * dt
      epicycle: if (z**4 > epsilon(1.0)) then
        x = 1.0 - cos(z)
        y = sin(z)
      else epicycle
        y = z**2
        x = 0.5 * y * (1.0 - y / 12.0)
        y = z * (1.0 - y / 6.0)
      endif epicycle
!
!  Find the velocity change of the gas.
!
      z = uy + dv_gas
      dux = -x * ux + y * epicycle_ratio * z
      duy = -x * z - y / epicycle_ratio * ux
!
    endsubroutine gas_epicycle
!***********************************************************************
    elemental subroutine get_nsh_solution(dv_gas, eps, ux, uy, vx, vy)
!
!  Finds the Nakagawa-Sekiya-Hayashi (1986) solution.
!
!  17-may-15/ccyang: coded.
!
!  Input Arguments
!      dv_gas  Velocity reduction of the gas due to radial pressure
!              gradient.
!      eps     Solid-to-gas mass ratio.
!
!  Output Arguments
!      ux      x-component of the velocity of the gas.
!      uy      y-component of the velocity of the gas.
!      vx      x-component of the velocity of the solids.
!      vy      y-component of the velocity of the solids.
!
      real, intent(in) :: dv_gas, eps
      real, intent(out) :: ux, uy, vx, vy
!
      real :: a, b
!
      a = 1.0 + eps
      b = dv_gas / (a**2 + taus2)
      vx = -twotaus * b
      vy = -a * b
      ux = -eps * vx
      uy = -(a + taus2) * b
!
    endsubroutine get_nsh_solution
!***********************************************************************
    elemental subroutine set_drag_equilibrium(cell)
!
!  Set the state of drag equilibrium to du and dv.
!
!  11-dec-15/ccyang: coded.
!
      type(pic), intent(inout) :: cell
!
      real, dimension(3) :: u, v
      integer :: i
      real :: eps
!
!  Find the equilibrium according to the drag type.
!
      u = 0.0
      v = 0.0
      eq: if (ldrag_on_par) then
        if (ldrag_on_gas) then
          eps = sum(cell%p%eps)
        else
          eps = 0.0
        endif
        call get_nsh_solution(dv_gas, eps, u(1), u(2), v(1), v(2))
        v = v + cell%u
      endif eq
!
!  Assign the solution to the velocity change.
!
      cell%du = u
      forall(i = 1:3) cell%p%dv(i) = v(i)
!
    endsubroutine set_drag_equilibrium
!***********************************************************************
!***********************************************************************
! LOCAL FUNCTIONS START BELOW HERE.
!***********************************************************************
!***********************************************************************
    function find_mp_swarm(eps_dtog) result(mp_swarm)
!
!  Finds the mass mp_swarm of each super-particle given the global dust-
!  to-gas density ratio eps_dtog.
!
!  13-sep-15/ccyang: coded.
!
      use EquationOfState, only: gamma, rho0, cs0
!
      real :: mp_swarm
      real, intent(in) :: eps_dtog
!
      real :: mass
!
!  Find the total gas mass.
!
      gasmass: if (gz_par_coeff > 0.0) then
!       For stratified gas
        mass = sqrt(twopi / gamma) * (cs0 / Omega) * product(lxyz(1:2), lactive_dimension(1:2))
      else gasmass
!       For uniform gas
        mass = product(lxyz, lactive_dimension)
      endif gasmass
      mass = rho0 * mass
!
!  Uniformly distribute the mass to the particles.
!
      mp_swarm = eps_dtog * mass / real(npar)
!
    endfunction find_mp_swarm
!***********************************************************************
    elemental real function particle_zaccel(zp) result(gz)
!
!  Returns the z-acceleration of a particle.
!
!  22-jun-15/ccyang: coded.
!
!  Currently, only supports -gz_par_coeff^2 * zp.
!
!  Input Argument:
!      zp  Vertical position of the particle.
!
      real, intent(in) :: zp
!
      gz = -gzpc * zp
!
    endfunction particle_zaccel
!***********************************************************************
endmodule Particles_drag