! $Id$
!
!** 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 :: ldensity = .true.
! CPARAM logical, parameter :: lanelastic = .false.
! CPARAM logical, parameter :: lboussinesq = .false.
!
! MVAR CONTRIBUTION 1
! MAUX CONTRIBUTION 0
!
! PENCILS PROVIDED rhos; rhos1
! PENCILS PROVIDED grhos(3); glnrhos(3); ugrhos
!
! PENCILS PROVIDED rho; rho1
! PENCILS PROVIDED grho(3); glnrho(3)
! PENCILS PROVIDED ekin
!
! PENCILS PROVIDED lnrho
! PENCILS PROVIDED ugrho; uglnrho
! PENCILS PROVIDED transprho
! PENCILS PROVIDED del2rho; del2lnrho; del6lnrho
! PENCILS PROVIDED hlnrho(3,3)
! PENCILS PROVIDED uij5glnrho(3); sglnrho(3)
! PENCILS PROVIDED totenergy_rel
!
!***************************************************************
module Density
!
  use Cparam
  use Cdata
  use General, only: keep_compiler_quiet
  use Messages
!
  use DensityMethods
!
  implicit none
!
  include 'density.h'
  integer :: pushpars2c, pushdiags2c  ! should be procedure pointer (F2003)
!
  character(len=labellen), dimension(ninit) :: initrho = 'nothing'
  real, dimension(ninit) :: amplrho = 0.0
  logical :: lconserve_mass = .false.
  logical :: ldiff_shock = .false.
  logical :: ldiff_hyper3_mesh = .false.
  logical :: lmassdiff_fix = .false.
  real :: density_floor = 0.0
  real :: diffrho_shock = 0.0
  real :: diffrho_hyper3_mesh = 0.0
  real, dimension(3) :: beta_glnrho_global=0.0, beta_glnrho_scaled=0.0
  logical :: lrelativistic_eos=.false.
!
  namelist /density_init_pars/ initrho, amplrho, beta_glnrho_global, lconserve_mass, lmassdiff_fix
!
  namelist /density_run_pars/ density_floor, diffrho_hyper3_mesh, diffrho_shock, lconserve_mass, lmassdiff_fix
!
!  Diagnostic Variables
!
  integer :: idiag_mass = 0     ! DIAG_DOC: $\int\rho\,d^3x$
  integer :: idiag_rhomin = 0   ! DIAG_DOC: $\min\left|\rho\right|$
  integer :: idiag_rhomax = 0   ! DIAG_DOC: $\max\left|\rho\right|$
  integer :: idiag_drhom = 0    ! DIAG_DOC: $\langle\Delta\rho/\rho_0\rangle$
  integer :: idiag_drho2m = 0   ! DIAG_DOC: $\langle\left(\Delta\rho/\rho_0\right)^2\rangle$
  integer :: idiag_drhorms = 0  ! DIAG_DOC: $\langle\Delta\rho/\rho_0\rangle_{rms}$
  integer :: idiag_drhomax = 0  ! DIAG_DOC: $\max\left|\Delta\rho/\rho_0\right|$
!
!  xy-averages
!
  integer :: idiag_drhomz = 0    ! XYAVG_DOC: $\langle\Delta\rho/\rho_0\rangle_{xy}$
  integer :: idiag_drho2mz = 0   ! XYAVG_DOC: $\langle\left(\Delta\rho/\rho_0\right)^2\rangle_{xy}$
!
!  xz-averages
!
  integer :: idiag_drhomy = 0    ! XZAVG_DOC: $\langle\Delta\rho/\rho_0\rangle_{xz}$
  integer :: idiag_drho2my = 0   ! XZAVG_DOC: $\langle\left(\Delta\rho/\rho_0\right)^2\rangle_{xz}$
!
!  yz-averages
!
  integer :: idiag_drhomx = 0    ! YZAVG_DOC: $\langle\Delta\rho/\rho_0\rangle_{yz}$
  integer :: idiag_drho2mx = 0   ! YZAVG_DOC: $\langle\left(\Delta\rho/\rho_0\right)^2\rangle_{yz}$
!
!  y-averages
!
  integer :: idiag_drhomxz = 0   ! YAVG_DOC: $\langle\Delta\rho/\rho_0\rangle_y$
  integer :: idiag_drho2mxz = 0  ! YAVG_DOC: $\langle\left(\Delta\rho/\rho_0\right)^2\rangle_y$
