! $Id$
!
!  This module is used both for the initial condition and during run time.
!  It contains dlnrhon_dt and init_lnrhon, among other auxiliary routines.
!
!** 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 :: lneutraldensity = .true.
!
! MVAR CONTRIBUTION 1
! MAUX CONTRIBUTION 0
!
! PENCILS PROVIDED lnrhon; rhon; rhon1; glnrhon(3); grhon(3);
! PENCILS PROVIDED unglnrhon; ungrhon; del2lnrhon; del2rhon; glnrhon2
! PENCILS PROVIDED del6lnrhon; del6rhon; snglnrhon(3); alpha; zeta
!
!***************************************************************
module NeutralDensity
!
  use Cparam
  use Cdata
  use Messages
!
  implicit none
!
  include 'neutraldensity.h'
!
  real :: kx_lnrhon=1.,ky_lnrhon=1.,kz_lnrhon=1.
  real :: ampllnrhon=0.,rhon_left=1.,rhon_right=1.
  real :: diffrhon=0.,diffrhon_hyper3=0.,diffrhon_shock=0.
  real :: lnrhon_const=0., rhon_const=1.
  real :: lnrhon_int=0.,lnrhon_ext=0.
  real :: lnrhon0,lnrhon_left,lnrhon_right,alpha,zeta
  real, dimension(3) :: diffrhon_hyper3_aniso=0.
  integer, parameter :: ndiff_max=4
  logical :: lmass_source=.false.,lcontinuity_neutral=.true.
  logical :: lupw_lnrhon=.false.,lupw_rhon=.false.
  logical :: ldiffn_normal=.false.,ldiffn_hyper3=.false.,ldiffn_shock=.false.
  logical :: ldiffn_hyper3lnrhon=.false.,ldiffn_hyper3_aniso=.false.
  logical :: lfreeze_lnrhonint=.false.,lfreeze_lnrhonext=.false.
  logical :: ldiffn_hyper3_polar=.false.
  logical :: lpretend_star,lramp_up
  real :: star_form_threshold=1.,star_form_exponent=1.5
  character (len=labellen), dimension(ninit) :: initlnrhon='nothing'
  character (len=labellen), dimension(ndiff_max) :: idiffn=''
  character (len=labellen) :: borderlnrhon='nothing'
  character (len=intlen) :: iinit_str
  integer :: iglobal_rhon=0
!
  namelist /neutraldensity_init_pars/ &
       ampllnrhon,initlnrhon,    &
       rhon_left,rhon_right,lnrhon_const,rhon_const, &
       idiffn,lneutraldensity_nolog,    &
       lcontinuity_neutral,lnrhon0,lnrhon_left,lnrhon_right, &
       alpha,zeta,kx_lnrhon,ky_lnrhon,kz_lnrhon,lpretend_star,&
       star_form_threshold,lramp_up,star_form_exponent
!
  namelist /neutraldensity_run_pars/ &
       diffrhon,diffrhon_hyper3,diffrhon_shock,   &
       lupw_lnrhon,lupw_rhon,idiffn,     &
       lnrhon_int,lnrhon_ext, &
       lfreeze_lnrhonint,lfreeze_lnrhonext,         &
       lnrhon_const,lcontinuity_neutral,borderlnrhon,    &
       diffrhon_hyper3_aniso,alpha,zeta,lpretend_star, &
       star_form_threshold,lramp_up,star_form_exponent
!
!  Diagnostic variables (needs to be consistent with reset list below).
!
  integer :: idiag_rhonm=0,idiag_rhon2m=0,idiag_lnrhon2m=0
  integer :: idiag_rhonmin=0,idiag_rhonmax=0,idiag_unglnrhonm=0
  integer :: idiag_lnrhonmphi=0,idiag_rhonmphi=0,idiag_dtnd=0
  integer :: idiag_rhonmz=0, idiag_rhonmy=0, idiag_rhonmx=0
  integer :: idiag_rhonmxy=0, idiag_rhonmr=0
  integer :: idiag_neutralmass=0
!
  contains
!***********************************************************************
    subroutine register_neutraldensity()
!
!  Initialise variables which should know that we solve the
!  compressible hydro equations: ilnrhon; increase nvar accordingly.
!
!  28-feb-07/wlad: adapted from density
!
      use FArrayManager
!
      if (.not.lcartesian_coords) &
          call fatal_error('register_neutraldensity','non cartesian '//&
           'not yet implemented in the neutrals module')
!
      if (lneutraldensity_nolog) then
        call farray_register_pde('rhon',irhon)
        ilnrhon=irhon
      else
        call farray_register_pde('lnrhon',ilnrhon)
      endif

!
!  Identify version number (generated automatically by SVN).
!
      if (lroot) call svn_id( &
           "$Id$")
!
    endsubroutine register_neutraldensity
!***********************************************************************
    subroutine initialize_neutraldensity()
