! $Id$
!
!  This module takes care of simple types of gravity, i.e. where
!    gx=gx(x) or gy=gy(y) or gz=gz(z)
!  Here the gravity master pencils gravx_xpencil, gravy_ypencil and
!  gravz_zpencil only need to be calculated once, and then these can
!  simply be added to the equations of motion again and again.
!
!** 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
!
! PENCILS PROVIDED gg(3); epot
!
!***************************************************************
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
!
  double precision, parameter :: g_B_cgs=6.172d20, g_D_cgs=3.086d21
  double precision :: g_B, g_D, g_B_factor=1.0, g_D_factor=1.0
  real :: gravitational_const=0., mass_cent_body=0.
  real, dimension(mx) :: gravx_xpencil=0.0, potx_xpencil=0.0
  real, dimension(my) :: gravy_ypencil=0.0, poty_ypencil=0.0
  real, dimension(mz) :: gravz_zpencil=0.0, potz_zpencil=0.0
  real, dimension(mx) :: xdep=0.0
  real, dimension(mz) :: zdep=0.0
  real, parameter :: g_A_cgs=4.4e-9, g_C_cgs=1.7e-9
  real :: gravx=0.0, gravy=0.0, gravz=0.0
  real :: kx_gg=1.0, ky_gg=1.0, kz_gg=1.0, gravz_const=1.0, reduced_top=1.0
  real :: xgrav=impossible, ygrav=impossible, zgrav=impossible
  real :: xinfty=0.0, yinfty=0.0, zinfty=impossible, zclip=impossible
  real :: dgravx=0.0, pot_ratio=1.0
  real :: z1=0.0, z2=1.0, xref=0.0, zref=impossible, sphere_rad=0.0, g_ref=0.0
  real :: nu_epicycle=1.0, nu_epicycle2=1.0
  real :: nux_epicycle=0.0, nux_epicycle2=0.0
  real :: g_A, g_C, g_A_factor=1.0, g_C_factor=1.0
  real :: cs0hs=0.0, H0hs=0.0
  real :: potx_const=0.0, poty_const=0.0, potz_const=0.0
  integer :: n_pot=10
  integer :: n_adjust_sphersym=0
  character (len=labellen) :: gravx_profile='zero',gravy_profile='zero', &
                              gravz_profile='zero'
!
!  Parameters used by other modules (only defined for other gravities)
!
  logical :: lnumerical_equilibrium=.false.
  logical :: lxyzdependence=.false.
  logical :: lcalc_zinfty=.false.
  logical :: lboussinesq_grav=.false.
  logical :: ladjust_sphersym=.false.
 
  real :: g0=0.0
  real :: lnrho_bot=0.0, lnrho_top=0.0, ss_bot=0.0, ss_top=0.0
  real :: kappa_x1=0.0, kappa_x2=0.0, kappa_z1=0.0, kappa_z2=0.0
  real :: grav_tilt=0.0, grav_amp=0.0
!
  namelist /grav_init_pars/ &
      gravx_profile, gravy_profile, gravz_profile, gravx, gravy, gravz, &
      xgrav, ygrav, zgrav, kx_gg, ky_gg, kz_gg, dgravx, pot_ratio, z1, z2, &
      nux_epicycle, nu_epicycle, xref, zref, g_ref, sphere_rad, lnrho_bot, lnrho_top, &
      ss_bot, ss_top, lgravx_gas, lgravx_dust, lgravy_gas, lgravy_dust, &
      lgravz_gas, lgravz_dust, xinfty, yinfty, zinfty, lxyzdependence, &
      lcalc_zinfty, kappa_x1, kappa_x2, kappa_z1, kappa_z2, reduced_top, &
      lboussinesq_grav, n_pot, cs0hs, H0hs, grav_tilt, grav_amp, &
      potx_const,poty_const,potz_const, zclip, n_adjust_sphersym, gravitational_const, &
      mass_cent_body, g_A_factor, g_C_factor, g_B_factor, g_D_factor
!
  namelist /grav_run_pars/ &
      gravx_profile, gravy_profile, gravz_profile, gravx, gravy, gravz, &
      xgrav, ygrav, zgrav, kx_gg, ky_gg, kz_gg, dgravx, pot_ratio, &
      nux_epicycle, nu_epicycle, xref, zref, g_ref, sphere_rad, &
      lgravx_gas, lgravx_dust, lgravy_gas, lgravy_dust, &
      lgravz_gas, lgravz_dust, xinfty, yinfty, zinfty, lxyzdependence, &
      lcalc_zinfty, kappa_x1, kappa_x2, kappa_z1, kappa_z2, reduced_top, &
      lboussinesq_grav, n_pot, grav_tilt, grav_amp, &
      potx_const,poty_const,potz_const, zclip, n_adjust_sphersym, gravitational_const, &
      mass_cent_body, g_A_factor, g_C_factor, g_B_factor, g_D_factor
!
!  Diagnostic variables for print.in
! (needs to be consistent with reset list below)
!
  integer :: idiag_epot=0          ! DIAG_DOC: $\left<\varrho \Phi_{\rm grav}
                                   ! DIAG_DOC: \right>$ \quad(mean potential
                                   ! DIAG_DOC: energy)
  integer :: idiag_epottot=0       ! DIAG_DOC: $\int_V\varrho \Phi_{\rm grav}
                                   ! DIAG_DOC: dV$ \quad(total potential
                                   ! DIAG_DOC: energy)
  integer :: idiag_ugm=0           ! DIAG_DOC: $\left<\uv \cdot \gv\right>$
!
! xy averaged diagnostics given in xyaver.in written every it1d timestep
!
  integer :: idiag_epotmz=0        ! XYAVG_DOC: $\left<\varrho \Phi_{\rm grav}
                                   ! XYAVG_DOC: \right>_{xy}$
  integer :: idiag_epotuzmz=0      ! XYAVG_DOC: $\left<\varrho \Phi_{\rm grav}
                                   ! XYAVG_DOC: u_z \right>_{xy}$
                                   ! XYAVG_DOC: \quad(potential energy flux)
!
! xz averaged diagnostics given in xzaver.in
!
  integer :: idiag_epotmy=0        ! XZAVG_DOC: $\left<\varrho \Phi_{\rm grav}
                                   ! XZAVG_DOC: \right>_{xz}$