!
!  z-averages
!
  integer :: idiag_drhomxy = 0   ! ZAVG_DOC: $\langle\Delta\rho/\rho_0\rangle_z$
  integer :: idiag_drho2mxy = 0  ! ZAVG_DOC: $\langle\left(\Delta\rho/\rho_0\right)^2\rangle_z$
  integer :: idiag_sigma = 0     ! ZAVG_DOC; $\Sigma\equiv\int\rho\,\mathrm{d}z$
!
!  Module Variables
!
  real, dimension(mz) :: rho0z = 0.0
  real, dimension(mz) :: dlnrho0dz = 0.0
  real :: mass0 = 0.0
!
!  Dummy Variables
!
  real, dimension(mz,1) :: lnrhomz
  logical :: lcalc_lnrhomean = .false.
  logical :: lupw_lnrho = .false.
!
  contains
!***********************************************************************
    subroutine register_density()
!
!  Register Density module.
!
!  28-feb-13/ccyang: coded.
!
      use FArrayManager, only: farray_register_pde
      use SharedVariables, only: put_shared_variable
!
!  Register relative density as dynamical variable rho.
!
      call farray_register_pde('rho', irho)
!
!  This module conflicts with Gravity module.
!
      if (lgrav) call fatal_error('register_density', 'conflicting with Gravity module. ')
!
!  This module does not consider logarithmic density.
!
      ldensity_nolog = .true.
!
!  Communicate lrelativistic_eos to entropy too.
!
      call put_shared_variable('lrelativistic_eos',lrelativistic_eos) 
!
!
!  Identify version number.
!
      if (lroot) call svn_id("$Id$")
!
    endsubroutine register_density
!***********************************************************************
    subroutine initialize_density(f)
!
!  Perform any post-parameter-read initialization i.e. calculate derived
!  parameters.
!
!  07-sep-14/ccyang: coded.
!
      use EquationOfState, only: select_eos_variable, get_stratz
      use SharedVariables, only: put_shared_variable
      use DensityMethods,  only: initialize_density_methods
!
      real, dimension(mx,my,mz,mfarray), intent(inout) :: f
!
      integer :: ierr
!
!  Tell the equation of state that we're here and what f variable we use.
!
      call select_eos_variable('rho', irho)
      call initialize_density_methods
!
!  Get density stratification.
!
      if (lstratz) call get_stratz(z, rho0z, dlnrho0dz)
!
!  Disable the force-free considerations.
!
      call put_shared_variable('lffree', .false., ierr)
      if (ierr /= 0) call fatal_error('initialize_density', 'cannot share variable lffree. ')
!
!  Check the switches.
!
      shock: if (diffrho_shock > 0.0) then
        ldiff_shock = .true.
        if (lroot) print *, 'initialize_density: shock mass diffusion; diffrho_shock = ', diffrho_shock
      endif shock
!
      hyper3_mesh: if (diffrho_hyper3_mesh > 0.0) then
        ldiff_hyper3_mesh = .true.
        if (lroot) print *, 'initialize_density: mesh hyper-diffusion; diffrho_hyper3_mesh = ', diffrho_hyper3_mesh
      endif hyper3_mesh
!
!  Get the total mass.
!
      if (lconserve_mass) mass0 = total_mass(f)
!
    endsubroutine initialize_density
!***********************************************************************
    subroutine init_lnrho(f)
!
!  Initializes the density field.
!
!  27-feb-13/ccyang: coded.
!
      use Initcond, only: gaunoise
!
      real, dimension(mx,my,mz,mfarray), intent(inout) :: f
!
      integer :: j
!
      init: do j = 1, ninit
        if (initrho(j) == 'nothing') cycle init
!
        init_cond: select case (initrho(j))
!
        case ('zero') init_cond
          f(:,:,:,irho) = 0.0
!
        case ('const') init_cond
          f(:,:,:,irho) = amplrho(j)
!
        case ('noise', 'gaussian') init_cond
          call gaunoise(amplrho(j), f, irho, irho)