!
!  Perform any post-parameter-read initialization i.e. calculate derived
!  parameters.
!
!  For compatibility with other applications, we keep the possibility
!  of giving diffrhon units of dxmin*cs0, but cs0 is not well defined general
!
!  28-feb-07/wlad: adapted
!
!
      use BorderProfiles, only: request_border_driving
      use FArrayManager
!
      integer :: i
      logical :: lnothing
!
!  Turn off continuity equation term for 0-D runs.
!
      if (nxgrid*nygrid*nzgrid==1) then
        lcontinuity_neutral=.false.
        print*, 'initialize_neutraldensity: 0-D run, turned off continuity equation'
      endif
!
!  Initialize dust diffusion
!
      ldiffn_normal=.false.
      ldiffn_shock=.false.
      ldiffn_hyper3=.false.
      ldiffn_hyper3lnrhon=.false.
      ldiffn_hyper3_aniso=.false.
      ldiffn_hyper3_polar=.false.
!
      lnothing=.false.
!
      do i=1,ndiff_max
        select case (idiffn(i))
        case ('normal')
          if (lroot) print*,'diffusion: div(D*grad(rhon))'
          ldiffn_normal=.true.
        case ('hyper3')
          if (lroot) print*,'diffusion: (d^6/dx^6+d^6/dy^6+d^6/dz^6)rhon'
          ldiffn_hyper3=.true.
        case ('hyper3lnrhon')
          if (lroot) print*,'diffusion: (d^6/dx^6+d^6/dy^6+d^6/dz^6)lnrhon'
          ldiffn_hyper3lnrhon=.true.
       case ('hyper3_aniso')
          if (lroot) print*,'diffusion: (Dx*d^6/dx^6 + Dy*d^6/dy^6 + Dz*d^6/dz^6)rhon'
          ldiffn_hyper3_aniso=.true.
        case ('hyper3_cyl','hyper3-cyl','hyper3_sph','hyper3-sph')
          if (lroot) print*,'diffusion: Dhyper/pi^4 *(Delta(rhon))^6/Deltaq^2'
          ldiffn_hyper3_polar=.true.
        case ('shock')
          if (lroot) print*,'diffusion: shock diffusion'
          ldiffn_shock=.true.
          call fatal_error("shock diffusion assumes the ions velocity",&
               "not yet implemented for neutrals")
        case ('')
          if (lroot .and. (.not. lnothing)) print*,'diffusion: nothing'
        case default
          write(unit=errormsg,fmt=*) 'initialize_neutraldensity: ', &
              'No such value for idiff(',i,'): ', trim(idiffn(i))
          call fatal_error('initialize_neutraldensity',errormsg)
        endselect
        lnothing=.true.
      enddo
!
!  If we're timestepping, die or warn if the the diffusion coefficient that
!  corresponds to the chosen diffusion type is not set.
!
      if (lrun) then
        if (ldiffn_normal.and.diffrhon==0.0) &
            call warning('initialize_neutraldensity', &
            'Diffusion coefficient diffrhon is zero!')
        if ( (ldiffn_hyper3 .or. ldiffn_hyper3lnrhon) &
            .and. diffrhon_hyper3==0.0) &
            call fatal_error('initialize_neutraldensity', &
            'Diffusion coefficient diffrhon_hyper3 is zero!')
        if ( (ldiffn_hyper3_aniso) .and.  &
             ((diffrhon_hyper3_aniso(1)==0. .and. nxgrid/=1 ).or. &
              (diffrhon_hyper3_aniso(2)==0. .and. nygrid/=1 ).or. &
              (diffrhon_hyper3_aniso(3)==0. .and. nzgrid/=1 )) ) &
            call fatal_error('initialize_neutraldensity', &
            'A diffusion coefficient of diffrhon_hyper3 is zero!')
        if (ldiffn_shock .and. diffrhon_shock==0.0) &
            call fatal_error('initialize_neutraldensity', &
            'diffusion coefficient diffrhon_shock is zero!')
      endif
!
!  Hyperdiffusion only works with (not log) density. One must either use
!  ldensity_nolog=T or work with GLOBAL =   global_nolog_density.
!
      if ((ldiffn_hyper3.or.ldiffn_hyper3_aniso) .and. (.not. lneutraldensity_nolog)) then
        if (lroot) print*,"initialize_neutraldensity: Creating global array for rhon to use hyperdiffusion"
        call farray_register_global('rhon',iglobal_rhon)
      endif
!
      if (lfreeze_lnrhonint) lfreeze_varint(ilnrhon) = .true.
      if (lfreeze_lnrhonext) lfreeze_varext(ilnrhon) = .true.
!
!  Tell the BorderProfiles module if we intend to use border driving, so
!  that the module can request the right pencils.
!
      if (borderlnrhon/='nothing') call request_border_driving(borderlnrhon)
