! $Id: gravity_r.f90,v 1.1 2018/08/24 15:48:10 wlyra Exp $
!
!  Radial gravity
!
!** 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 :: lgrav = .true.
!
! MVAR CONTRIBUTION 0
! MAUX CONTRIBUTION 0
! MGLOBAL CONTRIBUTION 3
!
! PENCILS PROVIDED gg(3)
!
!***************************************************************
module Gravity
!
  use Cparam
  use Cdata
  use General, only: keep_compiler_quiet
  use Messages
!
  implicit none
!
  include 'gravity.h'
!
  interface potential
    module procedure potential_global
    module procedure potential_penc
    module procedure potential_point
  endinterface
!
  interface acceleration
    module procedure acceleration_penc
    module procedure acceleration_penc_1D
    module procedure acceleration_point
  endinterface
!
  ! coefficients for potential
  real, dimension (5,ninit) :: cpot=0. !=(/ 0., 0., 0., 0., 0. /)
  real, dimension (16,ninit) :: cpot2=0.
  real, dimension(ninit) :: g01=0.,rpot=0.
  real :: lnrho_bot,lnrho_top,ss_bot,ss_top
  real :: gravz_const=1.,reduced_top=1.
  real :: g0=0.
  real, target :: g1=0.,rp1=1.0,rp1_smooth=impossible
  real, target :: gsum=0.
  real :: rp1_smooth1,frac_smooth=1.0
  real :: r0_pot=0.,r1_pot1=0.    ! peak radius for smoothed potential
  real :: n_pot=10,n_pot1=10   ! exponent for smoothed potential
  real :: qgshear=1.5  ! (global) shear parameter
                       !     1.5 for Keplerian disks, 1.0 for galaxies
!
  character (len=labellen), dimension(ninit) :: ipotential='zero'
!
  ! variables for compatibility with grav_z (used by Entropy and Density):
  real :: z1,z2,zref,zgrav,gravz=0.,zinfty
  real :: nu_epicycle=1.0
  real, target :: t_ramp_mass=impossible
  real :: t1_ramp_mass
  character (len=labellen) :: gravz_profile='zero'
  character (len=labellen) :: ipotential_secondary='plummer'
  character (len=labellen) :: iramp_function='linear'
!
  logical :: lnumerical_equilibrium=.false.
  logical :: lgravity_gas=.true.
  logical :: lgravity_neutrals=.true.
  logical :: lgravity_dust=.true.
  logical :: lindirect_terms=.true.
  logical, target :: lramp_mass=.false.
  integer :: iglobal_gg=0
  logical, target :: lsecondary_wait=.false.
  logical :: lcoriolis_force_gravity=.true.
  logical :: lcentrifugal_force_gravity=.true.
  real, target :: t_start_secondary = -impossible
!
  namelist /grav_init_pars/ &
      ipotential,g0,r0_pot,r1_pot1,n_pot,n_pot1,lnumerical_equilibrium, &
      qgshear,lgravity_gas,g01,rpot,gravz_profile,gravz,nu_epicycle, &
      lgravity_neutrals,g1,rp1_smooth,lindirect_terms,lramp_mass,t_ramp_mass, &
      ipotential_secondary,lsecondary_wait,t_start_secondary,iramp_function, &
      lcoriolis_force_gravity,lcentrifugal_force_gravity,frac_smooth
!
  namelist /grav_run_pars/ &
      ipotential,g0,r0_pot,n_pot,lnumerical_equilibrium, &
      qgshear,lgravity_gas,g01,rpot,gravz_profile,gravz,nu_epicycle, &
      lgravity_neutrals,g1,rp1_smooth,lindirect_terms,lramp_mass,t_ramp_mass, &
      ipotential_secondary,lsecondary_wait,t_start_secondary,iramp_function, &
      lcoriolis_force_gravity,lcentrifugal_force_gravity,frac_smooth
!
  integer :: idiag_torque=0
!
  contains
!***********************************************************************
    subroutine register_gravity
!
!  initialise gravity flags
!
!  10-jan-02/wolf: coded
!
!  Identify version number.
!
      use FArrayManager
      if (lroot) call svn_id("$Id: gravity_r.f90,v 1.1 2018/08/24 15:48:10 wlyra Exp $")
!
      lgravr=.true.
      lgravr_gas =.true.
      lgravr_dust=.true.
      lgravr_neutrals=.true.
!
!  Register global_gg, so we can later retrieve gravity via get_global.
!  Ensure the reserved f-array slots are initialized to zero, so we can add
!  ninit different gravity fields (done in initialize_gravity).
!
      if (.not.lnumerical_equilibrium) then
        if (iglobal_gg==0) call farray_register_global('global_gg',iglobal_gg,vector=3)
      endif
!
    endsubroutine register_gravity
!***********************************************************************
    subroutine initialize_gravity(f)
!
!  Set up cpot according to the value of ipotential, and initialize the
!  global variable gg (gravity field).
!  Needed by both start.f90 and run.f90
!
!  16-jul-02/wolf: coded
!  22-nov-02/tony: renamed
!  15-mar-07/wlad: made it coordinate-independent
!
      use Sub, only: poly, step, get_radial_distance
      use Mpicomm
      use SharedVariables, only: put_shared_variable
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3) :: gg_mn=0.0
      real, dimension (nx)   :: g_r,rr_mn,rr_sph,rr_cyl,pot
      logical :: lpade=.true. ! set to false for 1/r potential
      integer :: j
!
!  Pre-calculate the angular frequency of the rotating frame in the case a
!  secondary stationary body is added to the simulation.
!
      if (lcorotational_frame) then 
        gsum=g0+g1
        Omega_corot = sqrt(gsum/rcorot**3)
        rp1=rcorot
      else
        gsum=g0
        if (g1/=0) call fatal_error("initialize_gravity",&
             "companion gravity coded only for corotational frame")
      endif
!
!  Smoothing radius: default to Hill radius
!
      if (rp1_smooth == impossible) &
           rp1_smooth   = frac_smooth * rp1 * (g1/3.)**(1./3)
      if (rp1_smooth /= 0) then
        rp1_smooth1  = 1./rp1_smooth
      else
        rp1_smooth1 = 0.
      endif