!
        case default init_cond
          call fatal_error('init_lnrho', 'unknown initial condition ' // trim(initrho(j)))
!
        endselect init_cond
!
      enddo init
!
    endsubroutine init_lnrho
!***********************************************************************
    subroutine pencil_criteria_density()
!
!  All pencils that the Density module depends on are specified here.
!
!  02-nov-13/ccyang: coded.
!
      lpenc_requested(i_rhos) = .true.
      lpenc_requested(i_uu) = .true.
      lpenc_requested(i_divu) = .true.
      lpenc_requested(i_ugrhos) = .true.
!
      shock: if (ldiff_shock) then
        lpenc_requested(i_shock) = .true.
        lpenc_requested(i_gshock) = .true.
        lpenc_requested(i_grhos) = .true.
      endif shock
!
      massdiff: if (lmassdiff_fix) then
        hydro: if (lhydro) then
          lpenc_requested(i_rhos1) = .true.
          lpenc_requested(i_uu) = .true.
        endif hydro
        if (lthermal_energy) lpenc_requested(i_u2) = .true.
      endif massdiff
!
!  Diagnostic Pencils
!
      if (idiag_mass /= 0 .or. idiag_rhomin /= 0 .or. idiag_rhomax /= 0) lpenc_diagnos(i_rho) = .true.
      if (idiag_sigma /= 0) lpenc_diagnos(i_rho) = .true.
!
    endsubroutine pencil_criteria_density
!***********************************************************************
    subroutine pencil_interdep_density(lpencil_in)
!
!  Interdependency among pencils from the Density module is specified here.
!
!  02-nov-13/ccyang: coded.
!
      logical, dimension(npencils) :: lpencil_in
!
      if (lpencil_in(i_rhos1)) lpencil_in(i_rhos) = .true.
!
      glnrhos: if (lpencil_in(i_glnrhos)) then
        lpencil_in(i_rhos1) = .true.
        lpencil_in(i_grhos) = .true.
      endif glnrhos
!
      ugrhos: if (lpencil_in(i_ugrhos)) then
        lpencil_in(i_grhos) = .true.
        lpencil_in(i_uu) = .true.
      endif ugrhos
!
      if (lpencil_in(i_rho)) lpencil_in(i_rhos) = .true.
!
      if (lpencil_in(i_rho1)) lpencil_in(i_rho) = .true.
!
      grho: if (lpencil_in(i_grho)) then
        lpencil_in(i_rho) = .true.
        lpencil_in(i_grhos) = .true.
      endif grho
!
      glnrho: if (lpencil_in(i_glnrho)) then
        lpencil_in(i_rho1) = .true.
        lpencil_in(i_grho) = .true.
      endif glnrho
!
      ekin: if (lpencil_in(i_ekin)) then
        lpencil_in(i_rho)=.true.
        lpencil_in(i_u2)=.true.
      endif ekin
!
    endsubroutine pencil_interdep_density
!***********************************************************************
    subroutine calc_pencils_density(f, p)
!
!  Calculate Density pencils.
!  Most basic pencils should come first, as others may depend on them.
!
!  02-nov-13/ccyang: coded.
!
      use Sub, only: grad, u_dot_grad
      use EquationOfState, only: cs20
!
      real, dimension(mx,my,mz,mfarray), intent(in) :: f
      type(pencil_case), intent(inout) :: p
!
! rhos: density scaled by stratification
!
      if (lpencil(i_rhos)) p%rhos = 1.0 + f(l1:l2,m,n,irho)
!
! rhos1
!
      if (lpencil(i_rhos1)) p%rhos1 = 1.0 / p%rhos
!
! grhos
!
      if (lpencil(i_grhos)) call grad(f, irho, p%grhos)
!
! glnrhos
!
      if (lpencil(i_glnrhos)) p%glnrhos = spread(p%rhos1,2,3) * p%grhos
!
! ugrhos
!
      if (lpencil(i_ugrhos)) call u_dot_grad(f, irho, p%grhos, p%uu, p%ugrhos)
!
! rho
!
      if (lpencil(i_rho)) p%rho = rho0z(n) * p%rhos
!
! rho1
!
      if (lpencil(i_rho1)) p%rho1 = 1.0 / p%rho
!
! grho
!
      grho: if (lpencil(i_grho)) then
        p%grho = rho0z(n) * p%grhos
        p%grho(:,3) = p%grho(:,3) + dlnrho0dz(n) * p%rho
      endif grho
!
! glnrho
!
      if (lpencil(i_glnrho)) p%glnrho = spread(p%rho1,2,3) * p%grho
!
! ekin
!
      if (lpencil(i_ekin)) p%ekin = 0.5 * p%rho * p%u2
!
!  Currently not required pencils.
!
      lnrho: if (lpencil(i_lnrho)) then
        call fatal_error('calc_pencils_density', 'lnrho is not available. ')
        p%lnrho = 0.0
      endif lnrho
!
      ugrho: if (lpencil(i_ugrho)) then
        call fatal_error('calc_pencils_density', 'ugrho is not available. ')
        p%ugrho = 0.0
      endif ugrho
!
      uglnrho: if (lpencil(i_uglnrho)) then
        call fatal_error('calc_pencils_density', 'uglnrho is not available. ')
        p%uglnrho = 0.0
      endif uglnrho
!
      transprho: if (lpencil(i_transprho)) then
        call fatal_error('calc_pencils_density', 'transprho is not available. ')
        p%transprho = 0.0
      endif transprho
!
      del2rho: if (lpencil(i_del2rho)) then
        call fatal_error('calc_pencils_density', 'del2rho is not available. ')
        p%del2rho = 0.0
      endif del2rho
!
      del2lnrho: if (lpencil(i_del2lnrho)) then
        call fatal_error('calc_pencils_density', 'del2lnrho is not available. ')
        p%del2lnrho = 0.0
      endif del2lnrho
!
      del6lnrho: if (lpencil(i_del6lnrho)) then
        call fatal_error('calc_pencils_density', 'del6lnrho is not available. ')
        p%del6lnrho = 0.0
      endif del6lnrho
!
      hlnrho: if (lpencil(i_hlnrho)) then
        call fatal_error('calc_pencils_density', 'hlnrho is not available. ')
        p%hlnrho = 0.0
      endif hlnrho
!
      uij5glnrho: if (lpencil(i_uij5glnrho)) then
        call fatal_error('calc_pencils_density', 'uij5glnrho is not available. ')
        p%uij5glnrho = 0.0
      endif uij5glnrho
!
      sglnrho: if (lpencil(i_sglnrho)) then
        call fatal_error('calc_pencils_density', 'sglnrho is not available. ')
        p%sglnrho = 0.0
      endif sglnrho
!
    endsubroutine calc_pencils_density
!***********************************************************************
    subroutine dlnrho_dt(f,df,p)
!
!  Continuity equation.
!
!  02-nov-13/ccyang: coded.
!  24-aug-14/ccyang: add mass diffusion correction.
!
      use Sub, only: identify_bcs, dot_mn, del2
      use Deriv, only: der6
      use Special, only: special_calc_density
      use Diagnostics
!
      real, dimension(mx,my,mz,mfarray), intent(in) :: f
      real, dimension(mx,my,mz,mvar), intent(inout) :: df
      type(pencil_case), intent(in) :: p
!
      real, dimension(nx) :: fdiff, del2rhos, penc, src_density,diffus_diffrho,diffus_diffrho3
      integer :: j
!
!  Start the clock for this procedure.
!
      call timing('dlnrho_dt', 'entered', mnloop=.true.)
!
!  Identify module and boundary conditions.
!
      if (headtt .or. ldebug) print *, 'dlnrho_dt: SOLVE the continuity equation. '
      if (headtt) call identify_bcs('rho', irho)
!
!  Find the rate of change.
!
      penc = p%uu(:,3) * dlnrho0dz(n)
      df(l1:l2,m,n,irho) = df(l1:l2,m,n,irho) - p%ugrhos - p%rhos * (p%divu + penc)
      if (lfirst .and. ldt) then
        src_density = abs(penc)
        maxsrc = max(maxsrc, src_density)
      endif
!
      fdiff = 0.0
!
!  Shock mass diffusion
!
      shock: if (ldiff_shock) then
        call del2(f, irho, del2rhos)
        call dot_mn(p%gshock, p%grhos, penc)
        fdiff = fdiff + diffrho_shock * (p%shock * del2rhos + penc)
        if (lfirst .and. ldt) diffus_diffrho = diffus_diffrho + diffrho_shock * dxyz_2 * p%shock
      endif shock
!
!  Mesh hyper-diffusion
!
      hyper3: if (ldiff_hyper3_mesh) then
        dir: do j = 1, 3
          call der6(f, irho, penc, j, ignoredx=.true.)
          fdiff = fdiff + diffrho_hyper3_mesh * penc * dline_1(:,j)
        enddo dir
        if (lfirst .and. ldt) diffus_diffrho3 = diffus_diffrho3 &
                                              + diffrho_hyper3_mesh * sum(dline_1,2)
      endif hyper3
!
      df(l1:l2,m,n,irho) = df(l1:l2,m,n,irho) + fdiff
!
!  Mass diffusion corrections.
!
      massdiff: if (lmassdiff_fix) then
        hydro: if (lhydro) then
          penc = fdiff * p%rhos1
          forall(j = iux:iuz) df(l1:l2,m,n,j) = df(l1:l2,m,n,j) - penc * p%uu(:,j-iuu+1)
        endif hydro
        energy: if (lthermal_energy) then
          df(l1:l2,m,n,ieth) = df(l1:l2,m,n,ieth) - 0.5 * rho0z(n) * fdiff * p%u2
        elseif (lenergy) then energy
          call fatal_error('dlnrho_dt', 'mass diffusion correction for this energy equation is not implemented. ')
        endif energy
      endif massdiff
!
!  Call optional user-defined calculations.
!
      if (lspecial) call special_calc_density(f, df, p)
!
!  2D averages
!
      avg_2d: if (l2davgfirst) then
        penc = f(l1:l2,m,n,irho)
        if (idiag_drhomxz /= 0) call ysum_mn_name_xz(penc, idiag_drhomxz)
        if (idiag_drho2mxz /= 0) call ysum_mn_name_xz(penc**2, idiag_drho2mxz)
        if (idiag_drhomxy /= 0) call zsum_mn_name_xy(penc, idiag_drhomxy)
        if (idiag_drho2mxy /= 0) call zsum_mn_name_xy(penc**2, idiag_drho2mxy)
        if (idiag_sigma /= 0) call zsum_mn_name_xy(p%rho, idiag_sigma, lint=.true.)
      endif avg_2d
!
!  1D averages
!
      avg_1d: if (l1davgfirst) then
        penc = f(l1:l2,m,n,irho)
!       xy-averages
        if (idiag_drhomz /= 0) call xysum_mn_name_z(penc, idiag_drhomz)
        if (idiag_drho2mz /= 0) call xysum_mn_name_z(penc**2, idiag_drho2mz)
!       xz-averages
        if (idiag_drhomy /= 0) call xzsum_mn_name_y(penc, idiag_drhomy)
        if (idiag_drho2my /= 0) call xzsum_mn_name_y(penc**2, idiag_drho2my)
!       yz-averages
        if (idiag_drhomx /= 0) call yzsum_mn_name_x(penc, idiag_drhomx)
        if (idiag_drho2mx /= 0) call yzsum_mn_name_x(penc**2, idiag_drho2mx)
      endif avg_1d
!
!  Diagnostics
!
      diagnos: if (ldiagnos) then
!
        rho: if (idiag_mass /= 0 .or. idiag_rhomin /= 0 .or. idiag_rhomax /= 0) then
          penc = p%rho
          if (idiag_mass /= 0) call integrate_mn_name(penc, idiag_mass)
          if (idiag_rhomin /= 0) call max_mn_name(-penc, idiag_rhomin, lneg=.true.)
          if (idiag_rhomax /= 0) call max_mn_name(penc, idiag_rhomax)
        endif rho
!
        drho: if (idiag_drhom /= 0 .or. idiag_drho2m /= 0 .or. idiag_drhorms /= 0 .or. idiag_drhomax /= 0) then
          penc = f(l1:l2,m,n,irho)
          if (idiag_drhom /= 0) call sum_mn_name(penc, idiag_drhom)
          if (idiag_drho2m /= 0) call sum_mn_name(penc**2, idiag_drho2m)
          if (idiag_drhorms /= 0) call sum_mn_name(penc**2, idiag_drhorms, lsqrt=.true.)
          if (idiag_drhomax /= 0) call max_mn_name(abs(penc), idiag_drhomax)
        endif drho
!
      endif diagnos
!
!  Stop the clock.
!
      call timing('dlnrho_dt', 'finished', mnloop=.true.)
!
    endsubroutine dlnrho_dt
!***********************************************************************
    subroutine calc_diagnostics_density(f,p)

      real, dimension (mx,my,mz,mfarray) :: f
      type(pencil_case) :: p

      call keep_compiler_quiet(f)
      call keep_compiler_quiet(p)

    endsubroutine calc_diagnostics_density
!***********************************************************************
    subroutine split_update_density(f)
!
!  Integrate operator split terms and/or perform post-time-step
!  processing.
!
!  18-jun-14/ccyang: coded.
!
      real, dimension(mx,my,mz,mfarray), intent(inout) :: f
!
      real :: a, a1
!
!  Conserve mass and momentum.
!
      consrv: if (lconserve_mass) then
        a = mass0 / total_mass(f)
        a1 = 1.0 / a
        zscan: do n = n1, n2
          yscan: do m = m1, m2
            f(l1:l2,m,n,irho) = a * f(l1:l2,m,n,irho) + (a - 1.0)
            if (lhydro) f(l1:l2,m,n,iux:iuz) = a1 * f(l1:l2,m,n,iux:iuz)
          enddo yscan
        enddo zscan
      endif consrv
!
    endsubroutine split_update_density
!***********************************************************************
    subroutine impose_density_floor(f)
!
!  Check negative density or impose a density floor.
!
!  27-feb-13/ccyang: coded.
!
      real, dimension(mx,my,mz,mfarray), intent(inout) :: f
!
      integer :: i, j, k
      real :: drhos_min
!
!  Impose density floor or
!
      dens_floor: if (density_floor > 0.0) then
        scan1: do k = n1, n2
          drhos_min = (density_floor - rho0z(k)) / rho0z(k)
          where(f(l1:l2,m1:m2,k,irho) < drhos_min) f(l1:l2,m1:m2,k,irho) = drhos_min
        enddo scan1
!
!  Trap any negative density.
!
      else dens_floor
        neg_dens: if (any(f(l1:l2,m1:m2,n1:n2,irho) <= -1.0)) then
          scan2z: do k = n1, n2
            scan2y: do j = m1, m2
              scan2x: do i = l1, l2
                if (f(i,j,k,irho) <= -1.0) print 10, rho0z(k) * (1.0 + f(i,j,k,irho)), x(i), y(j), z(k)
                10 format (1x, 'Negative density ', es13.6, ' is detected at x = ', es10.3, ', y = ', es10.3, ', and z = ', es10.3)
              enddo scan2x
            enddo scan2y
          enddo scan2z
          call fatal_error_local('impose_density_floor', 'detected negative density. ')
        endif neg_dens
!
      endif dens_floor
!
    endsubroutine impose_density_floor
!***********************************************************************
    subroutine read_density_init_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=density_init_pars, IOSTAT=iostat)