!
    endsubroutine initialize_neutraldensity
!***********************************************************************
    subroutine read_neutraldensity_init_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=neutraldensity_init_pars, IOSTAT=iostat)
!
    endsubroutine read_neutraldensity_init_pars
!***********************************************************************
    subroutine write_neutraldensity_init_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=neutraldensity_init_pars)
!
    endsubroutine write_neutraldensity_init_pars
!***********************************************************************
    subroutine read_neutraldensity_run_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=neutraldensity_run_pars, IOSTAT=iostat)
!
    endsubroutine read_neutraldensity_run_pars
!***********************************************************************
    subroutine write_neutraldensity_run_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=neutraldensity_run_pars)
!
    endsubroutine write_neutraldensity_run_pars
!***********************************************************************
    subroutine init_lnrhon(f)
!
!  initialise lnrhon; called from start.f90
!
!  28-feb-07/wlad: adapted
!
      use General, only: itoa,complex_phase
      use Initcond
      use InitialCondition, only: initial_condition_lnrhon
      use General, only: notanumber
      use EquationOfState, only: cs20, cs2bot,cs2top
!
      real, dimension (mx,my,mz,mfarray) :: f
!
      real :: zbot, ztop
      logical :: lnothing
      integer :: j
!
!  define bottom and top height
!
      zbot=xyz0(3)
      ztop=xyz0(3)+Lxyz(3)
!
!  Set default values for sound speed at top and bottom.
!  These may be updated in one of the following initialization routines.
!
      cs2top=cs20; cs2bot=cs20
!
!  different initializations of lnrhon (called from start).
!  If initrhon does't match, f=0 is assumed (default).
!
      lnothing=.true.

      do j=1,ninit

         if (initlnrhon(j)/='nothing') then

            lnothing=.false.

            iinit_str=itoa(j)

            select case (initlnrhon(j))
!
! some one-liners from density
!
            case ('zero', '0'); f(:,:,:,ilnrhon)=0.
            case ('const_lnrhon'); f(:,:,:,ilnrhon)=lnrhon_const
            case ('const_rhon'); f(:,:,:,ilnrhon)=log(rhon_const)
            case ('constant'); f(:,:,:,ilnrhon)=log(rhon_left)
            case ('scale-ions')
              if (ldensity_nolog) then
                f(:,:,:,ilnrhon)=log(rhon_const)+log(f(:,:,:,ilnrho))
              else
                f(:,:,:,ilnrhon)=log(rhon_const)+f(:,:,:,ilnrho)
              endif
            case ('sinwave-z'); call sinwave(ampllnrhon,f,ilnrhon,kz=kz_lnrhon)
            case ('gaussian-noise')
               if (lnrhon_left /= 0.) f(:,:,:,ilnrhon)=lnrhon_left
               call gaunoise(ampllnrhon,f,ilnrhon,ilnrhon)
            case default
               !
               !  Catch unknown values
               !
               write(unit=errormsg,fmt=*) 'No such value for initlnrhon(' &
                    //trim(iinit_str)//'): ',trim(initlnrhon(j))
               call fatal_error('init_lnrhon',errormsg)
!
            endselect

            if (lroot) print*,'init_lnrhon: initlnrhon(' &
                 //trim(iinit_str)//') = ',trim(initlnrhon(j))

         endif

      enddo
!
!  Interface for user's own initial conditon
!
      if (linitial_condition) call initial_condition_lnrhon(f)
!
      if (lnothing.and.lroot) print*,'init_lnrhon: nothing'
!
!  If unlogarithmic density considered, take exp of lnrhon resulting from
!  initlnrhon
!
      if (lneutraldensity_nolog) f(:,:,:,irhon)=exp(f(:,:,:,ilnrhon))
!
!  sanity check
!
      if (notanumber(f(l1:l2,m1:m2,n1:n2,ilnrhon))) then
        call error('init_lnrhon', 'Imaginary density values')
      endif
!
    endsubroutine init_lnrhon
!***********************************************************************
    subroutine pencil_criteria_neutraldensity()
!
!  All pencils that the NeutralDensity module depends on are specified here.
!
!  28-feb-07/wlad: adapted
!
!  always needed for ionization and recombination
!
      lpenc_requested(i_rho)  =.true.
      lpenc_requested(i_rhon) =.true.
      lpenc_requested(i_alpha)=.true.
      lpenc_requested(i_zeta)=.true.
!
      if (.not.lneutraldensity_nolog) then
        lpenc_requested(i_rho1) =.true.
        lpenc_requested(i_rhon1)=.true.
      endif
!
      if (lpretend_star) then
        lpenc_requested(i_rho1)=.true.
        lpenc_requested(i_uu)=.true.
        lpenc_requested(i_rcyl_mn1)=.true.
      endif