!
! yz averaged diagnostics given in yzaver.in
!
  integer :: idiag_epotmx=0        ! YZAVG_DOC: $\left<\varrho \Phi_{\rm grav}
                                   ! YZAVG_DOC: \right>_{yz}$
  integer :: idiag_epotuxmx=0      ! YZAVG_DOC: $\left<\varrho \Phi_{\rm grav}
                                   ! YZAVG_DOC: u_x \right>_{yz}$
                                   ! YZAVG_DOC: \quad(potential energy flux)
!
! z averaged diagnostics given in zaver.in
!
  integer :: idiag_epotmxy=0       ! ZAVG_DOC: $\left<\varrho \Phi_{\rm grav}
                                   ! ZAVG_DOC: \right>_{z}$
  integer :: idiag_epotuxmxy=0     ! ZAVG_DOC: $\left<\varrho \Phi_{\rm grav}
                                   ! ZAVG_DOC: u_x \right>_{z}$
                                   ! ZAVG_DOC: \quad(potential energy flux)
!
! work variables
!
  real, dimension(mx) :: gravx_xpencil_0
  real ::G4pi
!  
  contains
!***********************************************************************
    subroutine register_gravity()
!
!  Initialise gravity variables (currently none).
!
!  12-nov-04/anders: coded
!
!  Identify version number (generated automatically by SVN).
!
      if (lroot) call svn_id( &
          '$Id$')
!
    endsubroutine register_gravity
!***********************************************************************
    subroutine initialize_gravity(f)
!
!  Calculate master pencils for gravity. These are put into gravity pencils
!  in the subroutine calc_pencils_grav.
!
!  12-nov-04/anders: coded, copied init conds from grav_x, grav_y and grav_y.
!   9-jun-15/MR: added parameter n_adjust_sphersym: if > 0 after each n_adjust_sphersym
!                timesteps the spherically symmetric part of gravity is adjusted
!                according to the actual density distribution.
!                Only in effect for spherical co-ordinates.
!  12-jun-15/MR: added (alternative) parameters gravitational_const and mass_cent_body for 
!                gravity adjustment.
!                For Kepler profile, gravitational_const is calculated from gravx and mass_cent_body.
!
      use SharedVariables, only: put_shared_variable, get_shared_variable
      use General, only: notanumber
      use Sub, only: cubic_step
      use Mpicomm, only: mpibcast_real, MPI_COMM_WORLD
!
      real, dimension(mx,my,mz,mfarray) :: f
      real, pointer :: cs20,mpoly,gamma
!
      real, dimension (mz) :: prof
      real :: ztop
!
!  Sanity check.
!
      if (lrun) then
        if (gravx_profile == 'zero' .and. &
            gravy_profile == 'zero' .and. &
            gravz_profile == 'zero') then
          call warning('initialize_gravity', &
              'You do not need gravity_simple for zero gravity...')
        endif
      endif
!
!  Possibility of specifying zref (if zinfty is not given).
!  In that case we need to know cs20 and mpoly from the EOS module.
!  We do this using shared variables.
!  Normally we give zinfty, but if this is not the case, we set it to zero.
!  zinfty has to be higher than the top of a polytropic atmosphere.
!
      if (lcalc_zinfty) then
        if (zinfty==impossible) then
          if (zref==impossible) then
            call fatal_error('initialize_gravity','zref=impossible')
          else
            call get_shared_variable('cs20',cs20,caller='initialize_gravity')
            call get_shared_variable('gamma',gamma)
            if (ldensity.and..not.lstratz) then
              call get_shared_variable('mpoly',mpoly)
            else
              if (lroot) call warning('initialize_eos','mpoly not obtained from density,'// &
                                      'set impossible')
              allocate(mpoly); mpoly=impossible
            endif
!
            zinfty=zref+cs20*(mpoly+1)/(-gamma*gravz)
            if (lroot) print*,'initialize_gravity: computed zinfty=',zinfty
          endif
        endif
      else
        if (zinfty==impossible) zinfty=0.
      endif
!
!  Different x-gravity profiles.
!
      if (gravx/=0) lgravx=.true.
!
      select case (gravx_profile)
!
      case ('zero')
        if (lroot) print*,'initialize_gravity: no x-gravity'
        gravx_xpencil=0.
        lgravx_gas=.false.
        lgravx_dust=.false.
!
      case ('const')
        if (lroot) print*,'initialize_gravity: constant x-grav=', gravx
        gravx_xpencil=gravx
        potx_xpencil=-gravx*(x-xinfty)
        call put_shared_variable('gravx', gravx, caller='initialize_gravity')
        call put_shared_variable('gravx_xpencil', gravx_xpencil)
!
      case ('const_tilt')
        gravx=grav_amp*sin(grav_tilt*pi/180.)
        if (lroot) print*,'initialize_gravity: constant gravx=', gravx
        gravx_xpencil=gravx
        potx_xpencil=-gravx*(x-xinfty)
        call put_shared_variable('gravx', gravx, caller='initialize_gravity')
!
! Linear gravity profile in x for planetary core dynamos
!
      case ('linear_x')
        if (lroot) print*,'initialize_gravity: linear x-grav, gravx=', gravx
        gravx_xpencil = -gravx*x
        potx_xpencil  = 0.5*gravx*(x**2-xinfty**2)
        call put_shared_variable('gravx', gravx, caller='initialize_gravity')
!
!  Linear gravity potential with additional z dependence.
!  Calculate zdep here, but don't multiply it onto gravx_xpencil
!  or potx_xpencil, so they will not be correct yet!
!
      case ('linear_zdep')      !  Linear gravity profile (for accretion discs)
        nux_epicycle2=nux_epicycle**2
        if (lroot) print*,'initialize_gravity: linear x-grav with z dep, nu=', &
          nux_epicycle,kappa_z1,kappa_z2
        zdep=(1.+kappa_z1*z+.5*(kappa_z1*z)**2)
        gravx_xpencil = -nux_epicycle2*x
        potx_xpencil=0.5*nux_epicycle2*(x**2-xinfty**2)
!
!  tanh profile
!  For isothermal EOS, we have 0=-cs2*dlnrho+gravx.
!  pot_ratio gives the resulting ratio in the density.
!AB: these potx_const values seem to be missing in corresponding entries for potx_xpoint.
!
      case ('tanh-pot')
        if (dgravx==0.) call fatal_error("initialize_gravity","dgravx=0 not OK")
        if (lroot) print*,'initialize_gravity: tanh x-grav, gravx=',gravx
        if (lroot) print*,'initialize_gravity: xgrav,dgravx=',xgrav,dgravx
        gravx=-alog(pot_ratio)/dgravx
        gravx_xpencil=gravx*.5/cosh((x-xgrav)/dgravx)**2
        potx_xpencil=-gravx*.5*(1.+tanh((x-xgrav)/dgravx))*dgravx + potx_const