!
    endsubroutine read_density_init_pars
!***********************************************************************
    subroutine write_density_init_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=density_init_pars)
!
    endsubroutine write_density_init_pars
!***********************************************************************
    subroutine read_density_run_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=density_run_pars, IOSTAT=iostat)
!
    endsubroutine read_density_run_pars
!***********************************************************************
    subroutine write_density_run_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=density_run_pars)
!
    endsubroutine write_density_run_pars
!***********************************************************************
    subroutine rprint_density(lreset, lwrite)
!
!  Reads and registers print parameters relevant for continuity equation.
!
!  27-feb-13/ccyang: coded.
!
      use Diagnostics, only: parse_name
      use FArrayManager, only: farray_index_append
!
      logical, intent(in) :: lreset
      logical, intent(in), optional :: lwrite
!
      logical :: lwr
      integer :: iname
!
      lwr = .false.
      if (present(lwrite)) lwr = lwrite
!
!  Reset everything in case of reset.
!
      reset: if (lreset) then
!       Diagnostics
        idiag_mass = 0
        idiag_rhomin = 0
        idiag_rhomax = 0
        idiag_drhom = 0
        idiag_drho2m = 0
        idiag_drhorms = 0
        idiag_drhomax = 0
!       xy-averages
        idiag_drhomz = 0
        idiag_drho2mz = 0