!
      if (lcontinuity_neutral) then
        lpenc_requested(i_divun)=.true.
        if (lneutraldensity_nolog) then
          lpenc_requested(i_ungrhon)=.true.
        else
          lpenc_requested(i_unglnrhon)=.true.
        endif
      endif
      if (ldiffn_shock) then
        lpenc_requested(i_shock)=.true.
        lpenc_requested(i_gshock)=.true.
        if (lneutraldensity_nolog) then
          lpenc_requested(i_grhon)=.true.
          lpenc_requested(i_del2rhon)=.true.
        else
          lpenc_requested(i_glnrhon)=.true.
          lpenc_requested(i_glnrhon2)=.true.
          lpenc_requested(i_del2lnrhon)=.true.
        endif
      endif
      if (ldiffn_normal) then
        if (lneutraldensity_nolog) then
          lpenc_requested(i_del2rhon)=.true.
        else
          lpenc_requested(i_glnrhon2)=.true.
          lpenc_requested(i_del2lnrhon)=.true.
        endif
      endif
      if (ldiffn_hyper3.or.ldiffn_hyper3_aniso) &
           lpenc_requested(i_del6rhon)=.true.
      if (ldiffn_hyper3.and..not.lneutraldensity_nolog) &
           lpenc_requested(i_rhon)=.true.
      if (ldiffn_hyper3lnrhon) &
           lpenc_requested(i_del6lnrhon)=.true.
      if (ldiffn_hyper3_polar.and..not.lneutraldensity_nolog) &
           lpenc_requested(i_rhon1)=.true.
!
      if (lmass_source) lpenc_requested(i_rcyl_mn)=.true.
!
      lpenc_diagnos2d(i_lnrhon)=.true.
      lpenc_diagnos2d(i_rhon)=.true.
!
      if (idiag_rhonm/=0 .or. idiag_rhonmz/=0 .or. idiag_rhonmy/=0 .or. &
           idiag_rhonmx/=0 .or. idiag_rhon2m/=0 .or. idiag_rhonmin/=0 .or. &
           idiag_rhonmax/=0 .or. idiag_rhonmxy/=0 .or. idiag_neutralmass/=0) &
           lpenc_diagnos(i_rhon)=.true.
      if (idiag_lnrhon2m/=0) lpenc_diagnos(i_lnrhon)=.true.
      if (idiag_unglnrhonm/=0) lpenc_diagnos(i_unglnrhon)=.true.
!
    endsubroutine pencil_criteria_neutraldensity
!***********************************************************************
    subroutine pencil_interdep_neutraldensity(lpencil_in)
!
!  Interdependency among pencils from the NeutralDensity module is
!    specified here.
!
!  28-feb-07/wlad: adapted
!
      logical, dimension(npencils) :: lpencil_in
!
      if (lneutraldensity_nolog) then
        if (lpencil_in(i_rhon1)) lpencil_in(i_rhon)=.true.
      else
        if (lpencil_in(i_rhon)) lpencil_in(i_rhon1)=.true.
      endif
      if (lpencil_in(i_unglnrhon)) then
        lpencil_in(i_uun)=.true.
        lpencil_in(i_glnrhon)=.true.
      endif
      if (lpencil_in(i_ungrhon)) then
        lpencil_in(i_uun)=.true.
        lpencil_in(i_grhon)=.true.
      endif
      if (lpencil_in(i_glnrhon2)) lpencil_in(i_glnrhon)=.true.
      if (lpencil_in(i_snglnrhon)) then
        lpencil_in(i_snij)=.true.
        lpencil_in(i_glnrhon)=.true.
      endif
!  The pencils glnrhon and grhon come in a bundle.
      if (lpencil_in(i_glnrhon) .and. lpencil_in(i_grhon)) then
        if (lneutraldensity_nolog) then
          lpencil_in(i_grhon)=.false.
        else
          lpencil_in(i_glnrhon)=.false.
        endif
      endif
!
    endsubroutine pencil_interdep_neutraldensity
!***********************************************************************
    subroutine calc_pencils_neutraldensity(f,p)
!
!  Calculate NeutralDensity pencils.
!  Most basic pencils should come first, as others may depend on them.
!
!  28-feb-07/wlad: adapted
!
      use Sub, only: grad,dot,dot2,u_dot_grad,del2,del6,multmv,g2ij
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx) :: tmp,smooth_step_threshold
      real :: alpha_time,ramping_period
      type (pencil_case) :: p
!
      integer :: i, mm, nn
!
      intent(inout) :: f,p
! lnrhon
      if (lpencil(i_lnrhon)) then
        if (lneutraldensity_nolog) then
          p%lnrhon=log(f(l1:l2,m,n,irhon))
        else
          p%lnrhon=f(l1:l2,m,n,ilnrhon)
        endif
      endif