!
      case ('sinusoidal')
        if (lroot) print*,'initialize_gravity: sinusoidal x-grav, gravx=',gravx
        gravx_xpencil = -gravx*sin(kx_gg*x)
        potx_xpencil  = -gravx/kx_gg*cos(kx_gg*x) + potx_const
!
      case ('Baker74')
        !abag added, need to make adaptable to any width/position
        if (lroot) print*,'initialize_gravity: Baker_74=',gravx
        gravx_xpencil =(tanh((x+pi/3.)/0.1)+tanh(-(x-pi/3.)/0.1))/2.*&
            gravx*sin(2*(x-pi/2.))
        potx_xpencil=(tanh((x+pi/3.)/0.1)+tanh(-(x-pi/3.)/0.1))/2.*&
            gravx*(.5*cos(2*(x-pi/2.))-0.5) + potx_const
!
      case ('kepler')
        if (lroot) print*,'initialize_gravity: kepler x-grav, gravx=',gravx
        gravx_xpencil=-gravx/x**2
        potx_xpencil=-gravx/x + potx_const
        g0=gravx
        call put_shared_variable('gravx', gravx, caller='initialize_gravity')
        call put_shared_variable('gravx_xpencil', gravx_xpencil)
!
      case ('kepler_2d')
        if (lroot) print*,'initialize_gravity: kepler_2d x-grav, gravx=',gravx
        gravx_xpencil=-gravx/x
        potx_xpencil=-gravx*alog(x) + potx_const
        g0=gravx
        call put_shared_variable('gravx', gravx, caller='initialize_gravity')
        call put_shared_variable('gravx_xpencil', gravx_xpencil)
!
!  Convection zone model, normalized to the bottom of the domain
!
      case ('CZbot1')
        if (lroot) print*,'initialize_gravity: CZbot1 x-grav, gravx=',gravx
        gravx_xpencil=-gravx/x**2
        potx_xpencil=-gravx*(1./x-1./xyz0(1)) + potx_const
        g0=gravx
        call put_shared_variable('gravx', gravx, caller='initialize_gravity')
        call put_shared_variable('gravx_xpencil', gravx_xpencil)
!
!  Convection zone model, normalized to the middle of the domain
!
      case ('CZmid1')
        if (lroot) print*,'initialize_gravity: CZmid1 x-grav, gravx=',gravx
        gravx_xpencil=-gravx/x**2
        potx_xpencil=-gravx*(1./x-2./(xyz0(1)+xyz1(1))) + potx_const
        g0=gravx
        call put_shared_variable('gravx', gravx, caller='initialize_gravity')
        call put_shared_variable('gravx_xpencil', gravx_xpencil)
!
      ! solid sphere (homogenous) gravity profile
      ! 'xref' determines the location of the sphere border relative to the lower box boundary.
      ! 'g_ref' is the gravity acceleration at xref (implies the mass of the sphere).
      ! 'sphere_rad' is the radius of the sphere (independant from box coordinates).
      ! Together, sphere_rad and xref describe, how much of the sphere lies inside the box.
      case ('solid_sphere')
        if (lroot) print *, 'initialize_gravity: solid sphere xref=', xref, ', g_ref=', g_ref, ', sphere_rad=', sphere_rad
        where (x > xref)
          gravx_xpencil = g_ref * sphere_rad**2 / (x-xref+sphere_rad)**2
        else where (x < xref-2*sphere_rad)
          gravx_xpencil = -g_ref * sphere_rad**2 / (x-xref+sphere_rad)**2
        else where
          gravx_xpencil = g_ref * (x-xref+sphere_rad) / sphere_rad
        end where
        potx_xpencil=-gravx_xpencil*(x-xinfty)
!
      case ('loop')
        if (lroot) print*,'initialize_gravity: 1D loop, gravx=',gravx
        gravx_xpencil=cos(x / (2*xyz0(1)+Lxyz(1)) * pi)
        gravx_xpencil=gravx * gravx_xpencil
        potx_xpencil =sin(x / (2*xyz0(1)+Lxyz(1)) * pi)
        potx_xpencil =gravx*(2*xyz0(1)+Lxyz(1))/pi * potx_xpencil + potx_const
!
      case ('half-loop')
        if (lroot) print*,'initialize_gravity: 1D half-loop, gravx=',gravx
        gravx_xpencil=cos(x / (2*xyz0(1)+2*Lxyz(1)) * pi)
        gravx_xpencil=gravx * gravx_xpencil
        potx_xpencil =sin(x / (2*xyz0(1)+2*Lxyz(1)) * pi)
        potx_xpencil =gravx*(2*xyz0(1)+2*Lxyz(1))/pi * potx_xpencil + potx_const
!
      case default
        if (lroot) print*, &
            'initialize_gravity: unknown gravx_profile ', gravx_profile
        call fatal_error('initialize_gravity','chosen gravx_profile not valid')
!
      endselect
!
!  Different y-gravity profiles.
!
      select case (gravy_profile)
!
      case ('zero')
        if (lroot) print*,'initialize_gravity: no y-gravity'
        gravy_ypencil=0.
        lgravy_gas=.false.
        lgravy_dust=.false.
!
      case ('const')
        if (lroot) print*,'initialize_gravity: constant y-grav=', gravy
        gravy_ypencil = gravy
        poty_ypencil  = -gravy*(y-yinfty)
!
      case ('sinusoidal')
        if (lroot) print*,'initialize_gravity: sinusoidal y-grav, gravy=', gravy
        gravy_ypencil = -gravy*sin(ky_gg*y)
        poty_ypencil  = -gravy/ky_gg*cos(ky_gg*y) + poty_const
!
      case ('kepler')
        if (lroot) print*,'initialize_gravity: kepler gravy, y-grav=', gravy
        gravy_ypencil = -gravy/y**2
        poty_ypencil  = -gravy/y + poty_const
!
      case default
        if (lroot) print*, &
            'initialize_gravity: unknown gravy_profile ', gravy_profile
        call fatal_error('initialize_gravity','chosen gravy_profile not valid')