!       xz-averages
        idiag_drhomy = 0
        idiag_drho2my = 0
!       yz-averages
        idiag_drhomx = 0
        idiag_drho2mx = 0
!       y-averages
        idiag_drhomxz = 0
        idiag_drho2mxz = 0
!       z-averages
        idiag_sigma = 0
        idiag_drhomxy = 0
        idiag_drho2mxy = 0
      endif reset
!
!  Check for diagnostic variables listed in print.in.
!
      whole: do iname = 1, nname
        call parse_name(iname, cname(iname), cform(iname), 'mass', idiag_mass)
        call parse_name(iname, cname(iname), cform(iname), 'rhomin', idiag_rhomin)
        call parse_name(iname, cname(iname), cform(iname), 'rhomax', idiag_rhomax)
        call parse_name(iname, cname(iname), cform(iname), 'drhom', idiag_drhom)
        call parse_name(iname, cname(iname), cform(iname), 'drho2m', idiag_drho2m)
        call parse_name(iname, cname(iname), cform(iname), 'drhorms', idiag_drhorms)
        call parse_name(iname, cname(iname), cform(iname), 'drhomax', idiag_drhomax)
      enddo whole
!
!  Check for xy-averages listed in xyaver.in.
!
      xyaver: do iname = 1, nnamez
        call parse_name(iname, cnamez(iname), cformz(iname), 'drhomz', idiag_drhomz)
        call parse_name(iname, cnamez(iname), cformz(iname), 'drho2mz', idiag_drho2mz)
      enddo xyaver