! rhon1 and rhon
      if (lneutraldensity_nolog) then
        if (lpencil(i_rhon)) p%rhon=f(l1:l2,m,n,irhon)
        if (lpencil(i_rhon1)) p%rhon1=1.0/p%rhon
      else
        if (lpencil(i_rhon1)) p%rhon1=exp(-f(l1:l2,m,n,ilnrhon))
        if (lpencil(i_rhon)) p%rhon=1.0/p%rhon1
      endif
! glnrhon and grhon
      if (lpencil(i_glnrhon).or.lpencil(i_grhon)) then
        if (lneutraldensity_nolog) then
          call grad(f,ilnrhon,p%grhon)
          if (lpencil(i_glnrhon)) then
            do i=1,3
              p%glnrhon(:,i)=p%grhon(:,i)/p%rhon
            enddo
          endif
        else
          call grad(f,ilnrhon,p%glnrhon)
          if (lpencil(i_grhon)) then
            do i=1,3
              p%grhon(:,i)=p%rhon*p%glnrhon(:,i)
            enddo
          endif
        endif
      endif
! unglnrhon
      if (lpencil(i_unglnrhon)) then
        if (lneutraldensity_nolog) then
          call dot(p%uun,p%glnrhon,p%unglnrhon)
        else
          call u_dot_grad(f,ilnrhon,p%glnrhon,p%uun,p%unglnrhon,UPWIND=lupw_lnrhon)
        endif
      endif
! ungrhon
      if (lpencil(i_ungrhon)) then
        if (lneutraldensity_nolog) then
          call u_dot_grad(f,ilnrhon,p%grhon,p%uun,p%ungrhon,UPWIND=lupw_rhon)
        else
          call dot(p%uun,p%grhon,p%ungrhon)
        endif
      endif
! glnrhon2
      if (lpencil(i_glnrhon2)) call dot2(p%glnrhon,p%glnrhon2)
! del2lnrhon
      if (lpencil(i_del2lnrhon)) then
        if (lneutraldensity_nolog) then
          if (headtt) then
            call fatal_error('calc_pencils_neutraldensity', &
                             'del2lnrhon not available for non-logarithmic mass density')
          endif
        else
          call del2(f,ilnrhon,p%del2lnrhon)
        endif
      endif
! del2rhon
      if (lpencil(i_del2rhon)) then
        if (lneutraldensity_nolog) then
          call del2(f,ilnrhon,p%del2rhon)
        else
          if (headtt) then
            call fatal_error('calc_pencils_neutraldensity', &
                             'del2rhon not available for logarithmic mass density')
          endif
        endif
      endif
! del6rhon
      if (lpencil(i_del6rhon)) then
        if (lneutraldensity_nolog) then
          call del6(f,ilnrhon,p%del6rhon)
        else
          if (lfirstpoint) then
            !
            ! Fill global rhon array using the ilnrhon data
!ajwm This won't work unless earlt_finalize is used... ?
            if (iglobal_rhon/=0) then
              do mm=1,my; do nn=1,mz
                f(:,mm,nn,iglobal_rhon) = exp(f(:,mm,nn,ilnrhon))
              enddo; enddo
            else
              call fatal_error("calc_pencils_neutraldensity",&
                       "A global rhon slot must be available for calculating del6rhon from lnrhon")
            endif
          endif
          if (iglobal_rhon/=0) call del6(f,iglobal_rhon,p%del6rhon)
        endif
      endif
! del6lnrhon
      if (lpencil(i_del6lnrhon)) then
        if (lneutraldensity_nolog) then
          if (headtt) then
            call fatal_error('calc_pencils_neutraldensity', &
                             'del6lnrhon not available for non-logarithmic mass density')
          endif
        else
          call del6(f,ilnrhon,p%del6lnrhon)
        endif
      endif
! snglnrhon
      if (lpencil(i_snglnrhon)) call multmv(p%snij,p%glnrhon,p%snglnrhon)
!
! ionization and recombination pencils
!
      p%zeta=zeta
      if (lpretend_star) then
        if (.not.lcylindrical_coords) &
             call fatal_error("lpretend_star",&
             "not implemented for other than cylindrical coordinates")
        !smooth transition over a ramping period of 5 orbits
        if (lramp_up) then
          ramping_period=2*pi*x(lpoint) !omega=v/r; v=1; 1/omega=r
          if (t <= ramping_period) then
            alpha_time=alpha*(sin((.5*pi)*(t/ramping_period))**2)
          else
            alpha_time=alpha
          endif
        else
          alpha_time=alpha
        endif
!
! Star formation rate
! These lines below recover d/dt(rho_star)=sfr_const*omega*rho_gas**1.5
!
! There is threshold that has to be smoothed. The star formation
! rate falls drastically after the threshold of 5-10 solar masses
! per cubic parsec. A arctangent smoothing over a tenth of this value
! is okay to avoid numerical disasters.
!
        tmp=(p%rho-star_form_threshold)/(.1*star_form_threshold)
        smooth_step_threshold=.5*(1+atan(tmp)*2*pi_1)

        !OO=p%uu(*,2)*p%rcyl_mn1
        p%alpha=(alpha_time*smooth_step_threshold)*&
             (p%uu(:,2)*p%rcyl_mn1)*p%rho1**(2-star_form_exponent)