!
      endselect
!
!  Different z-gravity profiles.
!  Set lgravz=T only when gravz_profile is not zero.
!
      if (gravz_profile/='zero') lgravz=.true.
!
      select case (gravz_profile)
!
      case ('zero')
        if (lroot) print*,'initialize_gravity: no z-gravity'
        gravz_zpencil=0.
        lgravz_gas=.false.
        lgravz_dust=.false.
!
      case ('const')
        if (lroot) print*,'initialize_gravity: constant gravz=', gravz
        gravz_zpencil=gravz
        if (zclip==impossible) then
          potz_zpencil=-gravz*(z-zinfty)
        else
          potz_zpencil=-gravz*min(z-zinfty,zclip-zinfty)
        endif
!
      case ('const_tilt')
        gravz=-grav_amp*cos(grav_tilt*pi/180.)
        if (lroot) print*,'initialize_gravity: constant gravz=', gravz
        gravz_zpencil=gravz
        potz_zpencil=-gravz*(z-zinfty)
!
      case ('tanh')
        if (lroot) print*,'initialize_gravity: tanh gravz=', gravz
        gravz_zpencil=-gravz*tanh(z/zref)
        potz_zpencil=gravz*zref*alog(cosh(z/zref))
!
      case ('boussinesq')
        if (lroot) print*,'initialize_gravity: boussinesq gravz=', gravz
        gravz_zpencil=gravz
        potz_zpencil=-gravz*(z-zinfty)
        lboussinesq_grav=.true.
!
      case ('const_zero')  !  Const. gravity acc. (but zero for z>zgrav)
        if (headtt) print*,'initialize_gravity: const_zero gravz=', gravz
        if (zgrav==impossible .and. lroot) &
            print*,'initialize_gravity: zgrav is not set!'
        do n=n1,n2
          if (z(n)<=zgrav) gravz_zpencil(n) = gravz
        enddo
!
      case ('linear')      !  Linear gravity profile (for accretion discs)
        nu_epicycle2=nu_epicycle**2
        if (lroot) print*,'initialize_gravity: linear z-grav, nu=', nu_epicycle
        gravz_zpencil=-nu_epicycle2*z
        potz_zpencil=0.5*nu_epicycle2*(z**2-zinfty**2)
!
      ! solid sphere (homogenous) gravity profile
      ! 'zref' determines the location of the sphere border relative to the lower box boundary.
      ! 'g_ref' is the gravity acceleration at zref (implies the mass of the sphere).
      ! 'sphere_rad' is the radius of the sphere (independant from box coordinates).
      ! Together, sphere_rad and zref describe, how much of the sphere lies inside the box.
      case ('solid_sphere')
        if (lroot) print *, 'initialize_gravity: solid sphere zref=', zref, ', g_ref=', g_ref, ', sphere_rad=', sphere_rad
        where (z > zref)
          gravz_zpencil = g_ref * sphere_rad**2 / (z-zref+sphere_rad)**2
        else where (z < zref-2*sphere_rad)
          gravz_zpencil = -g_ref * sphere_rad**2 / (z-zref+sphere_rad)**2
        else where
          gravz_zpencil = g_ref * (z-zref+sphere_rad) / sphere_rad
        end where
        potz_zpencil = -gravz_zpencil * (z - zinfty)
!
      case ('spherical')
        nu_epicycle2=nu_epicycle**2
        if (lroot) print*,'initialize_gravity: spherical z-grav, nu, z1=', nu_epicycle, z1
        gravz_zpencil=-nu_epicycle2*z/(1.0+(z/z1)**2)
        potz_zpencil=0.5*nu_epicycle2*z1**2*alog(1.0+(z/z1)**2) + potz_const
!
!  Linear gravity potential with additional x dependence.
!  Calculate xdep here, but don't multiply it onto gravz_zpencil
!  or potz_zpencil, so they will not be correct yet!
!
      case ('linear_xdep')
        nu_epicycle2=nu_epicycle**2
        if (lroot) print*,'initialize_gravity: linear z-grav with x dep, nu=', &
          nu_epicycle,kappa_x1,kappa_x2
        xdep=(1.+kappa_x1*x+.5*(kappa_x1*x)**2)
        gravz_zpencil = -nu_epicycle2*z
        potz_zpencil=0.5*nu_epicycle2*(z**2-zinfty**2)
!
      case ('linear_smoothed')
        nu_epicycle2=nu_epicycle**2
        if (lroot) print*,'initialize_gravity: linear z-grav, '// &
                         'smoothed to zero at top/bottom, nu=', nu_epicycle
        prof = 1. + (z/zref)**(2*n_pot)
        gravz_zpencil = -nu_epicycle2*z/prof**(1./n_pot+1.)
        potz_zpencil = 0.5*nu_epicycle2*z**2/prof**(1./n_pot) + potz_const
!
      case ('sinusoidal')
        if (lroot) print*,'initialize_gravity: sinusoidal z-grav, gravz=', gravz
        gravz_zpencil = -gravz*sin(kz_gg*z)
        potz_zpencil = -gravz/kz_gg*cos(kz_gg*z) + potz_const
!
      case ('kepler')
        if (lroot) print*,'initialize_gravity: kepler z-grav, gravz=', gravz
        gravz_zpencil = -gravz/z**2
        potz_zpencil  = -gravz/z + potz_const
!
      case ('Ferriere')
!
!  Set up physical units.
!
      if (unit_system=='cgs') then
        g_A = g_A_factor*g_A_cgs/unit_velocity*unit_time
        g_B = g_B_factor*g_B_cgs/unit_length
        g_C = g_C_factor*g_C_cgs/unit_velocity*unit_time
        g_D = g_D_factor*g_D_cgs/unit_length
      else if (unit_system=='SI') then
        call fatal_error('initialize_gravity','SI unit conversions not inplemented')
      endif
!
!  Gravity profile from K. Ferriere, ApJ 497, 759, 1998, eq (34)
!  at solar radius.  (for interstellar runs)
!
!  nb: 331.5 is conversion factor: 10^-9 cm/s^2 -> kpc/Gyr^2)  (/= 321.1 ?!?)
!AB: These numbers should be inserted in the appropriate units.
!AB: As it is now, it can never make much sense.
        gravz_zpencil = -(g_A*z/sqrt(z**2+g_B**2) + g_C*z/g_D)