!
!  Share variables related to corotational frame to modules that may need it
!
      call put_shared_variable('gsum',gsum,caller='initialize_gravity')
      call put_shared_variable('g0',g0)
      call put_shared_variable('g1',g1)
      call put_shared_variable('rp1',rp1)
      call put_shared_variable('rp1_smooth',rp1_smooth)
      call put_shared_variable('lramp_mass',lramp_mass)
      call put_shared_variable('t_ramp_mass',t_ramp_mass)
      call put_shared_variable('lsecondary_wait',lsecondary_wait)
      call put_shared_variable('t_start_secondary',t_start_secondary)
!
!  Shortcut for optimization
!
      if (t_ramp_mass/=impossible) t1_ramp_mass=1./t_ramp_mass
!
!  for lpade=.true. set coefficients for potential (coefficients a0, a2, a3,
!  b2, b3) for the rational approximation
!
!              a_0   +   a_2 r^2 + a_3 r^3
!    Phi(r) = ---------------------------------------
!               1    +   b_2 r^2 + b_3 r^3 + a_3 r^4
!
      if (lnumerical_equilibrium) then
!
        if (lroot) then
          print*,'initialize_gravity: numerical exact equilibrium'
          print*,' - gravity will be calculated in density module'
        endif
!
      else

        f(l1:l2,m1:m2,n1:n2,iglobal_gg:iglobal_gg+2) = 0.
!
        do j=1,1
!
          lpade=.true.
!
          select case (ipotential(j))
!
          case ('zero')           ! zero potential
            if (lroot) print*, 'initialize_gravity: zero gravity potential'
            cpot(:,j) = 0.
!
          case ('solar')          ! solar case
            if (lroot) print*, 'initialize_gravity: solar gravity potential'
            cpot(:,j) = (/ 5.088, -4.344, 61.36, 10.91, -13.93 /)
!
          case ('M5-dwarf')       ! M5 dwarf
            if (lroot) print*, 'initialize_gravity: M5 dwarf gravity potential'
            cpot(:,j) = (/ 2.3401, 0.44219, 2.5952, 1.5986, 0.20851 /)
!
!  Experimental flattened potential for M5 star. Similar to M5-dwarf but
!  includes 6th and 12th powers of r to obtain pot~const for r \gtrsim 1.2.
!  Here a 15th order polynomial fit is used to construct the potential.
!
          case ('M5-flat')        ! M5 flat
            if (lroot) print*, 'initialize_gravity: M5 flat gravity potential'
            cpot2(:,j) = (/ 2.33932201e+00, 1.65736990e-01, -9.66534918e+00, 9.93311842e+01, &
                           -7.84474588e+02, 3.80654218e+03, -1.20467561e+04, 2.60329111e+04, &
                           -3.94606849e+04, 4.25603487e+04, -3.27771332e+04, 1.78753221e+04, &
                           -6.73781490e+03, 1.66874150e+03, -2.44250835e+02, 1.60061840e+01 /)
            lpade=.false.
!
          case ('M2-sgiant')       ! M super giant
            if (lroot) print*, 'initialize_gravity: M super giant gravity potential'
            cpot(:,j) = (/ 1.100, 0.660, 2.800, 1.400, 0.100 /)
!
          case ('A7-star')       ! Ap star
            if (lroot) print*, 'initialize_gravity: A star gravity potential'
            cpot(:,j) = (/ 4.080, -3.444, 15.2000, 11.2000, -12.1000 /)
!
          case ('A0-star')       ! A0 star
            if (lroot) print*, 'initialize_gravity: A0 star gravity potential'
            !cpot(:,j) = (/ 4.7446,  -1.9456,  0.6884,  4.8007, 1.79877 /)
            cpot(:,j) = (/ 4.3641,  -1.5612,  0.4841, 4.0678, 1.2548 /)
!
          case ('simple')         ! simple potential for tests
            if (lroot) print*, 'initialize_gravity: very simple gravity potential'
            cpot(:,j) =  (/ 1., 0., 0., 1., 0. /)
!
          case ('simple-2')       ! another simple potential for tests
            if (lroot) print*, 'initialize_gravity: simple gravity potential'
            cpot(:,j) =  (/ 1., 1., 0., 1., 1. /)
!
          case ('smoothed-newton')
            if (lroot) print*,'initialize_gravity: smoothed 1/r potential'
            lpade=.false.
!
          case ('no-smooth','newton','newtonian')
            if (lroot) print*,'initialize_gravity: non-smoothed newtonian gravity'
            lpade=.false.
!
          case ('varying-q')
            if (lroot) print*,'initialize_gravity: shear with Omega proto r^-q, q=',qgshear
            lpade=.false.
!
          case ('varying-q-smooth')
            if (lroot) &
                 print*,'initialize_gravity: shear with smoothed Omega proto r^-q, q=',qgshear
            lpade=.false.
!            
          case ('dark-matter-halo')
            if (lroot) &
                 print*,'initialize_gravity: arc-tangent potential generated by a dark matter halo'
            lpade=.false.
!
          case ('light-matter')
            if (lroot) &
                 print*,'initialize_gravity: potential generated by an exponential baryonic disk'
            lpade=.false.
!
            ! geodynamo
          case ('geo-kws-approx')     ! approx. 1/r potential between r=.5 and r=1
            if (lroot) print*, 'initialize_gravity: approximate 1/r potential'
            cpot(:,j) = (/ 0., 2.2679, 0., 0., 1.1697 /)
!
          case ('geo-benchmark')      ! for geodynamo benchmark runs
            if (lroot) print*, 'initialize_gravity: gravity linear in radius'
            cpot(:,j) = (/ 0., -.5, 0., 0., 0. /)
!
          case ('geo-kws')
            if (lroot) print*, 'initialize_gravity: '//&
                 'smoothed 1/r potential in spherical shell'
            if (r0_pot < epsi) print*, 'WARNING: grav_r: r0_pot is too small.'//&
                 'Can be set in grav_r namelists.'
            lpade=.false.
            ! end geodynamo
!
          case default
!
!  Catch unknown values
!
            if (lroot) print*, 'initialize_gravity: '//&
                 'No such value for ipotential: ', trim(ipotential(j))
            call stop_it("")
!
          endselect
!
          do n=n1,n2
            do m=m1,m2