!
      else
        p%alpha=alpha
      endif
!
    endsubroutine calc_pencils_neutraldensity
!***********************************************************************
    subroutine dlnrhon_dt(f,df,p)
!
!  continuity equation
!  calculate dlnrhon/dt = - un.gradlnrhon - divun
!
!  28-feb-07/wlad: adapted
!
      use Deriv, only: der6
      use Diagnostics
      use Mpicomm, only: stop_it
      use Sub, only: identify_bcs, del6fj, dot_mn
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx,my,mz,mvar) :: df
      type (pencil_case) :: p
!
      intent(in)  :: f,p
      intent(out) :: df
!
      real, dimension (nx) :: fdiff, gshockglnrhon, gshockgrhon, tmp, diffus_diffrhon, diffus_diffrhon3
      integer :: j
!
!  identify module and boundary conditions
!
      if (headtt.or.ldebug) print*,'dlnrhon_dt: SOLVE dlnrhon_dt'
      if (headtt) call identify_bcs('lnrhon',ilnrhon)
!
!  continuity equation
!
      if (lcontinuity_neutral) then
        if (lneutraldensity_nolog) then
          df(l1:l2,m,n,irhon)   = df(l1:l2,m,n,irhon)   - p%ungrhon   - p%rhon*p%divun
        else
          df(l1:l2,m,n,ilnrhon) = df(l1:l2,m,n,ilnrhon) - p%unglnrhon - p%divun
        endif
      endif
!
!  Ionization and recombination
!
      if (lneutraldensity_nolog) then
         df(l1:l2,m,n,irhon)   = df(l1:l2,m,n,irhon)   - p%zeta*p%rhon        + p%alpha*p%rho**2
         df(l1:l2,m,n,ilnrho ) = df(l1:l2,m,n,ilnrho ) + p%zeta*p%rhon        - p%alpha*p%rho**2
      else
         df(l1:l2,m,n,ilnrhon) = df(l1:l2,m,n,ilnrhon) - p%zeta               + p%alpha*p%rho**2*p%rhon1
         df(l1:l2,m,n,ilnrho ) = df(l1:l2,m,n,ilnrho ) + p%zeta*p%rhon*p%rho1 - p%alpha*p%rho
      endif
!
!  Mass diffusion
!
      fdiff=0.0
      diffus_diffrhon=0.; diffus_diffrhon3=0.
!
      if (ldiffn_normal) then  ! Normal diffusion operator
        if (lneutraldensity_nolog) then
          fdiff = fdiff + diffrhon*p%del2rhon
        else
          fdiff = fdiff + diffrhon*(p%del2lnrhon+p%glnrhon2)
        endif
        if (lfirst.and.ldt) diffus_diffrhon=diffus_diffrhon+diffrhon*dxyz_2
        if (headtt) print*,'dlnrhon_dt: diffrhon=', diffrhon
      endif
!
!  Hyper diffusion
!
      if (ldiffn_hyper3) then
        if (lneutraldensity_nolog) then
          fdiff = fdiff + diffrhon_hyper3*p%del6rhon
        else
          fdiff = fdiff + 1/p%rhon*diffrhon_hyper3*p%del6rhon
        endif
        if (lfirst.and.ldt) diffus_diffrhon3=diffus_diffrhon3+diffrhon_hyper3*dxyz_6
        if (headtt) print*,'dlnrhon_dt: diffrhon_hyper3=', diffrhon_hyper3
      endif
!
      if (ldiffn_hyper3_aniso) then
        if (lneutraldensity_nolog) then
          call del6fj(f,diffrhon_hyper3_aniso,ilnrhon,tmp)
          fdiff = fdiff + tmp
          if (lfirst.and.ldt) diffus_diffrhon3=diffus_diffrhon3+ &
               diffrhon_hyper3_aniso(1)*dline_1(:,1)**6 + &
               diffrhon_hyper3_aniso(2)*dline_1(:,2)**6 + &
               diffrhon_hyper3_aniso(3)*dline_1(:,3)**6
          if (headtt) &
               print*,'dlnrhon_dt: diffrhon_hyper3=(Dx,Dy,Dz)=',&
               diffrhon_hyper3_aniso
        else
          call stop_it("anisotropic hyperdiffusion not implemented for lnrhon")
        endif
      endif