!
      case ('Galactic-hs')
        if (lroot) print*,'Galactic hydrostatic equilibrium gravity profile'
        if (lroot.and.(cs0hs==0.or.H0hs==0)) &
            call fatal_error('initialize-gravity', &
            'Set cs0hs and H0hs in grav_init_pars!')
        gravz_zpencil = -z*(cs0hs/H0hs)**2/sqrt(1 + (z/H0hs)**2)
!
      case ('reduced_top')
        if (lroot) print*,'initialize_gravity: reduced, gravz=',gravz
        if (zgrav==impossible.and.lroot) print*,'zgrav is not set!'
        ztop = xyz0(3)+Lxyz(3)
        prof = cubic_step(z,(zgrav+ztop)/2,(ztop-zgrav)/2)
        gravz_zpencil = (1 - prof*(1-reduced_top))*gravz
!
      case default
        if (lroot) print*, &
            'initialize_gravity: unknown gravz_profile ', gravz_profile
        call fatal_error('initialize_gravity','chosen gravz_profile not valid')
!
      endselect
!
!  Sanity check.
!
      if (notanumber(gravx_xpencil)) &
        call fatal_error('initialize_gravity','found NaN or +/-Inf in gravx_xpencil')
      if (notanumber(gravy_ypencil)) &
        call fatal_error('initialize_gravity','found NaN or +/-Inf in gravy_ypencil')
      if (notanumber(gravz_zpencil)) &
        call fatal_error('initialize_gravity','found NaN or +/-Inf in gravz_zpencil')
!
      call put_shared_variable('nu_epicycle',nu_epicycle,caller='initialize_gravity')
      call put_shared_variable('gravz_zpencil',gravz_zpencil)
      if (lreference_state) call put_shared_variable('gravx', gravx)
!
      if (n_adjust_sphersym>0 .and. lspherical_coords) then

        if (gravx_profile=='kepler' .and. mass_cent_body/=0.) then
          if (gravitational_const/=0.) then
            if (lroot) call warning('initialize_gravity', 'central body mass is ignored')
          else
            if (ipx==0) gravitational_const = -gravx_xpencil(l1)*x(l1)**2/mass_cent_body
            if (nprocx>1) call mpibcast_real(gravitational_const,comm=MPI_COMM_WORLD)
          endif
        else
          if (gravitational_const<=0.) &
            call fatal_error('initialize_gravity', 'positive gravitational constant needed')
        endif

        G4pi=gravitational_const*4.*pi
        ladjust_sphersym=.true.
        gravx_xpencil_0 = gravx_xpencil

      endif
  
      call keep_compiler_quiet(f)
!
    endsubroutine initialize_gravity
!***********************************************************************
    subroutine gravity_sphersym(f)
!
!  Adjusts spherically symmetric part of gravitational acceleration
!  according to density.
!
!  9-jun-15/MR: coded
!
      use Sub, only: meanyz
      use Mpicomm, only: mpirecv_real, mpisend_real

      real, dimension(mx,my,mz,mfarray) :: f

      integer :: l, la
      real :: integ
      real, dimension(nx) :: rhomean

      call meanyz(f,ilnrho,rhomean,lexp=.not.ldensity_nolog)

      if (ipx==0) then
        integ=0.; la=l1+1          ! -> rectangle rule applied in integration below
      else
        call mpirecv_real(integ,xlneigh,iproc)
        la=l1
      endif

      do l=la,l2
        integ = integ+x(l)**2*rhomean(l-l1+1)*(x(l)-x(l-1))
        gravx_xpencil(l) = gravx_xpencil_0(l) - G4pi*integ/x(l)**2
      enddo

      if (ipx<nprocx-1) call mpisend_real(integ,xuneigh,xuneigh)

!if (ipy==0 .and. ipz==0) print'(a,38(f6.3,",",1x))', 'gravx_xpencil=', gravx_xpencil
    endsubroutine gravity_sphersym
!***********************************************************************
    subroutine set_consistent_gravity(ginput,gtype,gprofile,lsuccess)
!
!  This subroutine checks, if the gravity paramters as type, profile and values
!  are set consistently with initial condition for example.
!
!  ginput     =     value for the gravity, GM    : 4, 10, 200
!  gtype      =     type of gravity              : 'gravx','gravy','gravz'
!  gprofile   =     profile of the gravity       : 'kepler','const'
!  lsuccess   =     switch, if it was successful : .true., .false.
!
!  13-jun-12/dhruba+joern: coded
!
      real :: ginput
      character (len=labellen) :: gtype,gprofile
      character (len=labellen) :: gprof
      logical :: lsuccess
      logical :: lconsistent=.true.
!
! check for consistency
!
      gprof=trim(gprofile)
      select case(trim(gtype))
        case('gravx')
          if (gprof/=gravx_profile) then
            lconsistent=.false.
            gravx_profile=gprof
          endif
          if (gravx/=ginput) then
            lconsistent=.false.
            gravx=ginput
          endif
        case('gravy')
          if (gprof/=gravy_profile) then
            lconsistent=.false.
            gravy_profile=gprof
          endif
          if (gravy/=ginput) then
            lconsistent=.false.
            gravy=ginput
          endif
        case('gravz')
          if (gprof/=gravz_profile) then
            lconsistent=.false.
            gravz_profile=gprof
          endif
          if (gravz/=ginput) then
            lconsistent=.false.
            gravz=ginput
          endif
      case default
        call fatal_error('set_consistent_gravity','gtype does not match any, aborting')
      endselect
      lsuccess=.true.
!
! gravity parameters set consistently.
!
    endsubroutine set_consistent_gravity
!***********************************************************************
    subroutine init_gg(f)
!
!  Initialise gravity; called from start.f90.
!
!  12-nov-04/anders: coded
!
      real, dimension (mx,my,mz,mfarray) :: f
!
!  Don't do anything
!
      call keep_compiler_quiet(f)
!
    endsubroutine init_gg
!***********************************************************************
    subroutine pencil_criteria_gravity()
!
!  All pencils that the Gravity module depends on are specified here.
!
!  20-nov-04/anders: coded
!  20-jan-15/MR: pencil request for rho1 added when reference_state is used
!
      lpenc_requested(i_gg)=.true.
      if (lanelastic) lpenc_requested(i_rho_anel)=.true.
      if (lreference_state) lpenc_requested(i_rho1)=.true.