!
!  rr_mn differs depending on the coordinate system
!
              call get_radial_distance(rr_sph,rr_cyl)
!
!  choose between sphere-in-box and cylindrical gravity
!
              if (lcylindrical_gravity) then
                rr_mn=rr_cyl
              else
! sphere in a box
                rr_mn=rr_sph
              endif
!
              if (lpade) then
!
                g_r = - rr_mn * poly( (/ 2*(cpot(1,j)*cpot(4,j)-cpot(2,j)), &
                     3*(cpot(1,j)*cpot(5,j)-cpot(3,j)), &
                     4*cpot(1,j)*cpot(3,j), &
                     cpot(5,j)*cpot(2,j)-cpot(3,j)*cpot(4,j), &
                     2*cpot(2,j)*cpot(3,j), &
                     cpot(3,j)**2  /), rr_mn) &
                     / poly( (/ 1., 0., cpot(4,j), cpot(5,j), &
                     cpot(3,j) /), rr_mn)**2
              else
                if (ipotential(j) == 'no-smooth'.or.&
                    ipotential(j) == 'newton'.or.&
                    ipotential(j) == 'newtonian') then
                  g_r=-g0/rr_mn**2
                elseif (ipotential(j) == 'varying-q') then
                  g_r=-g0/rr_mn**(2*qgshear-1)
                elseif (ipotential(j) == 'varying-q-smooth') then
                  g_r=-g0*rr_mn/(rr_mn**2+r0_pot**2)**qgshear
                elseif (ipotential(j) == 'dark-matter-halo') then
                  g_r=-g01(j)*(1-rpot(j)/rr_mn*atan2(rr_mn,rpot(j)))/rr_mn
                elseif (ipotential(j) == 'light-matter') then
                  !approximation of the bessel functions potential
                  !of Freeman 1970. Reference: Persic et al 1996
                  g_r=-4.134e-4*g01(j)*(.5*rr_mn/rpot(j))**1.22/&
                       ((.5*rr_mn/rpot(j))**2+1.502)**1.43/rr_mn
                elseif (ipotential(j) == 'M5-flat') then
                  !experimental flattened potential for M5 star
                   g_r = poly( (/   cpot2(2,j),  2*cpot2(3,j),  3*cpot2(4,j),  4*cpot2(5,j),  &
                                  5*cpot2(6,j),  6*cpot2(7,j),  7*cpot2(8,j),  8*cpot2(9,j),  &
                                  9*cpot2(10,j),10*cpot2(11,j),11*cpot2(12,j),12*cpot2(13,j), &
                                 13*cpot2(14,j),14*cpot2(15,j),15*cpot2(16,j) /), rr_mn)
                else
!
!  smoothed 1/r potential in a spherical shell
!  r0_pot is the smoothing radius, and n_pot the smoothing exponent
!
                  !g_r=-g0*rr_mn**(n_pot-1) &
                  !     *(rr_mn**n_pot+r0_pot**n_pot)**(-1./n_pot-1.)
                  pot = -g0*(1. + (r1_pot1*rr_mn)**n_pot1)**(1.0/n_pot1) &
                           /(rr_mn**n_pot + r0_pot**n_pot)**(1.0/n_pot)
                  g_r = pot*(rr_mn**(n_pot-1)/(rr_mn**n_pot+r0_pot**n_pot) &
                           - r1_pot1*(r1_pot1*rr_mn)**(n_pot1-1) &
                                    /(1 + (r1_pot1*rr_mn)**n_pot1))
                endif
              endif
!
              call get_gravity_field(g_r,gg_mn,rr_mn)
              f(l1:l2,m,n,iglobal_gg:iglobal_gg+2) = &
                  f(l1:l2,m,n,iglobal_gg:iglobal_gg+2) + gg_mn
!
            enddo
          enddo
!
        enddo
      endif
!
!  Add a simple vertical gravity as well.
!
      select case (gravz_profile)
!
      case ('zero')
!
      case ('const')
        if (lroot) print*,'initialize_gravity: constant gravz=', gravz
        do n=n1,n2; do m=m1,m2
          f(l1:l2,m,n,iglobal_gg+2)=f(l1:l2,m,n,iglobal_gg+2)+gravz
        enddo; enddo
!
      case ('linear')
        if (lroot) print*,'initialize_gravity: linear z-grav, nu=', nu_epicycle
        do n=n1,n2; do m=m1,m2
          f(l1:l2,m,n,iglobal_gg+2)=f(l1:l2,m,n,iglobal_gg+2)-nu_epicycle**2*z(n)
        enddo; enddo
!
      case default
        if (lroot) print*, &
            'initialize_gravity: unknown gravz_profile ', gravz_profile
        call fatal_error('initialize_gravity','chosen gravz_profile not valid')
!
      endselect
!
    endsubroutine initialize_gravity
!***********************************************************************
    subroutine read_gravity_init_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=grav_init_pars, IOSTAT=iostat)
!
    endsubroutine read_gravity_init_pars
!***********************************************************************
    subroutine write_gravity_init_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=grav_init_pars)
!
    endsubroutine write_gravity_init_pars
!***********************************************************************
    subroutine read_gravity_run_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=grav_run_pars, IOSTAT=iostat)
!
    endsubroutine read_gravity_run_pars
!***********************************************************************
    subroutine write_gravity_run_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=grav_run_pars)
!
    endsubroutine write_gravity_run_pars
!***********************************************************************
    subroutine init_gg(f)
!
!  initialise gravity; called from start.f90
!  10-jan-02/wolf: coded
!  24-nov-02/tony: renamed from init_grav for consistancy (i.e. init_[variable name])
!
      real, dimension (mx,my,mz,mfarray) :: f
!
! Not doing anything (this might change if we decide to save gg to a file)
!
      call keep_compiler_quiet(f)
!
    endsubroutine init_gg
!***********************************************************************
    subroutine pencil_criteria_gravity
!
!  All pencils that the Gravity module depends on are specified here.
!
!  20-11-04/anders: coded
!
      if (lcorotational_frame) lpenc_requested(i_uu)=.true.
      if (idiag_torque/=0)     lpenc_diagnos(i_rho) =.true.
!
    endsubroutine pencil_criteria_gravity
!***********************************************************************
    subroutine pencil_interdep_gravity(lpencil_in)