!
      if (ldiffn_hyper3_polar) then
        do j=1,3
          call der6(f,ilnrhon,tmp,j,IGNOREDX=.true.)
          if (.not.lneutraldensity_nolog) tmp=tmp*p%rhon1
          fdiff = fdiff + diffrhon_hyper3*pi4_1*tmp*dline_1(:,j)**2
        enddo
        if (lfirst.and.ldt) &
             diffus_diffrhon3=diffus_diffrhon3+diffrhon_hyper3*pi4_1*dxmin_pencil**4
        if (headtt) print*,'dlnrhon_dt: diffrhon_hyper3=', diffrhon_hyper3
      endif
!
      if (ldiffn_hyper3lnrhon) then
        if (.not. lneutraldensity_nolog) then
          fdiff = fdiff + diffrhon_hyper3*p%del6lnrhon
        endif
        if (lfirst.and.ldt) diffus_diffrhon3=diffus_diffrhon3+diffrhon_hyper3*dxyz_6
        if (headtt) print*,'dlnrhon_dt: diffrhon_hyper3=', diffrhon_hyper3
      endif
!
!  Shock diffusion
!
      if (ldiffn_shock) then
        if (lneutraldensity_nolog) then
          call dot_mn(p%gshock,p%grhon,gshockgrhon)
          df(l1:l2,m,n,ilnrhon) = df(l1:l2,m,n,ilnrhon) + &
              diffrhon_shock*p%shock*p%del2rhon + diffrhon_shock*gshockgrhon
        else
          call dot_mn(p%gshock,p%glnrhon,gshockglnrhon)
          df(l1:l2,m,n,ilnrhon) = df(l1:l2,m,n,ilnrhon) + &
              diffrhon_shock*p%shock*(p%del2lnrhon+p%glnrhon2) + &
              diffrhon_shock*gshockglnrhon
        endif
        if (lfirst.and.ldt) &
            diffus_diffrhon=diffus_diffrhon+diffrhon_shock*p%shock*dxyz_2
        if (headtt) print*,'dlnrhon_dt: diffrhon_shock=', diffrhon_shock
      endif
!
!  Add diffusion term to continuity equation
!
      if (lneutraldensity_nolog) then
        df(l1:l2,m,n,irhon)   = df(l1:l2,m,n,irhon)   + fdiff
      else
        df(l1:l2,m,n,ilnrhon) = df(l1:l2,m,n,ilnrhon) + fdiff
      endif
!
      if (lfirst.and.ldt) then
        if (headtt.or.ldebug) &
          print*,'dlnrhon_dt: max(diffus_diffrhon) =', maxval(diffus_diffrhon)
        maxdiffus=max(maxdiffus,diffus_diffrhon)
        maxdiffus3=max(maxdiffus3,diffus_diffrhon3)
      endif
!
!  Apply border profile
!
      if (lborder_profiles) call set_border_neutraldensity(f,df,p)
!
!  2d-averages
!  Note that this does not necessarily happen with ldiagnos=.true.
!
      if (l2davgfirst) then
        call zsum_mn_name_xy(p%rhon,idiag_rhonmxy)
        call phisum_mn_name_rz(p%lnrhon,idiag_lnrhonmphi)
        call phisum_mn_name_rz(p%rhon,idiag_rhonmphi)
      endif
!
!  1d-averages. Happens at every it1d timesteps, NOT at every it1
!
      if (l1davgfirst) then
         call phizsum_mn_name_r(p%rhon,idiag_rhonmr)
         call xysum_mn_name_z(p%rhon,idiag_rhonmz)
         call yzsum_mn_name_x(p%rhon,idiag_rhonmx)
         call xzsum_mn_name_y(p%rhon,idiag_rhonmy)
      endif
!
!  Calculate density diagnostics
!
      if (ldiagnos) then
        call sum_mn_name(p%rhon,idiag_rhonm)
        call sum_lim_mn_name(p%rhon,idiag_neutralmass,p)
        if (idiag_rhonmin/=0) &
            call max_mn_name(-p%rhon,idiag_rhonmin,lneg=.true.)
        call max_mn_name(p%rhon,idiag_rhonmax)
        if (idiag_rhon2m/=0)    call sum_mn_name(p%rhon**2,idiag_rhon2m)
        if (idiag_lnrhon2m/=0)  call sum_mn_name(p%lnrhon**2,idiag_lnrhon2m)
        call sum_mn_name(p%unglnrhon,idiag_unglnrhonm)
        if (idiag_dtnd/=0) &
            call max_mn_name(diffus_diffrhon/cdtv,idiag_dtnd,l_dt=.true.)
      endif
!
    endsubroutine dlnrhon_dt
!***********************************************************************
    subroutine set_border_neutraldensity(f,df,p)
!
!  Calculates the driving term for the border profile
!  of the lnrhon variable.
!
!  28-jul-06/wlad: coded
!
      use BorderProfiles, only: border_driving
      use EquationOfState, only: cs0,cs20
!
      real, dimension(mx,my,mz,mfarray) :: f
      type (pencil_case) :: p
      real, dimension(mx,my,mz,mvar) :: df
      real, dimension(nx) :: f_target !,OO_sph,OO_cyl,cs,theta