!
!  Check for xz-averages listed in xzaver.in.
!
      xzaver: do iname = 1, nnamey
        call parse_name(iname, cnamey(iname), cformy(iname), 'drhomy', idiag_drhomy)
        call parse_name(iname, cnamey(iname), cformy(iname), 'drho2my', idiag_drho2my)
      enddo xzaver
!
!  Check for yz-averages listed in yzaver.in.
!
      yzaver: do iname = 1, nnamex
        call parse_name(iname, cnamex(iname), cformx(iname), 'drhomx', idiag_drhomx)
        call parse_name(iname, cnamex(iname), cformx(iname), 'drho2mx', idiag_drho2mx)
      enddo yzaver
!
!  Check for y-averages listed in yaver.in.
!
      yaver: do iname = 1, nnamexz
        call parse_name(iname, cnamexz(iname), cformxz(iname), 'drhomxz', idiag_drhomxz)
        call parse_name(iname, cnamexz(iname), cformxz(iname), 'drho2mxz', idiag_drho2mxz)
      enddo yaver
!
!  Check for z-averages listed in zaver.in.
!
      zaver: do iname = 1, nnamexy
        call parse_name(iname, cnamexy(iname), cformxy(iname), 'drhomxy', idiag_drhomxy)
        call parse_name(iname, cnamexy(iname), cformxy(iname), 'drho2mxy', idiag_drho2mxy)
        call parse_name(iname, cnamexy(iname), cformxy(iname), 'sigma', idiag_sigma)
      enddo zaver