!
!  Interdependency among pencils from the Gravity module is specified here.
!
!  20-11-04/anders: coded
!
      logical, dimension(npencils) :: lpencil_in
!
      if ((lhydro.and.lgravity_gas).or.&
          (lneutralvelocity.and.lgravity_neutrals).or.&
          (ldustvelocity.and.lgravity_dust)) &
          lpencil_in(i_gg)=.true.
!
      call keep_compiler_quiet(lpencil_in)
!
    endsubroutine pencil_interdep_gravity
!***********************************************************************
    subroutine calc_pencils_gravity(f,p)
!
!  Calculate Gravity pencils.
!  Most basic pencils should come first, as others may depend on them.
!
!  12-nov-04/anders: coded
!
      use FArrayManager, only: farray_use_global
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3):: ggp
      type (pencil_case) :: p
      integer, pointer :: iglobal_gg
!
      intent(in) :: f
!
      call farray_use_global('global_gg',iglobal_gg)
!
      if (lpencil(i_gg)) then
        call farray_use_global('global_gg',iglobal_gg)
        p%gg = f(l1:l2,m,n,iglobal_gg:iglobal_gg+2)
!
!  If there is a secondary body whose mass is changing in time (ramped-up), then 
!  this gravity is not added to the global array. It is re-calculated instead. 
!
        if (g1/=0 .and. (lramp_mass.or.lsecondary_wait)) then
          call secondary_body_gravity(ggp)
          p%gg = p%gg + ggp
        endif
      endif
!
    endsubroutine calc_pencils_gravity
!***********************************************************************
    subroutine duu_dt_grav(f,df,p)
!
!  add duu/dt according to gravity
!
!  10-jan-02/wolf: coded
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx,my,mz,mvar) :: df
      type (pencil_case) :: p
      integer :: k
!
! if statement for testing purposes
!
      if (lgravity_gas) then
        df(l1:l2,m,n,iux:iuz) = df(l1:l2,m,n,iux:iuz) + p%gg
      endif
!
      if (lneutralvelocity.and.lgravity_neutrals) then
        df(l1:l2,m,n,iunx:iunz) = df(l1:l2,m,n,iunx:iunz) + p%gg
      endif
!
      if (ldustvelocity.and.lgravity_dust) then
        do k=1,ndustspec
          df(l1:l2,m,n,iudx(k):iudz(k)) = df(l1:l2,m,n,iudx(k):iudz(k)) + p%gg
        enddo
      endif
!
!  Indirect term for binary systems with origin at the primary.
!
      if (lcorotational_frame) call indirect_plus_inertial_terms(df,p)
!
      if (ldiagnos) then
         if (idiag_torque/=0) call calc_torque(p)
      endif
!
      call keep_compiler_quiet(f)
!
    endsubroutine duu_dt_grav
!***********************************************************************
    subroutine gravity_after_boundary(f)
!
!  For actions outside mn-loop. 
!
!  9-jun-15/MR: coded
!
      real, dimension(mx,my,mz,mfarray) :: f

      call keep_compiler_quiet(f)

    endsubroutine gravity_after_boundary
!***********************************************************************
    subroutine indirect_plus_inertial_terms(df,p)
!
!  Add the indirect terms from the motion of the primary in a non-inertial frame,
!  plus the Coriolis and centrifugal terms from the motion of the secondary.
!
!  The indirect potential is Phi = GMp/rp**3 r dot rp
!
!  Mp is planet mass and rp planet position. For rp=(1,0,0), the potential is
!
!    Phi = gp * x = gp * r*cos(phi)
!
!  with gp = GMp/rp**2. So the resulting acceleration is, in cylindrical: 
!
!    grad = -       dPhi/drad = - gp * cos(phi)
!    gphi = - 1/r * dPhi/dphi =   gp * sin(phi)
!    gzed = -       dPhi/dzed =   0.
!
!  In spherical coordinates, x = r*sin(tht)*cos(phi), so 
!
!    grad = -                  dPhi/drad = - gp * sin(tht)*cos(phi)
!    gtht = - 1/ r           * dPhi/dtht = - gp * cos(tht)*cos(phi)
!    gphi = - 1/[r*sin(tht)] * dPhi/dphi =   gp *          sin(phi)  
!
!  20-jul-14/wlad: coded
!
      real, dimension (mx,my,mz,mvar) :: df
      real, dimension (nx) :: rrcyl_mn
      real :: c2,s2,g2
      type (pencil_case) :: p
!
      g2 = g1/rp1**2
!
      if (lramp_mass) call rampup_secondary_mass(g2)
!
!  Do not allow secondary before t_start_secondary
!
      if (lsecondary_wait .and. t <= t_start_secondary) then
        g2 = 0. 
      endif
!
      if (lcylindrical_coords) then
!
!  Indirect terms
!
        if (lindirect_terms) then
          df(l1:l2,m,n,iux) = df(l1:l2,m,n,iux) - g2*cos(y(m))
          df(l1:l2,m,n,iuy) = df(l1:l2,m,n,iuy) + g2*sin(y(m))
        endif
!
!  Coriolis force
!
        if (lcoriolis_force_gravity) then
          c2 = 2*Omega_corot
          df(l1:l2,m,n,iux) = df(l1:l2,m,n,iux) + c2*p%uu(:,2)
          df(l1:l2,m,n,iuy) = df(l1:l2,m,n,iuy) - c2*p%uu(:,1)
        endif
!
!  Centrifugal force
!
        if (lcentrifugal_force_gravity) &
             df(l1:l2,m,n,iux) = df(l1:l2,m,n,iux) + x(l1:l2)*Omega_corot**2
!
      else if (lspherical_coords) then
!
!  Indirect terms
!
        if (lindirect_terms) then
          df(l1:l2,m,n,iux) = df(l1:l2,m,n,iux) - g2*sinth(m)*cosph(n)
          df(l1:l2,m,n,iuy) = df(l1:l2,m,n,iuy) - g2*costh(m)*cosph(n)
          df(l1:l2,m,n,iuz) = df(l1:l2,m,n,iuz) + g2*         sinph(n)
        endif
!
!  Coriolis force
!
        if (lcoriolis_force_gravity) then
          c2 = 2*Omega_corot*costh(m)
          s2 = 2*Omega_corot*sinth(m)