!
      if (idiag_epot/=0 .or. idiag_epot/=0 .or. idiag_epotmx/=0 .or. &
          idiag_epotmy/=0 .or. idiag_epotmz/=0 .or. &
          idiag_epottot/=0) lpenc_diagnos(i_epot)=.true.
      if (idiag_epotuxmx/=0 .or. idiag_epotuzmz/=0) then
        lpenc_diagnos(i_epot)=.true.
        lpenc_diagnos(i_uu)=.true.
      endif
      if (idiag_ugm/=0) then
        lpenc_diagnos(i_uu)=.true.
        lpenc_diagnos(i_gg)=.true.
      endif
      if (idiag_epotmxy/=0) then
        lpenc_diagnos2d(i_epot)=.true.
      endif
      if (idiag_epotuxmxy/=0) then
        lpenc_diagnos2d(i_epot)=.true.
        lpenc_diagnos2d(i_uu)=.true.
      endif
!
    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 (lpencil_in(i_epot)) lpencil_in(i_rho)=.true.
!
    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
!
      real, dimension (mx,my,mz,mfarray) :: f
      type (pencil_case) :: p
!
      intent(in) :: f
      intent(inout) :: p
!
      if (lpencil(i_gg)) then
        p%gg(:,1) = gravx_xpencil(l1:l2)
        p%gg(:,2) = gravy_ypencil(m)
        p%gg(:,3) = gravz_zpencil(n)
      endif
!
      if (lpencil(i_epot)) p%epot=p%rho* &
          (potx_xpencil(l1:l2)+poty_ypencil(m)+potz_zpencil(n))
!
      call keep_compiler_quiet(f)
!
    endsubroutine calc_pencils_gravity
!***********************************************************************
    subroutine duu_dt_grav(f,df,p)
!
!  Add gravitational acceleration to gas and dust.
!
!  The special option lboussinesq_grav=T is applicable when |z|/H  << 1.
!  However, in the present formulation the resulting equations,
!  du/dt = -lnrho, and dlnrho/dt=-du/dz, lead to an instability
!  with the growth rate lambda = (1+i)*sqrt(k/2).
!
!  12-nov-04/anders: coded
!   5-dec-06/petri: added Boussinesq approximation
!  20-jan-15/MR: changes for use of reference state
!
      use Diagnostics
      use SharedVariables, only: get_shared_variable
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx,my,mz,mvar) :: df
      type (pencil_case) :: p
!
      intent(in) :: f,p
      intent(inout) :: df
!
      integer :: k
      real, dimension(nx,3) :: gg
      real, dimension(nx) :: refac
      real, dimension(:,:), pointer :: reference_state
!
!  Add gravity acceleration on gas.
!
      if (lhydro) then
        if (lboussinesq_grav) then
          if (lentropy) then
            if (lgravx_gas) df(l1:l2,m,n,iux)=df(l1:l2,m,n,iux)+p%ss*p%gg(:,1)
            if (lgravy_gas) df(l1:l2,m,n,iuy)=df(l1:l2,m,n,iuy)+p%ss*p%gg(:,2)
            if (lgravz_gas) df(l1:l2,m,n,iuz)=df(l1:l2,m,n,iuz)+p%ss*p%gg(:,3)
          else
            if (headtt) print*,'duu_dt_grav: lboussinesq_grav w/o lentropy not ok!'
          endif
        else if (lanelastic) then
! Now works for the linear anelastic formulation only
                if (lgravx_gas) df(l1:l2,m,n,iux)=df(l1:l2,m,n,iux)+p%gg(:,1)*&
                                p%rho_anel
                if (lgravy_gas) df(l1:l2,m,n,iuy)=df(l1:l2,m,n,iuy)+p%gg(:,2)*&
                                p%rho_anel
                if (lgravz_gas) df(l1:l2,m,n,iuz)=df(l1:l2,m,n,iuz)+p%gg(:,3)*&
                                 (-p%ss)
!                                p%rho_anel
        else
          if (lgravx_gas) gg(:,1)=p%gg(:,1)
          if (lgravy_gas) gg(:,2)=p%gg(:,2)
          if (lgravz_gas) gg(:,3)=p%gg(:,3)
!
! When reference state is used, gravity needs a correction factor rho'/rho=1-rho_0/rho.
! Note that the pencil case contains always the total quantities.
!
          if (lreference_state) then
            call get_shared_variable('reference_state',reference_state,caller='duu_dt_grav')  ! shouldn't be within the mn loop
            refac=1.-reference_state(:,iref_rho)*p%rho1
            if (lgravx_gas) gg(:,1)=gg(:,1)*refac
            if (lgravy_gas) gg(:,2)=gg(:,2)*refac
            if (lgravz_gas) gg(:,3)=gg(:,3)*refac
          endif
!
          if (lxyzdependence) then
            if (lgravx_gas) df(l1:l2,m,n,iux)=df(l1:l2,m,n,iux)+gg(:,1)*zdep(n)
            if (lgravy_gas) df(l1:l2,m,n,iuy)=df(l1:l2,m,n,iuy)+gg(:,2)
            if (lgravz_gas) df(l1:l2,m,n,iuz)=df(l1:l2,m,n,iuz)+gg(:,3)*xdep(l1:l2)
          else
            if (lgravx_gas) df(l1:l2,m,n,iux)=df(l1:l2,m,n,iux)+gg(:,1)
            if (lgravy_gas) df(l1:l2,m,n,iuy)=df(l1:l2,m,n,iuy)+gg(:,2)
            if (lgravz_gas) df(l1:l2,m,n,iuz)=df(l1:l2,m,n,iuz)+gg(:,3)
          endif
        endif
      endif
!
!  Add gravity acceleration on dust.
!
      if (ldustvelocity) then
        do k=1,ndustspec
          if (lgravx_dust) &
              df(l1:l2,m,n,iudx(k)) = df(l1:l2,m,n,iudx(k)) + p%gg(:,1)
          if (lgravy_dust) &
              df(l1:l2,m,n,iudy(k)) = df(l1:l2,m,n,iudy(k)) + p%gg(:,2)
          if (lgravz_dust) &
              df(l1:l2,m,n,iudz(k)) = df(l1:l2,m,n,iudz(k)) + p%gg(:,3)
        enddo
      endif