!
!  check for those quantities for which we want video slices
!
      if (lwrite_slices) then 
        where(cnamev=='rho'.or.cnamev=='drho') cformv='DEFINED'
      endif
!
!  Write column where which density variable is stored.
!
      if (lwr) then
        call farray_index_append('ilnrho',0)
        call farray_index_append('irho',irho)
      endif
!
    endsubroutine rprint_density
!***********************************************************************
    subroutine get_slices_density(f,slices)
!
!  Write slices for animation of Density variables.
!
!  27-feb-13/ccyang: coded.
!
      use Slices_methods, only: assign_slices_scal

      real, dimension(mx,my,mz,mfarray), intent(in) :: f
      type(slice_data), intent(inout) :: slices
!
      var: select case (trim(slices%name))
!
      case ('drho');  call assign_slices_scal(slices,f,irho)
!
      case ('rho') var
        if (lwrite_slice_yz)  slices%yz  = (1.+f(ix_loc,m1:m2, n1:n2,  irho))*spread(rho0z(n1:n2),1,ny)
        if (lwrite_slice_xz)  slices%xz  = (1.+f(l1:l2, iy_loc,n1:n2,  irho))*spread(rho0z(n1:n2),1,nx)
        if (lwrite_slice_xz2) slices%xz2 = (1.+f(l1:l2, iy2_loc,n1:n2,  irho))*spread(rho0z(n1:n2),1,nx)
        if (lwrite_slice_xy)  slices%xy  = (1.+f(l1:l2, m1:m2, iz_loc, irho))*rho0z(iz_loc)
        if (lwrite_slice_xy2) slices%xy2 = (1.+f(l1:l2, m1:m2, iz2_loc,irho))*rho0z(iz2_loc)
        if (lwrite_slice_xy3) slices%xy3 = (1.+f(l1:l2, m1:m2, iz3_loc,irho))*rho0z(iz3_loc)
        if (lwrite_slice_xy4) slices%xy4 = (1.+f(l1:l2, m1:m2, iz4_loc,irho))*rho0z(iz4_loc)
        slices%ready = .true.
!
      endselect var
!
    endsubroutine get_slices_density
!***********************************************************************
    subroutine dynamical_diffusion(uc)
!
!  Dynamically set mass diffusion coefficient given fixed mesh Reynolds number.
!
!  28-feb-13/ccyang: coded
!
!  Input Argument
!      uc
!          Characteristic velocity of the system.
!
      real, intent(in) :: uc
!
!  Hyper-diffusion coefficient
!
      if (ldiff_hyper3_mesh) diffrho_hyper3_mesh = pi5_1 * uc / re_mesh / sqrt(3.0)
!
    endsubroutine dynamical_diffusion
!***********************************************************************
!***********************************************************************
!
!  LOCAL ROUTINES GO BELOW HERE.
!
!***********************************************************************
!***********************************************************************
    real function total_mass(f)
!
!  Gets the total mass.
!
!  18-jun-14/ccyang: coded
!
      use Mpicomm, only: mpiallreduce_sum
!
      real, dimension(mx,my,mz,mfarray), intent(inout) :: f
!
      real :: mass_loc, a