!
          df(l1:l2,m,n,iux) = df(l1:l2,m,n,iux) + s2*p%uu(:,3)
          df(l1:l2,m,n,iuy) = df(l1:l2,m,n,iuy) + c2*p%uu(:,3)
          df(l1:l2,m,n,iuz) = df(l1:l2,m,n,iuz) - c2*p%uu(:,2) - s2*p%uu(:,1)
        endif
!
!  Centrifugal force
!
        if (lcentrifugal_force_gravity) then
          rrcyl_mn=x(l1:l2)*sinth(m)
          df(l1:l2,m,n,iux) = df(l1:l2,m,n,iux) + rrcyl_mn*sinth(m)*Omega_corot**2
          df(l1:l2,m,n,iuy) = df(l1:l2,m,n,iuy) + rrcyl_mn*costh(m)*Omega_corot**2
        endif
!
      endif
!
    endsubroutine indirect_plus_inertial_terms
!***********************************************************************
    subroutine rampup_secondary_mass(gp)
!
!  Increase the mass of the secondary gradually, linearly or 
!  sinusoidally, to prevent too strong an impulse in the
!  beginning of the simulation.
!
!  15-nov-21/wlyra: coded
!
      real :: gp

      if (t<=t_ramp_mass) then
!
        select case (iramp_function)
!
        case ('linear')
           gp = gp * t*t1_ramp_mass
        case ('sinusoidal')
           gp = gp * (sin(t*t1_ramp_mass*.5*pi))**2
        case default
           call fatal_error("rampup_secondary_mass",&
                "no ramping function selected")
        endselect
      endif
!
    endsubroutine rampup_secondary_mass
!***********************************************************************
    subroutine potential_global(pot,pot0)
!
!  gravity potential; version called by init_hydro, which operates on
!  full global coordinate arrays
!
!  16-jul-02/wolf: coded
!
      use Sub,     only: poly
      use Mpicomm, only: stop_it
!
      real, dimension (:,:,:) :: pot
      real, optional :: pot0           ! potential at r=0
      real, dimension (size(pot,1),size(pot,2),size(pot,3)) :: rr
      integer :: j
!
      intent(out) :: pot,pot0
!
      if (lcylindrical_gravity) &
           call stop_it("gravity_r: potential global not implemented "//&
           "for cylindrical gravity approximation")
!
!  remove this if you are sure rr is already calculated elsewhere
!
      if     (coord_system=='cartesian') then
        do n=n1,n2; do m=m1,m2
          rr(l1:l2,m,n)=sqrt(x(l1:l2)**2+y(m)**2+z(n)**2)
        enddo; enddo
      elseif (coord_system=='cylindric') then
        do n=n1,n2; do m=m1,m2
          rr(l1:l2,m,n)=sqrt(x(l1:l2)**2+z(n)**2)
        enddo; enddo
      elseif (coord_system=='spherical') then
        do n=n1,n2; do m=m1,m2
          rr(l1:l2,m,n)=x(l1:l2)
        enddo; enddo
      endif
!
      pot=0.
      if (present(pot0)) pot0=0.
      do j=1,ninit
        select case (ipotential(j))
!
        case ('geo-kws','smoothed-newton')
          !pot = pot -g0*(rr**n_pot+r0_pot**n_pot)**(-1.0/n_pot)
          pot = pot - g0*(1. + (r1_pot1*rr)**n_pot1)**(1.0/n_pot1) &
                        /(rr**n_pot + r0_pot**n_pot)**(1.0/n_pot)
          if (present(pot0)) pot0=pot0-g0/r0_pot
!
        case ('no-smooth','newton','newtonian')
          pot = pot -g0/rr
        case ('M5-flat')
           pot = pot - poly( (/ cpot2(1,j),  cpot2(2,j),  cpot2(3,j),  cpot2(4,j),  &
                                cpot2(5,j),  cpot2(6,j),  cpot2(7,j),  cpot2(8,j),  &
                                cpot2(9,j),  cpot2(10,j), cpot2(11,j), cpot2(12,j), &
                                cpot2(13,j), cpot2(14,j), cpot2(15,j), cpot2(16,j) /), rr)
          if (present(pot0)) pot0 = pot0 - cpot2(1,j)
        case default
          pot = pot - poly((/cpot(1,j), 0., cpot(2,j), cpot(3,j)/), rr) &
               / poly((/1., 0., cpot(4,j), cpot(5,j), cpot(3,j)/), rr)
          if (present(pot0)) pot0 = pot0-cpot(1,j)
!
        endselect
      enddo
!
    endsubroutine potential_global
!***********************************************************************
    subroutine potential_penc(xmn,ymn,zmn,pot,pot0,grav,rmn)
!
!  Gravity potential along one pencil
!
!  21-jan-02/wolf: coded
!  29-aug-07/wlad: allow arrays of arbitrary dimension
!
      use Sub,     only: poly
      use Mpicomm, only: stop_it
!
      real, dimension (:) :: pot
      real, dimension (size(pot)) :: rad
      real, optional :: ymn,zmn,pot0
      real, optional, dimension (size(pot)) :: xmn,rmn
      real, optional, dimension (size(pot),3) :: grav
      integer :: j
!
      intent(in)  :: xmn,ymn,zmn,rmn
      intent(out) :: pot,pot0,grav
!
      if (present(rmn)) then
        rad = rmn
      else
        if (present(xmn) .and. present(ymn) .and. present(zmn)) then
!
          if (.not.lcartesian_coords) &
              call stop_it("gravity_r: potential_penc with xmn,ymn,zmn is "//&
              "not yet implemented for non-cartesian coordinates. Fix the call "//&
              "to use radial distance instead")
!
          rad = sqrt(xmn**2+ymn**2+zmn**2)
        else
          call stop_it("POTENTIAL_PENC: Need to specify either x,y,z or r.")
        endif
      endif
!
      pot=0.
      if (present(pot0)) pot0=0.
      do j=1,ninit
        select case (ipotential(j))
!
        case ('geo-kws','smoothed-newton')
          !pot=pot-g0*(rmn**n_pot+r0_pot**n_pot)**(-1.0/n_pot)
          pot = pot - g0*(1. + (r1_pot1*rmn)**n_pot1)**(1.0/n_pot1) &
                        /(rmn**n_pot + r0_pot**n_pot)**(1.0/n_pot)
          if (present(pot0)) pot0=pot0-g0/r0_pot