!
!  Gravity diagnostics.
!
      if (ldiagnos) then
        if (idiag_epot/=0) call sum_mn_name(p%epot,idiag_epot)
        if (idiag_ugm/=0) &
          call sum_mn_name(p%uu(:,1)*p%gg(:,1) + &
          p%uu(:,2)*p%gg(:,2) + &
          p%uu(:,3)*p%gg(:,3),idiag_ugm)
        if (idiag_epottot/=0) call integrate_mn_name(p%epot,idiag_epottot)
      endif
!
!  Gravity 1-D diagnostics.
!
      if (l1davgfirst) then
        call yzsum_mn_name_x(p%epot,idiag_epotmx)
        call yzsum_mn_name_x(p%epot*p%uu(:,1),idiag_epotuxmx)
        call xzsum_mn_name_y(p%epot,idiag_epotmy)
        call xysum_mn_name_z(p%epot,idiag_epotmz)
        call xysum_mn_name_z(p%epot*p%uu(:,3),idiag_epotuzmz)
      endif
!
!  Gravity 2-D diagnostics.
!
      if (l2davgfirst) then
        call zsum_mn_name_xy(p%epot,idiag_epotmxy)
        call zsum_mn_name_xy(p%epot*p%uu(:,1),idiag_epotuxmxy)
      endif
!
      call keep_compiler_quiet(f)
!
    endsubroutine duu_dt_grav
!***********************************************************************
    subroutine gravity_after_boundary(f)
!
!  For actions outside mn-loop. 
!  At the moment only adjustment of spherically symmetric gravity.
!
!  9-jun-15/MR: coded
!
      real, dimension(mx,my,mz,mfarray) :: f

      if (lfirst.and.ladjust_sphersym) then
        if ( mod(it, n_adjust_sphersym) == 0 ) call gravity_sphersym(f)
      endif

    endsubroutine gravity_after_boundary
!***********************************************************************
    subroutine potential_global(pot,pot0)
!
!  Calculates gravity potential globally
!
!  13-nov-04/anders: coded
!
      real, dimension (:,:,:) :: pot
      real, optional :: pot0
!
      call fatal_error('potential_global','this subroutine has been '// &
          'deprecated for gravity_simple')
!
      call keep_compiler_quiet(pot)
      if (present(pot0)) call keep_compiler_quiet(pot0)
!
    endsubroutine potential_global
!***********************************************************************
    subroutine potential_penc(xmn,ymn,zmn,pot,pot0,grav,rmn)
!
!  Calculates gravity potential on a pencil.
!
!  13-nov-04/anders: coded
!
      real, dimension (:) :: pot
      real, optional :: ymn,zmn,pot0
      real, optional, dimension (size(pot)) :: xmn,rmn
      real, optional, dimension (size(pot),3) :: grav
!
      intent(in) :: xmn,ymn,zmn,rmn
      intent(out) :: pot
!
!  Calculate potential from master pencils defined in initialize_gravity.
!
      if (lxyzdependence) then
        pot=potx_xpencil(l1:l2)*zdep(n) &
           +poty_ypencil(m) &
           +potz_zpencil(n)*xdep(l1:l2)
      else
        pot=potx_xpencil(l1:l2)+poty_ypencil(m)+potz_zpencil(n)
      endif
!
      if (present(xmn)) call keep_compiler_quiet(xmn)
      if (present(ymn)) call keep_compiler_quiet(ymn)
      if (present(zmn)) call keep_compiler_quiet(zmn)
      if (present(pot0)) call keep_compiler_quiet(pot0)
      if (present(grav)) call keep_compiler_quiet(grav)
      if (present(rmn)) call keep_compiler_quiet(rmn)
!
    endsubroutine potential_penc
!***********************************************************************
    subroutine potential_point(x,y,z,r, pot,pot0, grav)
!
!  Calculates gravity potential in one point.
!
!  13-nov-04/anders: coded
!  24-oct-06/bing: added constant gravity profiles
!
      real :: pot
      real, optional :: x,y,z,r
      real, optional :: pot0,grav
      real :: potx_xpoint,poty_ypoint,potz_zpoint
      real :: prof,xdep,zdep
!
      potx_xpoint=0.0
      poty_ypoint=0.0
      potz_zpoint=0.0