!
      mass_loc = 0.0
      zscan: do n = n1, n2
        a = rho0z(n) * zprim(n)
        do m = m1, m2
          mass_loc = mass_loc + sum(a * yprim(m) * xprim(l1:l2) * (1.0 + f(l1:l2,m,n,irho)))
        enddo
      enddo zscan
      call mpiallreduce_sum(mass_loc, total_mass)
!
    endfunction total_mass
!***********************************************************************
!***********************************************************************
!
!  DUMMY BUT PUBLIC ROUTINES GO BELOW HERE.
!
!***********************************************************************
!***********************************************************************
    subroutine anelastic_after_mn(f, p, df, mass_per_proc)
!
!  Dummy
!
      real, dimension(mx,my,mz,mfarray), intent(in) :: f
      real, dimension(mx,my,mz,mvar), intent(in) :: df
      real, dimension(1), intent(in) :: mass_per_proc
      type(pencil_case), intent(in) :: p
!
      call keep_compiler_quiet(f, df)
      call keep_compiler_quiet(p)
      call keep_compiler_quiet(mass_per_proc)
!
    endsubroutine anelastic_after_mn
!***********************************************************************
    subroutine boussinesq(f)
!
!  dummy routine for the Boussinesq approximation
!
      real, dimension(mx,my,mz,mfarray), intent(in) :: f
!
      call keep_compiler_quiet(f)
!
    endsubroutine boussinesq
!***********************************************************************
    subroutine density_after_boundary(f)
!
      real, dimension(mx,my,mz,mfarray), intent(in) :: f
!
      call keep_compiler_quiet(f)
!
  endsubroutine density_after_boundary
!***********************************************************************
    subroutine density_before_boundary(f)
!
      real, dimension(mx,my,mz,mfarray), intent(in) :: f
!
      call keep_compiler_quiet(f)
!
    endsubroutine density_before_boundary
!***********************************************************************
    subroutine get_init_average_density(f, init_average_density)
!
!  10-dec-09/piyali: added to pass initial average density
!
    real, dimension(mx,my,mz,mfarray), intent(in) :: f
    real, intent(in) :: init_average_density
!
      call keep_compiler_quiet(f)
      call keep_compiler_quiet(init_average_density)
!
    endsubroutine get_init_average_density
!***********************************************************************
    subroutine get_slices_pressure(f, slices)
!
      real, dimension(mx,my,mz,mfarray), intent(in) :: f
      type(slice_data), intent(in) :: slices
!
      call keep_compiler_quiet(f)
      call keep_compiler_quiet(slices%ready)
!
    endsubroutine get_slices_pressure
!***********************************************************************
    real function mean_density(f)
!
!  Calculate mean density of the whole box.
!
!  06-mar-15/ccyang: adapted.
!
      use Mpicomm, only: mpiallreduce_sum
!
      real, dimension(mx,my,mz,mfarray), intent(in) :: f
!
      integer :: m, n
      real :: tmp
!
      mean_density = 0.0
!
      zscan: do n = n1, n2
        tmp = 0.0
        yscan: do m = m1, m2
          tmp = tmp + sum((1.0 + f(l1:l2,m,n,irho)) * dVol_x(l1:l2)) * dVol_y(m)
        enddo yscan
        mean_density = mean_density + tmp * rho0z(n) * dVol_z(n)
      enddo zscan
!
      mean_density = mean_density / box_volume
!
      comm: if (ncpus > 1) then
        call mpiallreduce_sum(mean_density, tmp)
        mean_density = tmp
      endif comm
!
    endfunction mean_density
!***********************************************************************
    subroutine update_char_vel_density(f)
!
!  Updates characteristic velocity for slope-limited diffusion.
!
!  21-oct-15/MR: coded
!
      real, dimension(mx,my,mz,mfarray), intent(INOUT) :: f
!
      if (lslope_limit_diff) call fatal_error('update_char_vel_density', 'not implemented')
!
    endsubroutine update_char_vel_density
!***********************************************************************
    subroutine impose_density_ceiling(f)
!
!  Dummy routine.
!
      real, dimension (mx,my,mz,mfarray), intent(inout) :: f

      call keep_compiler_quiet(f)

    endsubroutine impose_density_ceiling
!***********************************************************************
    subroutine density_after_timestep(f,df,dtsub)
!
      real, dimension(mx,my,mz,mfarray) :: f
      real, dimension(mx,my,mz,mvar) :: df
      real :: dtsub
!
      call keep_compiler_quiet(f,df)
      call keep_compiler_quiet(dtsub)
!
    endsubroutine density_after_timestep
!***********************************************************************

endmodule Density