!
        case ('no-smooth','newton','newtonian')
          pot=pot-g0/rmn
!
        case ('M5-flat')
           pot = pot - poly( (/ cpot2(1,j),  cpot2(2,j),  cpot2(3,j),  cpot2(4,j),  &
                                cpot2(5,j),  cpot2(6,j),  cpot2(7,j),  cpot2(8,j),  &
                                cpot2(9,j),  cpot2(10,j), cpot2(11,j), cpot2(12,j), &
                                cpot2(13,j), cpot2(14,j), cpot2(15,j), cpot2(16,j) /), rad)
          if (present(pot0)) pot0 = pot0 - cpot2(1,j)
        case default
          pot = pot - poly((/cpot(1,j), 0., cpot(2,j), cpot(3,j)/), rad) &
               / poly((/1., 0., cpot(4,j), cpot(5,j), cpot(3,j)/), rad)
          if (present(pot0)) pot0=pot0-cpot(1,j)
!
        endselect
      enddo
!
      if (present(grav)) then
        call not_implemented("potential_penc", "optional argument grav")
        call keep_compiler_quiet(grav)
      endif
!
    endsubroutine potential_penc
!***********************************************************************
    subroutine potential_point(x,y,z,r, pot,pot0, grav)
!
!  Gravity potential in one point
!
!  20-dec-03/wolf: coded
!
      use Sub,     only: poly
      use Mpicomm, only: stop_it
!
      real :: pot,rad
      real, optional :: x,y,z,r
      real, optional :: pot0,grav
      integer :: j
!
      intent(in)  :: x,y,z,r
      intent(out) :: pot,pot0,grav
!
      if (present(r)) then
        rad = r
      else
        if (present(x) .and. present(y) .and. present(z)) then
!
          if (.not.lcartesian_coords) &
               call stop_it("gravity_r: potential_point with x,y,z is "//&
               "not yet implemented for non-cartesiand coordinates. Fix the call "//&
               "to  use radial distance instead")
!
          rad = sqrt(x**2+y**2+z**2)
        else
          call stop_it("Need to specify either x,y,z or r in potential_point()")
        endif
      endif
!
      pot=0.
      if (present(pot0)) pot0=0.
      do j=1,ninit
        select case (ipotential(j))
!
        case ('geo-kws','smoothed-newton')
          !pot=pot-g0*(rad**n_pot+r0_pot**n_pot)**(-1.0/n_pot)
          pot = pot - g0*(1. + (r1_pot1*rad)**n_pot1)**(1.0/n_pot1) &
                        /(rad**n_pot + r0_pot**n_pot)**(1.0/n_pot)
          if (present(pot0)) pot0=pot0-g0/r0_pot
!
        case ('no-smooth','newton','newtonian')
          pot=pot-g0/r
!
        case default
          pot = pot- poly((/cpot(1,j), 0., cpot(2,j), cpot(3,j)/), rad) &
               / poly((/1., 0., cpot(4,j), cpot(5,j), cpot(3,j)/), rad)
          if (present(pot0)) pot0=pot0-cpot(1,j)
!
        case ('M5-flat')
           pot = pot - poly( (/ cpot2(1,j),  cpot2(2,j),  cpot2(3,j),  cpot2(4,j),  &
                                cpot2(5,j),  cpot2(6,j),  cpot2(7,j),  cpot2(8,j),  &
                                cpot2(9,j),  cpot2(10,j), cpot2(11,j), cpot2(12,j), &
                                cpot2(13,j), cpot2(14,j), cpot2(15,j), cpot2(16,j) /), rad)
          if (present(pot0)) pot0 = pot0 - cpot2(1,j)
        endselect
      enddo
!
      if (present(grav)) then
        call not_implemented("potential_point", "optional argument grav")
        call keep_compiler_quiet(grav)
      endif
!
    endsubroutine potential_point
!***********************************************************************
    subroutine acceleration_penc(gg)
!
!  Calculates gravitational acceleration on a pencil
!
!  21-apr-07/tobi: adapted from potential_penc
!
      use Messages, only: fatal_error
!
      real, dimension (:,:), intent (out) :: gg
!
!  Calculate acceleration from master pencils defined in initialize_gravity
!
      call fatal_error("acceleration_penc","Not implemented")
      call keep_compiler_quiet(gg)
!
    endsubroutine acceleration_penc
!***********************************************************************
    subroutine acceleration_penc_1D(g_r)
!
!  Gravitational acceleration along one pencil
!
!  Analogous to potential, but for the radial acceleration.
!   useful for coding initial condition with centrifugal balance
!
!  21-aug-07/wlad: coded
!
      use Mpicomm,only: stop_it
      use Sub,    only: get_radial_distance
!
      real, dimension (:) :: g_r
      real, dimension(size(g_r)) :: rr_mn,rr_sph,rr_cyl,pot
      integer :: j
!
      intent(out) :: g_r
!
      call get_radial_distance(rr_sph,rr_cyl)
      if (lcylindrical_gravity) then
        rr_mn=rr_cyl
      else
        rr_mn=rr_sph
      endif
!
      g_r=0.
      do j=1,ninit
        select case (ipotential(j))
        case ('no-smooth','newton','newtonian')
          g_r=g_r -g0/rr_mn**2
!
        case ('smoothed-newton')
          !g_r=g_r -g0*rr_mn**(n_pot-1) &
          !     *(rr_mn**n_pot+r0_pot**n_pot)**(-1./n_pot-1.)
          pot = -g0*(1. + (r1_pot1*rr_mn)**n_pot1)**(1.0/n_pot1) &
                   /(rr_mn**n_pot + r0_pot**n_pot)**(1.0/n_pot)
          g_r = pot*(rr_mn**(n_pot-1)/(rr_mn**n_pot + r0_pot**n_pot) &
                   - r1_pot1*(r1_pot1*rr_mn)**(n_pot1-1) &
                            /(1 + (r1_pot1*rr_mn)**n_pot1))
!
        case ('varying-q')
          g_r=g_r -g0/rr_mn**(2*qgshear-1)