!      real :: r0_pot=0.1
!
      select case (borderlnrhon)
!
      case ('zero','0')
         f_target=0.
      case ('constant')
         f_target=lnrhon_const
      case ('stratification')
         !OO_sph = sqrt((r_mn**2 + r0_pot**2)**(-1.5))
         !OO_cyl = sqrt((rcyl_mn**2 + r0_pot**2)**(-1.5))
         !cs = OO_cyl*rcyl_mn*cs0
         !f_target=lnrhon_const - 0.5*(theta/cs0)**2
         f_target=(p%rcyl_mn-p%r_mn)/(cs20*p%r_mn)
      case ('nothing')
         if (lroot.and.ip<=5) &
              print*,"set_border_neutraldensity: borderlnrhon='nothing'"
      case default
         write(unit=errormsg,fmt=*) &
              'set_border_neutraldensity: No such value for borderlnrhon: ', &
              trim(borderlnrhon)
         call fatal_error('set_border_neutraldensity',errormsg)
      endselect
!
      if (lneutraldensity_nolog) f_target=exp(f_target)
!
      if (borderlnrhon /= 'nothing') then
        call border_driving(f,df,p,f_target,ilnrhon)
      endif
!
    endsubroutine set_border_neutraldensity
!***********************************************************************
    subroutine rprint_neutraldensity(lreset,lwrite)
!
!  reads and registers print parameters relevant for compressible part
!
!  28-feb-07/wlad: adapted
!
      use Diagnostics, only: parse_name
      use FArrayManager, only: farray_index_append
!
      logical :: lreset
      logical, optional :: lwrite
!
      integer :: iname, inamex, inamey, inamez, inamexy, irz, inamer
      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_rhonm=0; idiag_rhon2m=0; idiag_lnrhon2m=0; idiag_unglnrhonm=0
        idiag_rhonmin=0; idiag_rhonmax=0; idiag_dtnd=0
        idiag_lnrhonmphi=0; idiag_rhonmphi=0
        idiag_rhonmz=0; idiag_rhonmy=0; idiag_rhonmx=0
        idiag_rhonmxy=0; idiag_rhonmr=0; idiag_neutralmass=0
        diffrhon=0.
      endif
!
!  iname runs through all possible names that may be listed in print.in
!
      if (lroot.and.ip<14) print*,'rprint_neutraldensity: run through parse list'
      do iname=1,nname
        call parse_name(iname,cname(iname),cform(iname),'rhonm',idiag_rhonm)
        call parse_name(iname,cname(iname),cform(iname),'rhon2m',idiag_rhon2m)
        call parse_name(iname,cname(iname),cform(iname),'rhonmin',idiag_rhonmin)
        call parse_name(iname,cname(iname),cform(iname),'rhonmax',idiag_rhonmax)
        call parse_name(iname,cname(iname),cform(iname),'lnrhon2m',idiag_lnrhon2m)
        call parse_name(iname,cname(iname),cform(iname),'unglnrhonm',idiag_unglnrhonm)
        call parse_name(iname,cname(iname),cform(iname),'dtnd',idiag_dtnd)
        call parse_name(iname,cname(iname),cform(iname),'neutralmass',idiag_neutralmass)
      enddo
!
!  check for those quantities for which we want xy-averages
!
      do inamez=1,nnamez
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'rhonmz',idiag_rhonmz)
      enddo
!
!  check for those quantities for which we want xz-averages
!
      do inamey=1,nnamey
        call parse_name(inamey,cnamey(inamey),cformy(inamey),'rhonmy',idiag_rhonmy)
      enddo
!
!  check for those quantities for which we want yz-averages
!
      do inamex=1,nnamex
        call parse_name(inamex,cnamex(inamex),cformx(inamex),'rhonmx',idiag_rhonmx)
      enddo
!
!  check for those quantities for which we want phiz-averages
!
      do inamer=1,nnamer
        call parse_name(inamer,cnamer(inamer),cformr(inamer),'rhonmr',idiag_rhonmr)
      enddo
!
!  check for those quantities for which we want z-averages
!
      do inamexy=1,nnamexy
        call parse_name(inamexy,cnamexy(inamexy),cformxy(inamexy),'rhonmxy',idiag_rhonmxy)
      enddo
!
!  check for those quantities for which we want phi-averages
!
      do irz=1,nnamerz
        call parse_name(irz,cnamerz(irz),cformrz(irz),&
            'lnrhonmphi',idiag_lnrhonmphi)
        call parse_name(irz,cnamerz(irz),cformrz(irz),'rhonmphi',idiag_rhonmphi)
      enddo
!
!  write column where which density variable is stored
!
      if (lwr) then
        call farray_index_append('ilnrhon',ilnrhon)
      endif
!
    endsubroutine rprint_neutraldensity
!***********************************************************************
endmodule NeutralDensity