!
!  Selection different potentials at one point.
!  These entries should match those for potx_xpencil.
!
      if (present(x)) then
        select case (gravx_profile)
        case ('zero')
          if (lroot) print*,'potential_point: no x-gravity'
        case ('const')
          potx_xpoint=-gravx*(x-xinfty)
        case ('linear_zdep')
          zdep=(1.+kappa_z1*z+.5*(kappa_z1*z)**2)
          potx_xpoint=0.5*(x**2-xinfty**2)*nux_epicycle**2*zdep
        case ('kepler')
          potx_xpoint=-gravx/x
        case ('CZbot1')
          potx_xpoint=-gravx*(1./x-1./xyz0(1))
        case ('CZmid1')
          potx_xpoint=-gravx*(1./x-2./(xyz0(1)+xyz1(1)))
        case default
          call fatal_error('potential_point', &
               'gravx_profile='//gravx_profile//' not implemented')
        endselect
      endif
!
      if (present(y)) then
        select case (gravy_profile)
        case ('zero')
          if (lroot) print*,'potential_point: no y-gravity'
        case ('const')
          poty_ypoint=-gravy*(y-yinfty)
        case default
          call fatal_error('potential_point', &
               'gravy_profile='//gravy_profile//' not implemented')
        endselect
      endif
!
      if (present(z)) then
        select case (gravz_profile)
        case ('zero')
          if (lroot) print*,'potential_point: no z-gravity'
        case ('const')
          if (zclip==impossible) then
            potz_zpoint=-gravz*(z-zinfty)
          else
            potz_zpoint=-gravz*max(z-zinfty,zclip-zinfty)
          endif
        case ('linear')
          potz_zpoint=0.5*(z**2-zinfty**2)*nu_epicycle**2
        case ('spherical')
          potz_zpoint=0.5*nu_epicycle**2*z1**2*alog(1.0+(z/z1)**2)
        case ('linear_xdep')
          xdep=(1.+kappa_x1*x+.5*(kappa_x1*x)**2)
          potz_zpoint=0.5*(z**2-zinfty**2)*nu_epicycle**2*xdep
        case ('linear_smoothed')
          prof = 1. + (z/zref)**(2*n_pot)
          potz_zpoint = 0.5*(nu_epicycle*z)**2/prof**(1./n_pot)
        case ('tanh')
          potz_zpoint=gravz*zref*alog(cosh(z/zref))
        case default
          call fatal_error('potential_point', &
               'gravz_profile='//gravz_profile//' not implemented')
        endselect
      endif
!
      pot = potx_xpoint + poty_ypoint + potz_zpoint
!
      if (present(r)) call keep_compiler_quiet(r)
      if (present(grav)) call keep_compiler_quiet(grav)
      if (present(pot0)) call keep_compiler_quiet(pot0)
!
    endsubroutine potential_point
!***********************************************************************
    subroutine acceleration_penc(gg)
!
!  Calculates gravitational acceleration on a pencil.
!
!  21-apr-07/tobi: adapted from potential_penc
!
      real, dimension (:,:), intent (out) :: gg
!
!  Calculate acceleration from master pencils defined in initialize_gravity.
!
      if (size(gg,2)/=3) then
        call fatal_error('acceleration_penc','Expecting a 3-vector pencil')
      endif
!
!  Note: the following would not yet work if lxyzdependence is set to true.
!
      select case (size(gg,1))
      case (nx)
        gg(:,1) = gravx_xpencil(l1:l2)
        !ABlater: gg(:,1) = gravx_xpencil(l1:l2)*zdep(n)
      case (mx)
        gg(:,1) = gravx_xpencil
        !ABlater: gg(:,1) = gravx_xpencil*zdep
      case default
        call fatal_error('acceleration_penc','Wrong pencil size.')
      endselect
!
      gg(:,2) = gravy_ypencil(m)
      gg(:,3) = gravz_zpencil(n)
!
      !ABlater: gg(:,3) = gravz_zpencil(n)*xdep(:)
!
    endsubroutine acceleration_penc
!***********************************************************************
    subroutine acceleration_penc_1D(gr)
!
!  Calculates gravitational acceleration on a pencil.
!
!  21-apr-07/tobi: adapted from potential_penc
!
      real, dimension (nx), intent (out) :: gr
!
!  Calculate acceleration from master pencils defined in initialize_gravity.
!
      call fatal_error('acceleration_penc_1D','Not implemented')
!
      call keep_compiler_quiet(gr)
!
    endsubroutine acceleration_penc_1D
!***********************************************************************
    subroutine acceleration_point(x,y,z,r,g_r)
!
!  Gravity in one point
!
!  18-nov-08/wlad: coded.
!
      real :: g_r
      real, optional :: x,y,z,r
!
      intent(in)  :: x,y,z,r
      intent(out) :: g_r
!
      call fatal_error('gravity_simple','acceleration_point not implemented')
!
      g_r=0.0
!
      call keep_compiler_quiet(present(x))
      call keep_compiler_quiet(present(y))
      call keep_compiler_quiet(present(z))
      call keep_compiler_quiet(present(r))
!
    endsubroutine acceleration_point
!***********************************************************************
    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 rprint_gravity(lreset,lwrite)
!
!  Reads and registers print parameters relevant for gravity advance.
!
!  12-jun-04/axel: adapted from grav_z
!
      use Diagnostics, only: parse_name
!
      logical :: lreset
      logical, optional :: lwrite
!
      integer :: iname, inamex, inamey, inamez, inamexy
      logical :: lwr
!
      lwr = .false.
      if (present(lwrite)) lwr=lwrite
!
!  Reset everything in case of reset.
!  (this needs to be consistent with what is defined above!)
!
      if (lreset) then
        idiag_epot=0; idiag_epottot=0; idiag_ugm=0
        idiag_epotmx=0; idiag_epotuxmx=0; idiag_epotmy=0; idiag_epotmz=0;
        idiag_epotmxy=0; idiag_epotuzmz=0; idiag_epotuxmxy=0
      endif
!
!  iname runs through all possible names that may be listed in print.in.
!
      do iname=1,nname
        call parse_name(iname,cname(iname),cform(iname),'epot',idiag_epot)
        call parse_name(iname,cname(iname),cform(iname),'epottot',idiag_epottot)
        call parse_name(iname,cname(iname),cform(iname),'ugm',idiag_ugm)
      enddo
!
!  Check for those quantities for which we want yz-averages.
!
      do inamex=1,nnamex
        call parse_name(inamex,cnamex(inamex),cformx(inamex),'epotmx', &
            idiag_epotmx)
        call parse_name(inamex,cnamex(inamex),cformx(inamex),'epotuxmx', &
            idiag_epotuxmx)
      enddo
!
!  Check for those quantities for which we want xz-averages.
!
      do inamey=1,nnamey
        call parse_name(inamey,cnamey(inamey),cformy(inamey),'epotmy', &
            idiag_epotmy)
      enddo
!
!  Check for those quantities for which we want xy-averages.
!
      do inamez=1,nnamez
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'epotmz', &
            idiag_epotmz)
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'epotuzmz', &
            idiag_epotuzmz)
      enddo
!
!  Check for those quantities for which we want z-averages.
!
      do inamexy=1,nnamexy
        call parse_name(inamexy,cnamexy(inamexy),cformxy(inamexy), &
            'epotmxy', idiag_epotmxy)
        call parse_name(inamexy,cnamexy(inamexy),cformxy(inamexy), &
            'epotuxmxy', idiag_epotuxmxy)
      enddo
!
!  Write column, idiag_XYZ, where our variable XYZ is stored.
!  IDL needs this even if everything is zero.
!
      if (lwr) then
        write(3,*) 'igg=',igg
      endif
!
    endsubroutine rprint_gravity
!***********************************************************************
    subroutine compute_gravity_star(f,wheat,luminosity,star_cte)
!
!  5-jan-10/boris: coded
!
      real, dimension (mx,my,mz,mfarray) :: f
      real :: wheat, luminosity, star_cte
!
      call keep_compiler_quiet(f)
      call keep_compiler_quiet(wheat,luminosity,star_cte)
!
    endsubroutine compute_gravity_star
!***********************************************************************
    subroutine get_xgravity(xgrav)
!
!  Used from the initial conditions
!
!  04-oct-10/bing: coded
!
      real, dimension(mx) :: xgrav
!
      xgrav = gravx_xpencil
!
    endsubroutine get_xgravity
!***********************************************************************
    logical function is_constant_zgrav()
!
!  15-apr-15/MR: coded
!
      is_constant_zgrav = gravx_profile=='zero'.and.gravy_profile=='zero'.and.gravz_profile=='const'

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