!
        case ('varying-q-smooth')
          g_r=g_r -g0*rr_mn/(rr_mn**2+r0_pot**2)**qgshear
!
        case ('dark-matter-halo')
          g_r=g_r -g01(j)*(1-rpot(j)/rr_mn*atan2(rr_mn,rpot(j)))/rr_mn
!
        case ('light-matter')
          g_r=g_r-4.134e-4*g01(j)*(.5*rr_mn/rpot(j))**1.22/&
               ((.5*rr_mn/rpot(j))**2+1.502)**1.43/rr_mn
!
        case ('zero')
          g_r=g_r
!
        case default
          if (lroot) print*, 'acceleration: '//&
               'No such value for ipotential: ', trim(ipotential(j))
          call stop_it("")
!
        endselect
      enddo
!
    endsubroutine acceleration_penc_1D
!***********************************************************************
    subroutine acceleration_point(x,y,z,r,g_r)
!
!  Gravitational acceleration in one point
!
!  Analogous to potential, but for the radial acceleration.
!   useful for coding initial condition with centrifugal balance
!
!  18-nov-08/wlad: coded
!
      use Mpicomm,only: stop_it
!
      real :: g_r,rad,pot
      real, optional :: x,y,z,r
      integer :: j
!
      intent(in)  :: x,y,z,r
      intent(out) :: g_r
!
      if (present(r)) then
        rad = r
      else
        if (present(x) .and. present(y) .and. present(z)) then
!
          if (.not.lcartesian_coords) &
              call stop_it("gravity_r: acceleration_point with x,y,z is "//&
              "not yet implemented for non-cartesiand coordinates. Fix  "//&
              "the call to  use radial distance instead")
!
          rad = sqrt(x**2+y**2+z**2)
        else
          call stop_it("Need to specify either x,y,z or r in acceleration_point()")
        endif
      endif
!
      g_r=0.
      do j=1,ninit
        select case (ipotential(j))
        case ('no-smooth','newton','newtonian')
          g_r=g_r -g0/rad**2
!
        case ('smoothed-newton')
          !g_r=g_r -g0*rad**(n_pot-1) &
          !     *(rad**n_pot+r0_pot**n_pot)**(-1./n_pot-1.)
          pot = -g0*(1. + (r1_pot1*rad)**n_pot1)**(1.0/n_pot1) &
                   /(rad**n_pot + r0_pot**n_pot)**(1.0/n_pot)
          g_r = pot*(rad**(n_pot-1)/(rad**n_pot + r0_pot**n_pot) &
                   - r1_pot1*(r1_pot1*rad)**(n_pot1-1) &
                            /(1 + (r1_pot1*rad)**n_pot1))
!
        case ('varying-q')
          g_r=g_r -g0/rad**(2*qgshear-1)
!
        case ('varying-q-smooth')
          g_r=g_r -g0*rad/(rad**2+r0_pot**2)**qgshear
!
        case ('dark-matter-halo')
          g_r=g_r -g01(j)*(1-rpot(j)/rad*atan2(rad,rpot(j)))/rad
!
        case ('light-matter')
          g_r=g_r-4.134e-4*g01(j)*(.5*rad/rpot(j))**1.22/&
               ((.5*rad/rpot(j))**2+1.502)**1.43/rad
!
        case ('zero')
          g_r=g_r
!
        case default
          if (lroot) print*, 'acceleration: '//&
               'No such value for ipotential: ', trim(ipotential(j))
          call stop_it("")
!
        endselect
      enddo
!
    endsubroutine acceleration_point
!***********************************************************************
    subroutine get_gravity_field(gr,gg_mn,rr_mn)
!
!  Calculate gravity field for different coordinate systems
!
!  15-mar-07/wlad: coded
!
      real, dimension(nx),intent(in) :: gr,rr_mn
      real, dimension(nx,3),intent(out) :: gg_mn
      real, dimension(nx,3) :: ggp
!
      if (coord_system=='cartesian') then
        gg_mn(:,1) = x(l1:l2)/rr_mn*gr
        gg_mn(:,2) = y(  m  )/rr_mn*gr
        gg_mn(:,3) = z(  n  )/rr_mn*gr
        if (lcylindrical_gravity) gg_mn(:,3)=0.
      elseif (coord_system=='cylindric') then
        gg_mn(:,1) = x(l1:l2)/rr_mn*gr
        gg_mn(:,2) = 0.
        gg_mn(:,3) = z(  n  )/rr_mn*gr
        if (lcylindrical_gravity) gg_mn(:,3)=0.
      elseif (coord_system=='spherical') then
        gg_mn(:,2)=0.
        gg_mn(:,3)=0.
        if (lcylindrical_gravity) then
          gg_mn(:,1) = gr*sinth(m)
          gg_mn(:,2) = gr*costh(m)
        else
          gg_mn(:,1) = gr
        endif
      endif
!
!  Add the gravity of a stationary secondary body (i.e., following reference frame)
!
      if (g1/=0 .and. (.not.(lramp_mass.or.lsecondary_wait))) then 
!
!  In this case, the mass is constant in time. Add to the global array.
!
        call secondary_body_gravity(ggp)
        gg_mn = gg_mn + ggp 
      endif
!
    endsubroutine get_gravity_field
!***********************************************************************
    subroutine secondary_body_gravity(ggp)
!
!  Add the gravity of a body fixed at position (rp1,0,0).
!
!  20-jul-14/wlad: coded
!  03-may-16/wlad: added boley smoothing      
!
      real, dimension(nx,3), intent(out) :: ggp
      real, dimension(nx) :: rr2_pm,gp
      real :: g2
      integer :: i
!
      if (lcylindrical_coords) then
        if (lcylindrical_gravity) then
          rr2_pm = x(l1:l2)**2 + rp1**2 -2*x(l1:l2)*rp1*cos(y(m))
        else
          rr2_pm = x(l1:l2)**2 + rp1**2 -2*x(l1:l2)*rp1*cos(y(m)) + z(n)**2
        endif
      elseif (lspherical_coords) then
        rr2_pm = x(l1:l2)**2 + rp1**2 - 2*x(l1:l2)*rp1*sinth(m)*cosph(n)
      else
        call fatal_error("secondary_body_gravity","not coded for Cartesian")
      endif
