! $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; 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., luse_as_ionization=.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', alpha_prescription='const'
  character (len=intlen) :: iinit_str
!
  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, &
       alpha_prescription, luse_as_ionization
!
  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, &
       alpha_prescription, luse_as_ionization
!
!  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, idiag_alprec=0
!
! Auxiliaries
!
  real :: alpha_time
  real, dimension(nx) :: diffus_diffrhon

  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 not_implemented('register_neutraldensity','non-Cartesian coords')
!
      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 General, only: itoa
      use BorderProfiles, only: request_border_driving
      use FArrayManager
!
      integer :: i
      logical :: lnothing
!
!  Turn off continuity equation term for 0-D runs.
!
      if (dimensionality==0) 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 (lneutraldensity_nolog) call fatal_error('initialize_neutraldensity', &
                             "'hyper3lnrhon' diffusion not meaningful for linear density")
          if (lroot) print*,'diffusion: (d^6/dx^6+d^6/dy^6+d^6/dz^6)lnrhon'
          ldiffn_hyper3lnrhon=.true.
       case ('hyper3_aniso')
          if (.not.lneutraldensity_nolog) &
            call not_implemented('initialize_neutraldensity',"anisotropic hyperdiffusion for log density")
          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 not_implemented("initialize_neutraldensity", &
                               "shock diffusion assumes ions velocity for neutrals")
        case ('')
          if (lroot .and. (.not. lnothing)) print*,'diffusion: nothing'
        case default
          call fatal_error('initialize_neutraldensity', &
               'No such value for idiff('//trim(itoa(i))//'): '//trim(idiffn(i)))
        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

      if (.not.lneutraldensity_nolog) then         ! for logarithmic density

        if (ldiffn_hyper3.or.ldiffn_hyper3_aniso) then
!
!  To obtain del6rhon with log density, one must work with a COMAUX slot in the f-array for linear density.
!
          if (lroot) print*,"initialize_neutraldensity: Creating comm. aux variable rhon for del6rhon ", &
                            "and setting exponentiation BC for rhon"
          if (irhon==0) call farray_register_auxiliary('rhon',irhon,communicated=.true.)

          if (.not.lperi(1)) then
            bcx12(irhon,:)='exp'
            fbcx(irhon,:) = float(ilnrhon)
          endif
          if (.not.lperi(2)) then
            bcy12(irhon,:)='exp'
            fbcy(irhon,:) = float(ilnrhon)
          endif
          if (.not.lperi(3)) then
            bcz12(irhon,:)='exp'
            fbcz(irhon,:) = float(ilnrhon)
          endif

        endif
      endif

      if (lpretend_star.and..not.lcylindrical_coords) &
        call not_implemented('initialize_neutraldensity', &
                             "lpretend_star for other than cylindrical coordinates")
!
      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.
!
      select case (borderlnrhon)
!
      case ('zero','0','constant','stratification')
         call request_border_driving(borderlnrhon)
      case ('nothing')
        if (lroot.and.ip<=5) print*,"initialize_neutraldensity: borderlnrhon='nothing'"
      case default
         call fatal_error('initialize_neutraldensity','no such borderlnrhon: '//trim(borderlnrhon))
      endselect
!
    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
!
      logical :: lnothing
      integer :: j
!
!  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
               call fatal_error('init_lnrhon','no such initlnrhon('// &
                                trim(iinit_str)//'): '//trim(initlnrhon(j)))
!
            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
!  initialisation above. (irhon=ilnrhon!)
!
      if (lneutraldensity_nolog) f(:,:,:,irhon)=exp(f(:,:,:,ilnrhon))
!
!  sanity check
!
      if (notanumber(f(l1:l2,m1:m2,n1:n2,ilnrhon))) call error('init_lnrhon','imaginary density values')
!
    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.
!
!  needed for temperature prescription
!
      if (alpha_prescription=='Temp_dep') then
        lpenc_requested(i_TT)=.true.
      endif

      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.
        endif
      endif
      if (ldiffn_normal) then
        if (lneutraldensity_nolog) then
          lpenc_requested(i_del2rhon)=.true.
        else
          lpenc_requested(i_glnrhon2)=.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.
      if (idiag_alprec/=0) lpenc_diagnos(i_alpha)=.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
      type (pencil_case) :: p
      intent(inout) :: f,p
!
      real, dimension (nx) :: tmp,smooth_step_threshold
      real :: unit_alpha=impossible
      integer :: i
!
! 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)
! del2rhon
      if (lpencil(i_del2rhon)) then
        call del2(f,ilnrhon,p%del2rhon)    ! can be either del2rhon or del2lnrhon
        if (.not.lneutraldensity_nolog) p%del2rhon=p%del2rhon+p%glnrhon2  ! here del2rhon/rhon
      endif
!
! del6rhon
! Note that irhon either points to the linear neutraldensity as a PDE variable
! for lneutraldensity_nolog, or as a COMAUX variable otherwise.
!
      if (lpencil(i_del6rhon).and.irhon>0) call del6(f,irhon,p%del6rhon)
! del6lnrhon
      if (lpencil(i_del6lnrhon)) call del6(f,ilnrhon,p%del6lnrhon)
! snglnrhon
      if (lpencil(i_snglnrhon)) call multmv(p%snij,p%glnrhon,p%snglnrhon)
!
! ionization and recombination pencils
!
      p%zeta=zeta
      if (lpretend_star) then
!
! Star formation rate
! These lines below recover d/dt(rho_star)=sfr_const*omega*rho_gas**1.5
!
! There is a threshold that has to be smoothed. The star formation
! rate falls drastically after the threshold of 5-10 solar masses
! per cubic parsec. An arctan 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
!
!  Choice of different alpha prescriptions
!
        select case (alpha_prescription)
        case ('const'); p%alpha=alpha
        case ('Temp_dep')
          if (lentropy) then
            unit_alpha=1./(unit_time*unit_density)
            p%alpha=unit_alpha*1.14e-19*4.309e6*p%TT**(-.6166)/(1.+.6703*p%TT**.53)
  !         if (ip<10) print*,'AXEL: p%alpha=',p%alpha
  !         if (ip<10) print*,'AXEL: unit_alpha=',unit_alpha
          else
            call fatal_error('calc_pencils_neutraldensity', &
              'no energy equation is used')
          endif
        case default
          call fatal_error('calc_pencils_neutraldensity', &
            'No such value for alpha_prescription')
        endselect
      endif
!
    endsubroutine calc_pencils_neutraldensity
!***********************************************************************
    subroutine neutraldensity_after_boundary(f)

      use General, only: keep_compiler_quiet

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

      if (lpretend_star) then
!
! Smooth transition over a ramping period of 5 orbits.
!
        if (lramp_up) then              ! MR: strange use of lpoint!
          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

      endif

      call keep_compiler_quiet(f)

    endsubroutine neutraldensity_after_boundary
!***********************************************************************
    subroutine neutraldensity_before_boundary(f)

      use General, only: keep_compiler_quiet

      real, dimension (mx,my,mz,mfarray) :: f
!
! Fill global rhon array using the ilnrhon data.
!
      if (.not.lneutraldensity_nolog.and.irhon/=0) &
        f(l1:l2,m1:m2,n1:n2,irhon) = exp(f(l1:l2,m1:m2,n1:n2,ilnrhon))

    endsubroutine neutraldensity_before_boundary
!***********************************************************************
    subroutine dlnrhon_dt(f,df,p)
!
!  continuity equation
!  calculate dlnrhon/dt = - un.gradlnrhon - divun
!
!  28-feb-07/wlad: adapted
!
      use Deriv, only: der6
      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_diffrhon3
      integer :: j
!
!  identify module and boundary conditions
!
      if (headtt.or.ldebug) print*,'dlnrhon_dt: SOLVE dlnrhon_dt'
      if (headtt) then
        call identify_bcs('lnrhon',ilnrhon)
        if (.not.lneutraldensity_nolog.and.irhon/=0) call identify_bcs('rhon',irhon)
      endif
!
!  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 the helium fraction Y_He=0.079 is to be included, we can redefine alpha correspondingly.
!
      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
         if (luse_as_ionization) then
           df(l1:l2,m,n,ilnrhon) = df(l1:l2,m,n,ilnrhon) - p%alpha*p%rhon1
         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
      endif
!
!  Normal mass diffusion
!
      if (dimensionality>0) then

        fdiff=0.0
        diffus_diffrhon=0.; diffus_diffrhon3=0.
!
        if (ldiffn_normal) then  ! Normal diffusion operator
          fdiff = fdiff + diffrhon*p%del2rhon
          if (lupdate_courant_dt) 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 (lupdate_courant_dt) 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 (lupdate_courant_dt) 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
          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 (lupdate_courant_dt) 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) fdiff = fdiff + diffrhon_hyper3*p%del6lnrhon
          if (lupdate_courant_dt) 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%del2rhon + diffrhon_shock*gshockglnrhon
          endif
          if (lupdate_courant_dt) 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 (lupdate_courant_dt) 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
      endif
!
!  Apply border profile
!
      if (lborder_profiles) call set_border_neutraldensity(f,df,p)
!
      call calc_diagnostics_neutraldens(p)

    endsubroutine dlnrhon_dt
!***********************************************************************
    subroutine calc_diagnostics_neutraldens(p)
! 
      use Diagnostics

      type (pencil_case) :: 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_mn_name(p%alpha,idiag_alprec)
        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 (dimensionality>0) then
          if (idiag_dtnd/=0)   call max_mn_name(diffus_diffrhon/cdtv,idiag_dtnd,l_dt=.true.)
        endif
      endif
!
    endsubroutine calc_diagnostics_neutraldens
!***********************************************************************
    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')
         return
      endselect
!
      if (lneutraldensity_nolog) f_target=exp(f_target)
      call border_driving(f,df,p,f_target,ilnrhon)
!
    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.; idiag_alprec=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)
        call parse_name(iname,cname(iname),cform(iname),'alprec',idiag_alprec)
      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