!
!  Select the potential smoothing for the secondary. 
!      
      g2=g1
      if (lramp_mass) call rampup_secondary_mass(g2)
!
      select case (ipotential_secondary)
!
      case ('plummer')
        gp = -g2*(rr2_pm+rp1_smooth**2)**(-1.5)
!           
      case ('boley')
!
!  Correct potential outside a sphere of radius rsmooth. Default to Hill sphere.
!
        do i=1,nx
          if (rr2_pm(i) .gt. rp1_smooth**2) then
            gp(i) = -g2*rr2_pm(i)**(-1.5)
          else
            gp(i) =  g2*(3*sqrt(rr2_pm(i))*rp1_smooth1 - 4)*rp1_smooth1**3
          endif
        enddo
!           
      case default
!
!  Catch unknown values
!           
        if (lroot) print*, 'secondary_body_gravity: '//&
             'No such value for ipotential: ', trim(ipotential_secondary)
        call fatal_error("","")
!
      endselect
!
      if (lsecondary_wait.and.(t<=t_start_secondary)) gp = 0.
!       
!  Set the acceleration
!
      if (lcylindrical_coords) then 
        ggp(:,1) =  gp * (x(l1:l2)-rp1*cos(y(m)))
        ggp(:,2) =  gp *           rp1*sin(y(m))
        ggp(:,3) =  gp *  z(  n  )
        if (lcylindrical_gravity) ggp(:,3)=0.
      else if (lspherical_coords) then
        ggp(:,1) =  gp * (x(l1:l2) - rp1*sinth(m)*cosph(n))
        ggp(:,2) = -gp *             rp1*costh(m)*cosph(n)
        ggp(:,3) =  gp *             rp1*         sinph(n)
      endif
!
    endsubroutine secondary_body_gravity
!***********************************************************************
    subroutine calc_torque(p)
!
!  Output torque diagnostic for nbody particle ks
!
!  24-jan-20/wlad : coded
!
      use Diagnostics
!
      type (pencil_case) :: p
      real, dimension(nx) :: torque
      real, dimension(nx) :: rr2,rpre
!
      if (lcartesian_coords) then
        call fatal_error("calc_torque","not coded for Cartesian")
      elseif (lcylindrical_coords) then
        rpre = x(l1:l2)*rp1*sin(y(m))
        rr2  = x(l1:l2)**2 + rp1**2 -2*x(l1:l2)*rp1*cos(y(m)) + z(n)**2
      elseif (lspherical_coords) then
        rpre = x(l1:l2)*rp1*sinth(m)*sinph(n)
        rr2  = x(l1:l2)**2 + rp1**2 - 2*x(l1:l2)*rp1*sinth(m)*cosph(n)
      else
        call fatal_error("calc_torque",&
             "the world is flat and we should never gotten here")
      endif
!
!  Define separate torques for gas and dust/particles
!
      torque = g1*p%rho*rpre*(rr2 + rp1_smooth**2)**(-1.5)
!
      call integrate_mn_name(torque,idiag_torque)
!
    endsubroutine calc_torque
!***********************************************************************
    subroutine rprint_gravity(lreset,lwrite)
!
!  reads and registers print parameters relevant for gravity advance
!  dummy routine
!
!  26-apr-03/axel: coded
!
      use FArrayManager, only: farray_index_append
      use Diagnostics
!
      logical :: lreset,lwr
      logical, optional :: lwrite
      integer :: iname
!
      lwr = .false.
      if (present(lwrite)) lwr=lwrite
!
      if (lreset) then
        idiag_torque=0
      endif

      do iname=1,nname
        call parse_name(iname,cname(iname),cform(iname),'torque',idiag_torque)
      enddo
!
      call keep_compiler_quiet(lreset)
!
    endsubroutine rprint_gravity
!***********************************************************************
    subroutine compute_gravity_star(f, wheat, luminosity, star_cte)
!
!  Compute the global gravity field for a star-in-a-box with a central heating
!  05-jan-2010/dintrans: coded
!
    use Sub, only: get_radial_distance, erfunc
!
    real, dimension (mx,my,mz,mfarray) :: f
    real :: wheat, luminosity, star_cte
    real, dimension (nx) :: g_r, rr_mn, rr_sph, rr_cyl, u_mn, lumi_mn
    real, dimension (nx,3) :: gg_mn=0.0
!
    do n=n1,n2
    do m=m1,m2
      call get_radial_distance(rr_sph, rr_cyl)
      rr_mn=rr_sph
      u_mn=rr_mn/sqrt(2.)/wheat
      if (nzgrid==1) then
        lumi_mn=luminosity*(1.-exp(-u_mn**2))
        g_r=-lumi_mn/(2.*pi*rr_mn)*star_cte
      else
        lumi_mn=luminosity*(erfunc(u_mn)-2.*u_mn/sqrt(pi)*exp(-u_mn**2))
        g_r=-lumi_mn/(4.*pi*rr_mn**2)*star_cte
      endif
      call get_gravity_field(g_r, gg_mn, rr_mn)
      f(l1:l2,m,n,iglobal_gg:iglobal_gg+2) = gg_mn
    enddo
    enddo
!
    endsubroutine compute_gravity_star
!***********************************************************************
    subroutine get_xgravity(xgrav)
!
!  Dummy routine for getting the gravity profile into the
!  intial conditions
!
!  04-oct-10/bing: coded
!
      real, dimension(nx) :: xgrav
!
      call keep_compiler_quiet(xgrav)
!
    endsubroutine get_xgravity
!***********************************************************************
    subroutine set_consistent_gravity(ginput,gtype,gprofile,lsuccess)
!
! Dummy routine
!
      real :: ginput
      character (len=labellen) :: gtype,gprofile
      logical :: lsuccess
!
      call keep_compiler_quiet(ginput)
      call keep_compiler_quiet(gtype,gprofile)
!
! This routine should never be called in the way it is written now.
!
      lsuccess=.false.
!
! gravity parameters set consistently.
!
    endsubroutine set_consistent_gravity
!***********************************************************************
    logical function is_constant_zgrav()
!
!  14-apr-15/MR: coded
!
      is_constant_zgrav=.false.

    endfunction is_constant_zgrav
!***********************************************************************
endmodule Gravity