! $Id$
!
! This module solves the Poisson equation for pressure for
! the linearized density
!
!** 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 = .false.
! CPARAM logical, parameter :: lanelastic = .true.
! CPARAM logical, parameter :: lboussinesq = .false.
!
! MVAR CONTRIBUTION 0
! MAUX CONTRIBUTION 5
! COMMUNICATED AUXILIARIES 5
!
! PENCILS PROVIDED glnrho(3); grho(3); gpp(3);
! PENCILS PROVIDED uglnrho; ugrho
! PENCILS PROVIDED glnrho2; del2lnrho; del2rho; del6lnrho; del6rho
! PENCILS PROVIDED hlnrho(3,3); sglnrho(3); uij5glnrho(3),transprho
! PENCILS PROVIDED ekin
! PENCILS PROVIDED rho; rho1; lnrho; glnrhos(3)
!
!***************************************************************
module Density
!
use Cdata
use EquationOfState, only: cs0, cs20, cs2top, cs2bot, get_gamma_etc, rho0, lnrho0, &
get_average_pressure, select_eos_variable
use General, only: keep_compiler_quiet
use Messages
use Sub
use Diagnostics
!
use Special
!
implicit none
!
include '../density.h'
!
real, dimension (ninit) :: ampllnrho=0.0, widthlnrho=0.1
real, dimension (ninit) :: rho_left=1.0, rho_right=1.0
real, dimension (ninit) :: amplrho=0.0, phase_lnrho=0.0, radius_lnrho=0.5
real, dimension (ninit) :: kx_lnrho=1.0, ky_lnrho=1.0, kz_lnrho=1.0
real, dimension (ninit) :: kxx_lnrho=0.0, kyy_lnrho=0.0, kzz_lnrho=0.0
real :: lnrho_const=0.0, rho_const=1.0
real :: cdiffrho=0.0, diffrho=0.0, diffrho_hyper3=0.0, diffrho_shock=0.0
real :: eps_planet=0.5, q_ell=5.0, hh0=0.0
real :: xblob=0., yblob=0., zblob=0.
real :: co1_ss=0.,co2_ss=0.,Sigma1=150.
real :: lnrho_int=0.,lnrho_ext=0.,damplnrho_int=0.,damplnrho_ext=0.
real :: wdamp=0.,density_floor=-1.0
real :: radial_percent_smooth=10.,rshift=0.0
real, dimension(3) :: diffrho_hyper3_aniso=0.
real, dimension(mz) :: lnrho_init_z=0.0,del2lnrho_init_z=0.0
real, dimension(mz) :: dlnrhodz_init_z=0.0, glnrho2_init_z=0.0
real, dimension(3) :: beta_glnrho_global=0.0, beta_glnrho_scaled=0.0
real, target :: plaw=0.0
real :: lnrho_z_shift=0.0
real, dimension (mz,1) :: lnrhomz
real, dimension (nz,3) :: glnrhomz
real :: powerlr=3.0, zoverh=1.5, hoverr=0.05
real :: init_average_density
real, target :: mpoly=impossible
integer, parameter :: ndiff_max=4
logical :: lcontinuity_gas=.true.
logical :: lupw_lnrho=.false.,lupw_rho=.false.
logical :: ldiff_normal=.false.,ldiff_hyper3=.false.,ldiff_shock=.false.
logical :: ldiff_hyper3lnrho=.false.,ldiff_hyper3_aniso=.false.
logical :: ldiff_hyper3_polar=.false.,lanti_shockdiffusion=.false.
logical :: lfreeze_lnrhoint=.false.,lfreeze_lnrhoext=.false.
logical :: lfreeze_lnrhosqu=.false.,lexponential_smooth=.false.
logical :: lrho_as_aux=.false., ldiffusion_nolog=.false.
logical :: lshare_plaw=.false.,lmassdiff_fix=.false.
logical :: lcheck_negative_density=.false.
logical :: lcalc_glnrhomean=.false., lcalc_lnrhomean=.false.
logical, pointer :: lanelastic_lin
!
character (len=labellen), dimension(ninit) :: initlnrho='nothing'
character (len=labellen) :: strati_type='lnrho_ss'
character (len=labellen), dimension(ndiff_max) :: idiff=''
character (len=intlen) :: iinit_str
complex :: coeflnrho=0.
!
integer :: iglobal_gg=0
integer :: niter=1
!
namelist /density_init_pars/ &
ampllnrho,initlnrho,widthlnrho, &
rho_left,rho_right,lnrho_const,rho_const,cs2bot,cs2top, &
radius_lnrho,eps_planet,xblob,yblob,zblob, &
b_ell,q_ell,hh0,rbound, &
mpoly,strati_type,beta_glnrho_global,radial_percent_smooth, &
kx_lnrho,ky_lnrho,kz_lnrho,amplrho,phase_lnrho,coeflnrho, &
kxx_lnrho, kyy_lnrho, kzz_lnrho, &
co1_ss,co2_ss,Sigma1,idiff,ldensity_nolog,lexponential_smooth,&
wdamp,plaw,lcontinuity_gas,density_floor,lanti_shockdiffusion,&
rshift,lrho_as_aux,ldiffusion_nolog,lnrho_z_shift, &
lshare_plaw, powerlr, zoverh, hoverr
!
namelist /density_run_pars/ &
cdiffrho,diffrho,diffrho_hyper3,diffrho_shock, &
cs2bot,cs2top,lupw_lnrho,lupw_rho,idiff, &
lnrho_int,lnrho_ext,damplnrho_int,damplnrho_ext, &
wdamp,lfreeze_lnrhoint,lfreeze_lnrhoext,beta_glnrho_global, &
lnrho_const,plaw,lcontinuity_gas, &
diffrho_hyper3_aniso,lfreeze_lnrhosqu,density_floor, &
lanti_shockdiffusion,lrho_as_aux,ldiffusion_nolog, &
lcheck_negative_density,lmassdiff_fix,niter
!
! diagnostic variables (need to be consistent with reset list below)
!
integer :: idiag_rhom=0 ! DIAG_DOC: $\left<\varrho\right>$
! DIAG_DOC: \quad(mean density)
integer :: idiag_rho2m=0 ! DIAG_DOC:
integer :: idiag_lnrho2m=0 ! DIAG_DOC:
integer :: idiag_drho2m=0 ! DIAG_DOC:
integer :: idiag_drhom=0 ! DIAG_DOC:
integer :: idiag_rhomin=0 ! DIAG_DOC:
integer :: idiag_rhomax=0 ! DIAG_DOC:
integer :: idiag_ugrhom=0 ! DIAG_DOC: $\left<\uv\cdot\nabla\varrho\right>$
integer :: idiag_uglnrhom=0 ! DIAG_DOC:
integer :: idiag_lnrhomphi=0 ! PHIAVG_DOC: $\left<\ln\varrho\right>_\varphi$
integer :: idiag_rhomphi=0 ! PHIAVG_DOC: $\left<\varrho\right>_\varphi$
integer :: idiag_dtd=0 ! DIAG_DOC:
integer :: idiag_rhomz=0 ! DIAG_DOC:
integer :: idiag_rhomy=0 ! DIAG_DOC:
integer :: idiag_rhomx=0 ! DIAG_DOC:
integer :: idiag_rhomxy=0 ! DIAG_DOC:
integer :: idiag_rhomxz=0 ! DIAG_DOC:
integer :: idiag_rhomr=0 ! DIAG_DOC:
integer :: idiag_totmass=0 ! DIAG_DOC:
integer :: idiag_mass=0 ! DIAG_DOC: $\int\varrho\,dV$
integer :: idiag_divrhoum=0 ! DIAG_DOC: $\left<\nabla\cdot(\varrho\uv)\right>$
integer :: idiag_divrhourms=0 ! DIAG_DOC: $\left|\nabla\cdot(\varrho\uv)\right|_{\rm rms}$
integer :: idiag_divrhoumax=0 ! DIAG_DOC: $\left|\nabla\cdot(\varrho\uv)\right|_{\rm max}$
!
real :: gamma, gamma_m1
contains
!***********************************************************************
subroutine register_density()
!
! Initialise variables which should know that we solve the
! compressible hydro equations: ilnrho; increase nvar accordingly.
!
! 4-jun-02/axel: adapted from hydro
!
use FArrayManager
!
call farray_register_auxiliary('pp',ipp,communicated=.true.)
call farray_register_auxiliary('rhs',irhs,vector=3,communicated=.true.)
! call farray_register_auxiliary('divu',idivu,communicated=.false.)
!
! Identify version number (generated automatically by CVS).
!
if (lroot) call svn_id( &
"$Id$")
!
endsubroutine register_density
!***********************************************************************
subroutine initialize_density(f)
!
! Perform any post-parameter-read initialization i.e. calculate derived
! parameters.
!
! For compatibility with other applications, we keep the possibility
! of giving diffrho units of dxmin*cs0, but cs0 is not well defined general
!
! 24-nov-02/tony: coded
! 31-aug-03/axel: normally, diffrho should be given in absolute units
!
use Deriv, only: der_pencil,der2_pencil
use FArrayManager, only: farray_register_auxiliary, farray_register_global
use Gravity, only: lnumerical_equilibrium
use Mpicomm, only: stop_it
use SharedVariables, only: get_shared_variable, put_shared_variable
use DensityMethods, only: initialize_density_methods
!
real, dimension (mx,my,mz,mfarray) :: f
integer :: i,ierr
logical :: lnothing
!
call get_shared_variable('lanelastic_lin',lanelastic_lin,caller='initialize_density')
!
if (lanelastic_lin) then
call farray_register_auxiliary('rho_b',irho_b,communicated=.true.)
else
call farray_register_auxiliary('rho',irho,communicated=.true.)
endif
!
! initialize cs2cool to cs20
! (currently disabled, because it causes problems with mdwarf auto-test)
! cs2cool=cs20
!
if (diffrho==0.) then
!
! Made to work by adding diffrho + cdiffrho to the rprint reset list.
!
diffrho=cdiffrho*dxmin*cs0
endif
!
! Turn off continuity equation.
!
lcontinuity_gas=.false.
print*, 'initialize_density: density_anelastic, turned off continuity equation'
!
! Initialize mass diffusion
!
ldiff_normal=.false.
ldiff_shock=.false.
ldiff_hyper3=.false.
ldiff_hyper3lnrho=.false.
ldiff_hyper3_aniso=.false.
ldiff_hyper3_polar=.false.
!
lnothing=.false.
!
do i=1,ndiff_max
select case (idiff(i))
case ('','none')
if (lroot .and. (.not. lnothing)) print*,'no mass diffusion'
case default
write(unit=errormsg,fmt=*) 'initialize_density: ', &
'You cannot have mass diffusion in anelastic approximation', trim(idiff(i))
call fatal_error('initialize_density',errormsg)
endselect
lnothing=.true.
enddo
!
! Tell the equation of state that we're here and what f variable we use
! DM+PC
! For anelastic case use pressure
call select_eos_variable('pp',ipp)
!
! Do not allow inconsistency between rho0 (from eos) and rho_const
! or lnrho0 and lnrho_const.
!
if (rho0/=rho_const) then
if (lroot) then
print*,"WARNING!"
print*,"inconsistency between the density constants from eos "
print*,"(rho0 or lnrho0) and the ones from the density module "
print*,"(rho_const or lnrho_const). It may damage your "
print*,"simulation if you are using them in different places. "
call warning("initialize_density","")
endif
endif
!
if (lnumerical_equilibrium) then
if (lroot) print*,'initializing global gravity in density'
call farray_register_global('gg',iglobal_gg,vector=3)
endif
call put_shared_variable('mpoly',mpoly)
call initialize_density_methods
!
call get_gamma_etc(gamma); gamma_m1=gamma-1.
endsubroutine initialize_density
!***********************************************************************
subroutine init_lnrho(f)
!
! Initialise logarithmic or non-logarithmic density.
!
! 7-nov-01/wolf: coded
! 28-jun-02/axel: added isothermal
! 15-oct-03/dave: added spherical shell (kws)
!
use General, only: itoa,complex_phase,notanumber
use Gravity, only: zref,z1,z2,gravz,nu_epicycle,potential, &
lnumerical_equilibrium
use Initcond
use IO
use Mpicomm
! use Selfgravity, only: rhs_poisson_const
use InitialCondition, only: initial_condition_lnrho
use Poisson, only: inverse_laplacian
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (nx) :: pot,prof
real, dimension (nx,ny,nz) :: psi
real, dimension (ninit) :: lnrho_left,lnrho_right
real :: lnrhoint,cs2int,pot0
real :: pot_ext,lnrho_ext,cs2_ext,tmp1,k_j2
real :: zbot,ztop,haut
real, dimension(1) :: mass_per_proc,pres_per_proc
real :: average_density
real, dimension (nx) :: r_mn,lnrho,TT,ss
real, pointer :: gravx
complex :: omega_jeans
integer :: j, ierr
logical :: lnothing
!
intent(inout) :: f
type (pencil_case) :: p
!
! 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 lnrho (called from start).
!
lnrho0 = log(rho0)
lnrho_left = log(rho_left)
lnrho_right = log(rho_right)
!
lnothing=.true.
!
do j=1,ninit
!
if (initlnrho(j)=='nothing') cycle
!
lnothing=.false.
!
iinit_str=itoa(j)
!
select case (initlnrho(j))
!
case ('const')
if (lroot) print*,'initialize anelastic: const'
do n=1,mz; do m=1,my
if (lanelastic_lin) then
f(1:mx,m,n,irho_b)=rho0
else
f(1:mx,m,n,irho)=rho0
endif
enddo; enddo
f(:,:,:,ipp)=0.
!
case ('test-poisson')
if (lroot) print*,'initialize anelastic: test-poisson'
do n=n1,n2; do m=m1,m2
f(l1:l2,m,n,ipp)=-2.*sin(x(l1:l2))*cos(z(n))
enddo; enddo
! call inverse_laplacian(f,f(l1:l2,m1:m2,n1:n2,ipp))
call inverse_laplacian_z(f,f(l1:l2,m1:m2,n1:n2,ipp))
print*, 'print results in binary file'
write(11) f(l1:l2,4,n1:n2,ipp)
!
case ('-ln(1+u2/2cs02)')
f(:,:,:,ilnrho) = -alog(1. &
+(f(:,:,:,iux)**2+f(:,:,:,iuy)**2+f(:,:,:,iuz)**2)/(2.*cs0**2))
!
case ('anelastic')
do m=1,my; do n=1,mz
if (lanelastic_lin) then
f(1:mx,m,n,ipp)=0.0
f(1:mx,m,n,irho_b)=rho0*exp(gamma*gravz*z(n)/cs20) ! Define the base state density
else
f(1:mx,m,n,irho)=rho0*exp(gamma*gravz*z(n)/cs20)
f(1:mx,m,n,ipp)=f(1:mx,m,n,irho)*cs20
endif
enddo; enddo
!
case ('polytropic_simple')
if (lanelastic_lin) then
call polytropic_simple(f)
f(:,:,:,ipp)=0.0
else
call fatal_error('init_lnrho','Not coded yet')
endif
!
case default
!
! Catch unknown values
!
write(unit=errormsg,fmt=*) 'No such value for initlnrho(' &
//trim(iinit_str)//'): ',trim(initlnrho(j))
call fatal_error('init_lnrho',errormsg)
!
endselect
!
! if the ipp f-array exists (e.g. in anelastic problems), set it
! (for now corresponding to an isothermal eos)
!
! if (ipp/=0.and.leos) f(:,:,:,ipp) = f(:,:,:,irho)*cs20
!
if (lroot) print*,'init_lnrho: initlnrho('//trim(iinit_str)//') = ', &
trim(initlnrho(j))
!
enddo ! End loop over initial conditions.
!
! Interface for user's own initial condition
!
if (linitial_condition) call initial_condition_lnrho(f)
!
if (lnothing.and.lroot) print*,'init_lnrho: nothing'
!
! check that cs2bot,cs2top are ok
! for runs with ionization or fixed ionization, don't print them
!
if (leos_ionization) then
cs2top=impossible
cs2bot=impossible
else
if (lroot) print*,'init_lnrho: cs2bot,cs2top=',cs2bot,cs2top
endif
!
! If unlogarithmic density considered, take exp of lnrho resulting from
! initlnrho
!
! if (ldensity_nolog) f(:,:,:,irho)=exp(f(:,:,:,ilnrho))
!
! sanity check
!
! if (notanumber(f(l1:l2,m1:m2,n1:n2,irho))) then
! call error('init_rho', 'Infinit density values')
! endif
!
endsubroutine init_lnrho
!**********************************************************************
subroutine density_after_boundary(f)
!
real, dimension (mx,my,mz,mfarray), intent(in) :: f
!
call keep_compiler_quiet(f)
!
endsubroutine density_after_boundary
!***********************************************************************
subroutine numerical_equilibrium(f)
!
! sets gravity gg in order to achieve an numerical exact equilbrium
! at t=0. This is only valid for the polytropic case, i.e.
!
! (1/rho) grad(P) = cs20 (rho/rho0)^(gamma-2) grad(rho)
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (nx) :: lnrho,cs2
real, dimension (nx,3) :: glnrho
real, dimension (nx,3) :: gg_mn
integer :: i,j,ilnrho
!
do m=m1,m2
do n=n1,n2
lnrho=f(l1:l2,m,n,ilnrho)
cs2=cs20*exp(gamma_m1*(lnrho-lnrho0))
call grad(f,ilnrho,glnrho)
do j=1,3
gg_mn(:,j)=cs2*glnrho(:,j)
enddo
f(l1:l2,m,n,iglobal_gg:iglobal_gg+2)=gg_mn
enddo
enddo
!
endsubroutine numerical_equilibrium
!***********************************************************************
subroutine pencil_criteria_density()
!
! All pencils that the Density module depends on are specified here.
!
! 19-11-04/anders: coded
!
lpenc_requested(i_rho)=.true.
! lpenc_requested(i_pp)=.true.
! lpenc_requested(i_lnrho)=.true.
! lpenc_requested(i_ugrho)=.true.
!
! lpenc_diagnos2d(i_lnrho)=.true.
lpenc_diagnos2d(i_rho)=.true.
!
! if (idiag_lnrho2m/=0) lpenc_diagnos(i_lnrho)=.true.
! if (idiag_ugrhom/=0) lpenc_diagnos(i_ugrho)=.true.
! if (idiag_uglnrhom/=0) lpenc_diagnos(i_uglnrho)=.true.
! if (idiag_divrhoum/=0.or.idiag_divrhourms/=0..or.idiag_divrhoumax/=0.) then
! lpenc_diagnos(i_rho)=.true.
!! lpenc_diagnos(i_uglnrho)=.true.
! lpenc_diagnos(i_ugrho)=.true.
! lpenc_diagnos(i_divu)=.true.
! endif
!
endsubroutine pencil_criteria_density
!***********************************************************************
subroutine pencil_interdep_density(lpencil_in)
!
! Interdependency among pencils from the Density module is specified here.
!
! 19-11-04/anders: coded
!
logical, dimension(npencils) :: lpencil_in
!
! if (ldensity_nolog) then
! if (lpencil_in(i_rho1)) lpencil_in(i_rho)=.true.
! else
! if (lpencil_in(i_rho)) lpencil_in(i_rho1)=.true.
! endif
! if (lpencil_in(i_uglnrho)) then
! lpencil_in(i_uu)=.true.
! lpencil_in(i_glnrho)=.true.
! endif
! if (lpencil_in(i_ugrho)) then
! lpencil_in(i_uu)=.true.
! lpencil_in(i_grho)=.true.
! endif
! if (lpencil_in(i_glnrho2)) lpencil_in(i_glnrho)=.true.
! if (lpencil_in(i_sglnrho)) then
! lpencil_in(i_sij)=.true.
! lpencil_in(i_glnrho)=.true.
! endif
! if (lpencil_in(i_uij5glnrho)) then
! lpencil_in(i_uij5)=.true.
! lpencil_in(i_glnrho)=.true.
! endif
! The pencils glnrho and grho come in a bundle.
! if (lpencil_in(i_glnrho) .and. lpencil_in(i_grho)) then
! if (ldensity_nolog) then
! lpencil_in(i_grho)=.false.
! else
! lpencil_in(i_glnrho)=.false.
! endif
! endif
!
endsubroutine pencil_interdep_density
!***********************************************************************
subroutine calc_pencils_density(f,p)
!
! Dummy routine copied from nodensity.f90
!
! 20-11-04/anders: coded
!
use EquationOfState, only: lnrho0, rho0
!
real, dimension (mx,my,mz,mfarray) :: f
type (pencil_case) :: p
!
intent(in) :: f
intent(inout) :: p
integer :: i, mm, nn, ierr,l, irhoxx
!
if (ldensity_nolog) call fatal_error('density_anelastic','working with lnrho')
!
if (lanelastic_lin) then
irhoxx=irho_b
else
irhoxx=irho
endif
!
p%rho=f(l1:l2,m,n,irhoxx)
p%rho1=1./p%rho
p%lnrho = log(p%rho)
! glnrho and grho
if (lpencil(i_grho)) call grad(f, irhoxx, p%grho)
if (lpencil(i_glnrho)) then
p%glnrho(:,1)=p%grho(:,1)*p%rho1
p%glnrho(:,2)=p%grho(:,2)*p%rho1
p%glnrho(:,3)=p%grho(:,3)*p%rho1
endif
! uglnrho
if (lpencil(i_uglnrho)) call dot(p%uu,p%glnrho,p%uglnrho)
! ugrho
if (lpencil(i_ugrho)) then
call grad(f,irhoxx,p%grho)
call dot(p%uu,p%grho,p%ugrho)
endif
! sglnrho
if (lpencil(i_sglnrho)) call multmv(p%sij,p%glnrho,p%sglnrho)
!
endsubroutine calc_pencils_density
!***********************************************************************
subroutine density_before_boundary(f)
!
! Actions to take before boundary conditions are set.
!
! 2-apr-08/anders: coded
!
real, dimension (mx,my,mz,mfarray), intent(inout) :: f
!
! if ( (.not.ldensity_nolog) .and. (irho/=0) ) &
! f(l1:l2,m1:m2,n1:n2,irho)=exp(f(l1:l2,m1:m2,n1:n2,ilnrho))
!
endsubroutine density_before_boundary
!***********************************************************************
subroutine dlnrho_dt(f,df,p)
!
! Dummy routine (taken from nodensity.f90 )
! 14-oct-09/dhruba: coded
!
use sub
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx,my,mz,mvar) :: df
type (pencil_case) :: p
!
intent(in) :: f,df,p
!!
! Identify module and boundary conditions.
!
if (headtt.or.ldebug) print*,'dlnrho_dt: not SOLVING dlnrho_dt in anelastic'
! if (headtt) call identify_bcs('lnrho',ilnrho)
if (headtt) call identify_bcs('pp',ipp)
if (headtt) call identify_bcs('rhs',irhs)
if (headtt) call identify_bcs('rhs',irhs+1)
if (headtt) call identify_bcs('rhs',irhs+2)
! if (headtt) call identify_bcs('idivu',idivu)
!
! Calculate density diagnostics
!
if (ldiagnos) then
if (idiag_rhom/=0) call sum_mn_name(p%rho,idiag_rhom)
if (idiag_totmass/=0) call sum_mn_name(p%rho,idiag_totmass,lint=.true.)
if (idiag_mass/=0) call integrate_mn_name(p%rho,idiag_mass)
if (idiag_rhomin/=0) &
call max_mn_name(-p%rho,idiag_rhomin,lneg=.true.)
if (idiag_rhomax/=0) call max_mn_name(p%rho,idiag_rhomax)
if (idiag_rho2m/=0) call sum_mn_name(p%rho**2,idiag_rho2m)
if (idiag_lnrho2m/=0) call sum_mn_name(p%lnrho**2,idiag_lnrho2m)
if (idiag_drho2m/=0) call sum_mn_name((p%rho-rho0)**2,idiag_drho2m)
if (idiag_drhom/=0) call sum_mn_name(p%rho-rho0,idiag_drhom)
if (idiag_ugrhom/=0) call sum_mn_name(p%ugrho,idiag_ugrhom)
if (idiag_divrhoum/=0) &
call sum_mn_name(p%rho*p%divu+p%ugrho,idiag_divrhoum)
if (idiag_divrhourms/=0) call sum_mn_name((p%rho*p%divu+p%ugrho)**2,idiag_divrhourms,lsqrt=.true.)
if (idiag_divrhoumax/=0) call max_mn_name(p%rho*p%divu+p%ugrho,idiag_divrhoumax)
!
! MR: diffus_diffrho is never set
!
! if (idiag_dtd/=0) &
! call max_mn_name(diffus_diffrho/cdtv,idiag_dtd,l_dt=.true.)
endif
!
endsubroutine dlnrho_dt
!***********************************************************************
subroutine split_update_density(f)
!
real, dimension(mx,my,mz,mfarray), intent(in) :: f
!
call keep_compiler_quiet(f)
!
endsubroutine split_update_density
!***********************************************************************
! Here comes a collection of different density stratification routines
!***********************************************************************
subroutine isothermal_density(f)
!
! Isothermal stratification (for lnrho and ss)
! This routine should be independent of the gravity module used.
! When entropy is present, this module also initializes entropy.
!
! Sound speed (and hence Temperature), and density (at infinity) are
! initialised to their respective reference values:
! sound speed: cs^2_0 from start.in
! density: rho0 = exp(lnrho0)
!
! 8-jul-02/axel: incorporated/adapted from init_lnrho
! 11-jul-02/axel: fixed sign; should be tmp=gamma*pot/cs20
! 02-apr-03/tony: made entropy explicit rather than using tmp/-gamma
! 11-jun-03/tony: moved entropy initialisation to separate routine
! to allow isothermal condition for arbitrary density
!
use Gravity
use SharedVariables, only: get_shared_variable
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (nx) :: pot,tmp
real, pointer :: cp
real :: cp1
!
! Stratification depends on the gravity potential
!
if (leos_idealgas) then
call get_shared_variable('cp',cp); cp1=1./cp
endif
if (lroot) print*,'isothermal_density: isothermal stratification'
if (gamma/=1.0) then
if ((.not. lentropy) .and. (.not. ltemperature)) &
call fatal_error('isothermal_density','for gamma/=1.0, you need entropy or temperature!');
endif
!
do n=n1,n2
do m=m1,m2
call potential(x(l1:l2),y(m),z(n),pot=pot)
tmp=-gamma*pot/cs20
f(l1:l2,m,n,ilnrho) = f(l1:l2,m,n,ilnrho) + lnrho0 + tmp
if (lentropy) f(l1:l2,m,n,iss) = f(l1:l2,m,n,iss) &
-gamma_m1*(f(l1:l2,m,n,ilnrho)-lnrho0)/gamma
if (ltemperature) f(l1:l2,m,n,ilnTT)=log(cs20*cp1/gamma_m1)
enddo
enddo
!
! cs2 values at top and bottom may be needed to boundary conditions.
! The values calculated here may be revised in the entropy module.
!
cs2bot=cs20
cs2top=cs20
!
endsubroutine isothermal_density
!***********************************************************************
subroutine polytropic_simple(f)
!
! Polytropic stratification (for lnrho and ss)
! This routine should be independent of the gravity module used.
!
! To maintain continuity with respect to the isothermal case,
! one may want to specify cs20 (=1), and so zinfty is calculated from that.
! On the other hand, for polytropic atmospheres it may be more
! physical to specify zinfty (=1), ie the top of the polytropic atmosphere.
! This is done if zinfty is different from 0.
!
! 8-jul-02/axel: incorporated/adapted from init_lnrho
!
use Gravity, only: gravz_profile,gravz,zinfty,zref,zgrav, &
potential,nu_epicycle
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (nx) :: pot,dlncs2,r_mn
real :: ggamma,ztop,zbot,zref2,pot_ext,lnrho_ref,ptop,pbot
!
! identifier
!
if (lroot) print*,'polytropic_simple: mpoly=',mpoly
!
! The following is specific only to cases with gravity in the z direction
! zref is calculated such that rho=rho0 and cs2=cs20 at z=zref.
! Note: gravz is normally negative!
!
if (lgravz) then
if (gravz_profile=='const') then
if (lroot.and.gravz==0.) print*,'polytropic_simple: divide by gravz=0'
zref=zinfty-(mpoly+1.)*cs20/(-gamma*gravz)
elseif (gravz_profile=='const_zero') then
if (lroot.and.gravz==0.) print*,'polytropic_simple: divide by gravz=0'
zref=zinfty-(mpoly+1.)*cs20/(-gamma*gravz)
elseif (gravz_profile=='linear') then
if (lroot.and.gravz==0.) print*,'polytropic_simple: divide by gravz=0'
zref2=zinfty**2-(mpoly+1.)*cs20/(0.5*gamma*nu_epicycle**2)
if (zref2<0) then
if (lroot) print*,'polytropic_simple: zref**2<0 is not ok'
zref2=0. !(and see what happens)
endif
zref=sqrt(zref2)
else
if (lroot) print*,'polytropic_simple: zref not prepared!'
endif
if (lroot) print*,'polytropic_simple: zref=',zref
!
! check whether zinfty lies outside the domain (otherwise density
! would vanish within the domain). At the moment we are not properly
! testing the lower boundary on the case of a disc (commented out).
!
ztop=xyz0(3)+Lxyz(3)
zbot=xyz0(3)
!-- if (zinfty<min(ztop,zgrav) .or. (-zinfty)>min(zbot,zgrav)) then
if (zinfty<min(ztop,zgrav)) then
if (lroot) print*,'polytropic_simple: domain too big; zinfty=',zinfty
!call stop_it( &
! 'polytropic_simply: rho and cs2 will vanish within domain')
endif
endif
!
! stratification Gamma (upper case in the manual)
!
ggamma=1.+1./mpoly
!
! polytropic sphere with isothermal exterior
! calculate potential at the stellar surface, pot_ext
!
if (lgravr) then
call potential(R=r_ext,POT=pot_ext)
cs2top=-gamma/(mpoly+1.)*pot_ext
lnrho_ref=mpoly*log(cs2top)-(mpoly+1.)
print*,'polytropic_simple: pot_ext=',pot_ext
do n=n1,n2; do m=m1,m2
r_mn=sqrt(x(l1:l2)**2+y(m)**2+z(n)**2)
call potential(x(l1:l2),y(m),z(n),pot=pot)
!
! density
! these formulae assume lnrho0=0 and cs0=1
!
where (r_mn > r_ext)
f(l1:l2,m,n,irho_b)=exp(lnrho_ref-gamma*pot/cs2top)
elsewhere
dlncs2=log(-gamma*pot/((mpoly+1.)*cs20))
f(l1:l2,m,n,irho_b)=exp(lnrho0+mpoly*dlncs2)
endwhere
!
! entropy
!
if (lentropy) then
where (r_mn > r_ext)
f(l1:l2,m,n,iss_b)=-(1.-1./gamma)*log(f(l1:l2,m,n,irho_b))+log(cs2top)/gamma
elsewhere
dlncs2=log(-gamma*pot/((mpoly+1.)*cs20))
f(l1:l2,m,n,iss_b)=mpoly*(ggamma/gamma-1.)*dlncs2
endwhere
endif
enddo; enddo
else
!
! cartesian case with gravity in the z direction
!
do n=n1,n2; do m=m1,m2
call potential(x(l1:l2),y(m),z(n),pot=pot)
dlncs2=log(-gamma*pot/((mpoly+1.)*cs20))
f(l1:l2,m,n,irho_b)=rho0*exp(mpoly*dlncs2)
! if (lentropy) f(l1:l2,m,n,iss_b)=mpoly*(ggamma/gamma-1.)*dlncs2
if (lentropy) f(l1:l2,m,n,iss_b)=(ggamma/gamma-1.)*log(f(l1:l2,m,n,irho_b)/rho0)
enddo; enddo
do n=n1,n2
write(*,*) f(l1,m1,n,irho_b)
end do
!
! cs2 values at top and bottom may be needed to boundary conditions.
! In spherical geometry, ztop is z at the outer edge of the box,
! so this calculation still makes sense.
!
call potential(xyz0(1),xyz0(2),ztop,pot=ptop)
cs2top=-gamma*ptop/(mpoly+1.)
!
! In spherical geometry ztop should never be used.
! Even in slab geometry ztop is not normally used.
!
call potential(xyz0(1),xyz0(2),zbot,pot=pbot)
cs2bot=-gamma*pbot/(mpoly+1.)
endif
!
endsubroutine polytropic_simple
!***********************************************************************
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 compressible part
!
! 3-may-02/axel: coded
! 27-may-02/axel: added possibility to reset list
!
use Diagnostics
use FArrayManager, only: farray_index_append
!
logical :: lreset
logical, optional :: lwrite
!
integer :: iname, inamex, inamey, inamez, inamexy, inamexz, 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_rhom=0; idiag_rho2m=0; idiag_lnrho2m=0
idiag_drho2m=0; idiag_drhom=0
idiag_ugrhom=0; idiag_uglnrhom=0
idiag_rhomin=0; idiag_rhomax=0; idiag_dtd=0
idiag_lnrhomphi=0; idiag_rhomphi=0
idiag_rhomz=0; idiag_rhomy=0; idiag_rhomx=0
idiag_rhomxy=0; idiag_rhomr=0; idiag_totmass=0
idiag_rhomxz=0; idiag_divrhoum=0; idiag_divrhourms=0; idiag_divrhoumax=0
cformv=''
endif
!
! iname runs through all possible names that may be listed in print.in
!
if (lroot.and.ip<14) print*,'rprint_density: run through parse list'
do iname=1,nname
call parse_name(iname,cname(iname),cform(iname),'rhom',idiag_rhom)
call parse_name(iname,cname(iname),cform(iname),'rho2m',idiag_rho2m)
call parse_name(iname,cname(iname),cform(iname),'drho2m',idiag_drho2m)
call parse_name(iname,cname(iname),cform(iname),'drhom',idiag_drhom)
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),'lnrho2m',idiag_lnrho2m)
call parse_name(iname,cname(iname),cform(iname),'ugrhom',idiag_ugrhom)
call parse_name(iname,cname(iname),cform(iname),'uglnrhom',idiag_uglnrhom)
call parse_name(iname,cname(iname),cform(iname),'dtd',idiag_dtd)
call parse_name(iname,cname(iname),cform(iname),'totmass',idiag_totmass)
call parse_name(iname,cname(iname),cform(iname),'mass',idiag_mass)
call parse_name(iname,cname(iname),cform(iname),'divrhoum',idiag_divrhoum)
call parse_name(iname,cname(iname),cform(iname),'divrhourms',idiag_divrhourms)
call parse_name(iname,cname(iname),cform(iname),'divrhoumax',idiag_divrhoumax)
!
! alternatively, use these shorter names: drurms and drumax,
! instead of divrhourms and divrhoumax
!
call parse_name(iname,cname(iname),cform(iname),'drurms',idiag_divrhourms)
call parse_name(iname,cname(iname),cform(iname),'drumax',idiag_divrhoumax)
enddo
!
! check for those quantities for which we want xy-averages
!
do inamez=1,nnamez
call parse_name(inamez,cnamez(inamez),cformz(inamez),'rhomz',idiag_rhomz)
enddo
!
! check for those quantities for which we want xz-averages
!
do inamey=1,nnamey
call parse_name(inamey,cnamey(inamey),cformy(inamey),'rhomy',idiag_rhomy)
enddo
!
! check for those quantities for which we want yz-averages
!
do inamex=1,nnamex
call parse_name(inamex,cnamex(inamex),cformx(inamex),'rhomx',idiag_rhomx)
enddo
!
! check for those quantities for which we want phiz-averages
!
do inamer=1,nnamer
call parse_name(inamer,cnamer(inamer),cformr(inamer),'rhomr',idiag_rhomr)
enddo
!
! check for those quantities for which we want z-averages
!
do inamexz=1,nnamexz
call parse_name(inamexz,cnamexz(inamexz),cformxz(inamexz),'rhomxz',idiag_rhomxz)
enddo
!
! check for those quantities for which we want z-averages
!
do inamexy=1,nnamexy
call parse_name(inamexy,cnamexy(inamexy),cformxy(inamexy),'rhomxy',idiag_rhomxy)
enddo
!
! check for those quantities for which we want phi-averages
!
do irz=1,nnamerz
call parse_name(irz,cnamerz(irz),cformrz(irz),&
'lnrhomphi',idiag_lnrhomphi)
call parse_name(irz,cnamerz(irz),cformrz(irz),'rhomphi',idiag_rhomphi)
enddo
!
! check for those quantities for which we want video slices
!
if (lwrite_slices) then
where(cnamev=='rho'.or.cnamev=='pp') cformv='DEFINED'
endif
!
! write column where which density variable is stored
!
if (lwr) then
call farray_index_append('i_rhom',idiag_rhom)
call farray_index_append('i_rho2m',idiag_rho2m)
call farray_index_append('i_drho2m',idiag_drho2m)
call farray_index_append('i_drhom',idiag_drhom)
call farray_index_append('i_rhomin',idiag_rhomin)
call farray_index_append('i_rhomax',idiag_rhomax)
call farray_index_append('i_lnrho2m',idiag_lnrho2m)
call farray_index_append('i_ugrhom',idiag_ugrhom)
call farray_index_append('i_uglnrhom',idiag_uglnrhom)
call farray_index_append('i_rhomz',idiag_rhomz)
call farray_index_append('i_rhomy',idiag_rhomy)
call farray_index_append('i_rhomx',idiag_rhomx)
call farray_index_append('i_rhomxy',idiag_rhomxy)
call farray_index_append('i_rhomxz',idiag_rhomxz)
call farray_index_append('nname',nname)
call farray_index_append('ilnrho',ilnrho)
call farray_index_append('irho',irho)
call farray_index_append('i_lnrhomphi',idiag_lnrhomphi)
call farray_index_append('i_rhomphi',idiag_rhomphi)
call farray_index_append('i_rhomr',idiag_rhomr)
call farray_index_append('i_dtd',idiag_dtd)
call farray_index_append('i_totmass',idiag_totmass)
call farray_index_append('i_mass',idiag_mass)
call farray_index_append('i_divrhoum',idiag_divrhoum)
call farray_index_append('i_divrhourms',idiag_divrhourms)
call farray_index_append('i_divrhoumax',idiag_divrhoumax)
endif
!
endsubroutine rprint_density
!***********************************************************************
subroutine get_init_average_density(f,init_average_density)
!
! 10-dec-09/piyali: added to pass initial average density
! equ.f90
!
use Diagnostics, only: integrate_mn,get_average_density
!
real, dimension (mx,my,mz,mfarray):: f
real, dimension(1) :: mass_per_proc
real :: init_average_density
intent(in):: f
!
do m=m1,m2
do n=n1,n2
call integrate_mn(exp(f(l1:l2,m,n,ilnrho)),mass_per_proc(1))
enddo
enddo
call get_average_density(mass_per_proc(1),init_average_density)
!
endsubroutine get_init_average_density
!***********************************************************************
subroutine get_slices_density(f,slices)
!
! Write slices for animation of Density variables.
!
! 26-jul-06/tony: coded
!
use Slices_methods, only: assign_slices_scal
!
real, dimension (mx,my,mz,mfarray) :: f
type (slice_data) :: slices
!
if (.not. lanelastic_lin) then
!
! Loop over slices
!
select case (trim(slices%name))
! Density.
case ('rho'); call assign_slices_scal(slices,f,irho)
!
! Logarithmic density (not implemented).
!
!case ('lnrho')
endselect
endif
!
endsubroutine get_slices_density
!***********************************************************************
subroutine get_slices_pressure(f,slices)
!
! Write slices for animation of Pressure variables.
!
! 26-jul-06/tony: coded
!
use Slices_methods, only: assign_slices_scal
!
real, dimension (mx,my,mz,mfarray) :: f
type (slice_data) :: slices
!
! Loop over slices
!
select case (trim(slices%name))
!
! Pressure.
!
case ('pp'); call assign_slices_scal(slices,f,ipp)
!
endselect
!
endsubroutine get_slices_pressure
!***********************************************************************
subroutine density_after_mn(f, df, mass_per_proc)
!
use Poisson, only: inverse_laplacian
use Mpicomm, only: initiate_isendrcv_bdry, finalize_isendrcv_bdry
use Boundcond, only: boundconds
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx,my,mz,mvar) :: df
! real, dimension (nx,ny,nz) :: pold
real, dimension (nx,3) :: gpp
real, dimension (nx) :: phi_rhs_pencil
real, dimension (1) :: mass_per_proc
real :: average_density,average_pressure,init_average_density
integer :: j, ju, l, i
!
! Set first the boundary conditions on rhs
!
call initiate_isendrcv_bdry(f,irhs,irhs+2)
call finalize_isendrcv_bdry(f,irhs,irhs+2)
call boundconds(f,irhs,irhs+2)
!
! Find the divergence of rhs
!
! pold(1:nx,1:ny,1:nz)=f(l1:l2,m1:m2,n1:n2,ipp)
do n=n1,n2
do m=m1,m2
call div(f,irhs,phi_rhs_pencil)
f(l1:l2,m,n,ipp)=phi_rhs_pencil
enddo
enddo
!
! get pressure from inverting the Laplacian
!
if (lperi(3)) then
call inverse_laplacian(f(l1:l2,m1:m2,n1:n2,ipp))
! refresh the ghost zones: periodic pressure
call initiate_isendrcv_bdry(f,ipp)
call finalize_isendrcv_bdry(f,ipp)
call boundconds(f,ipp,ipp)
if (.not. lanelastic_lin) then
call get_average_density(mass_per_proc(1),average_density)
call get_average_pressure(init_average_density,average_density,average_pressure)
f(:,:,:,ipp)=f(:,:,:,ipp)+average_pressure
endif
else
call inverse_laplacian_z(f,f(l1:l2,m1:m2,n1:n2,ipp))
!
! fill the ghost zones:
! In the vertical direction: dP/dz=0 (symmetric)
!
m=4; n=n1
do i=1,nghost
f(l1:l2,4,n1-i,ipp)= f(l1:l2,4,n1+i,ipp)
enddo
n=n2
do i=1,nghost
f(l1:l2,4,n2+i,ipp)=f(l1:l2,4,n2-i,ipp)
enddo
! Bc in the horizontal direction: periodic
f(1:l1-1,:,:,ipp) = f(l2i:l2,:,:,ipp)
f(l2+1:mx,:,:,ipp) = f(l1:l1i,:,:,ipp)
endif
!
! Add the pressure gradient term to the NS equation
!
do n=n1,n2; do m=m1,m2
call grad(f,ipp,gpp)
do j=1,3
ju=j+iuu-1
if (lanelastic_lin) then
! df(l1:l2,m,n,ju)=df(l1:l2,m,n,ju)-gpp(:,j)/f(l1:l2,m,n,irho_b) &
! + gamma*f(l1:l2,m,n,ipp)/(f(l1:l2,m,n,irho_b)*p%cs2)*p%gg(:,j)
df(l1:l2,m,n,ju)=df(l1:l2,m,n,ju)-gpp(:,j)/f(l1:l2,m,n,irho_b)
else
df(l1:l2,m,n,ju)=df(l1:l2,m,n,ju)-gpp(:,j)/f(l1:l2,m,n,irho)
endif
enddo
!
if (.not. lanelastic_lin) then
f(l1:l2,m,n,irho)=f(l1:l2,m,n,ipp)*cs20
! call initiate_isendrcv_bdry(f,ilnrho)
! call finalize_isendrcv_bdry(f,ilnrho)
! call boundconds(f,ilnrho,ilnrho)
endif
enddo; enddo
!
endsubroutine density_after_mn
!***********************************************************************
subroutine inverse_laplacian_z(f,phi)
!
! Solve the Poisson equation by Fourier transforming in the xy-plane and
! solving the discrete matrix equation in the z-direction.
!
! 02-mar-2012/dintrans+dubuffet: coded
!
use Fourier, only: fourier_transform_xy, kx_fft, ky_fft
use General, only: tridag
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (nx,ny,nz) :: phi, b1, d2uz, b2
real, dimension (nx) :: tmpx
real, dimension (nz) :: a_tri, b_tri, c_tri, r_tri, u_tri
real :: k2
integer :: ikx, iky, j, k
logical :: err
!
! Identify version.
!
if (lroot .and. ip<10) call svn_id( &
'$Id$')
!
! test with the pressure
!
! do n=n1,n2
! do m=m1,m2
!! call der2(f,iuz,tmpx,3)
! call der(f,iuz,tmpx,3)
! d2uz(:,m-nghost,n-nghost)=1.e-3*tmpx
! enddo
! enddo
!
! The right-hand-side of the Poisson equation is purely real.
!
b1 = 0.0
b2 = 0.0
!
! Forward transform (to k-space).
!
call fourier_transform_xy(phi,b1)
! call fourier_transform_xy(d2uz,b2)
!
! Solve for discrete z-direction with zero density above and below z-boundary.
!
do iky=1,ny
do ikx=1,nx
if ((kx_fft(ikx)==0.0) .and. (ky_fft(iky)==0.0)) then
phi(ikx,iky,:) = 0.0
else
k2=kx_fft(ikx)**2+ky_fft(iky)**2
a_tri=1.0/dz**2
b_tri=-2.0/dz**2-k2
c_tri=1.0/dz**2
r_tri=phi(ikx,iky,:)
! P = 0 (useful to test the Poisson solver)
! b_tri(1)=1. ; c_tri(1)=0. ; r_tri(1)=0.
! b_tri(nz)=1. ; a_tri(nz)=0. ; r_tri(nz)=0.
! dP/dz = 0
c_tri(1)=c_tri(1)+a_tri(1)
a_tri(nz)=a_tri(nz)+c_tri(nz)
! delta P = 0 outside
! b_tri(1)=b_tri(1)+2./dz**2/(2.+2*sqrt(k2)*dz+k2*dz**2)
! b_tri(nz)=b_tri(nz)+2./dz**2/(2.+2*sqrt(k2)*dz+k2*dz**2)
! dP/dz = mu d2 uz /dz2
! c_tri(1)=c_tri(1)+a_tri(1)
! r_tri(1)=r_tri(1)+2./dz*d2uz(ikx,iky,1)
! a_tri(nz)=a_tri(nz)+c_tri(nz)
! r_tri(nz)=r_tri(nz)-2./dz*d2uz(ikx,iky,nz)
! P = mu d uz /dz
! b_tri(1)=1.
! c_tri(1)=0.
! r_tri(1)=d2uz(ikx,iky,1)
! b_tri(nz)=1.
! a_tri(nz)=0.
! r_tri(nz)=d2uz(ikx,iky,nz)
!
call tridag(a_tri,b_tri,c_tri,r_tri,u_tri,err)
phi(ikx,iky,:)=u_tri
endif
enddo
enddo
!
do iky=1,ny
do ikx=1,nx
if ((kx_fft(ikx)==0.0) .and. (ky_fft(iky)==0.0)) then
b1(ikx,iky,:) = 0.0
else
k2=kx_fft(ikx)**2+ky_fft(iky)**2
a_tri=1.0/dz**2
b_tri=-2.0/dz**2-k2
c_tri=1.0/dz**2
r_tri=b1(ikx,iky,:)
! P = 0 (useful to test the Poisson solver)
! b_tri(1)=1. ; c_tri(1)=0. ; r_tri(1)=0.
! b_tri(nz)=1. ; a_tri(nz)=0. ; r_tri(nz)=0.
! dP/dz = 0
c_tri(1)=c_tri(1)+a_tri(1)
a_tri(nz)=a_tri(nz)+c_tri(nz)
! delta P = 0 outside
! b_tri(1)=b_tri(1)+2./dz**2/(2.+2*sqrt(k2)*dz+k2*dz**2)
! b_tri(nz)=b_tri(nz)+2./dz**2/(2.+2*sqrt(k2)*dz+k2*dz**2)
! dP/dz = mu d2 uz /dz2
! c_tri(1)=c_tri(1)+a_tri(1)
! r_tri(1)=r_tri(1)+2./dz*b2(ikx,iky,1)
! a_tri(nz)=a_tri(nz)+c_tri(nz)
! r_tri(nz)=r_tri(nz)-2./dz*b2(ikx,iky,nz)
! P = mu d uz /dz
! b_tri(1)=1.
! c_tri(1)=0.
! r_tri(1)=b2(ikx,iky,1)
! b_tri(nz)=1.
! a_tri(nz)=0.
! r_tri(nz)=b2(ikx,iky,nz)
!
call tridag(a_tri,b_tri,c_tri,r_tri,u_tri,err)
b1(ikx,iky,:)=u_tri
endif
enddo
enddo
!
! Inverse transform (to real space).
!
call fourier_transform_xy(phi,b1,linv=.true.)
!
endsubroutine inverse_laplacian_z
!***********************************************************************
subroutine dynamical_diffusion(uc)
!
! Dynamically set mass diffusion coefficient given fixed mesh Reynolds number.
!
! 27-jul-11/nils: dummy
!
! Input Argument
! uc
! Characteristic velocity of the system.
!
real, intent(in) :: uc
!
call keep_compiler_quiet(uc)
call fatal_error('dynamical_diffusion',&
'This subroutine is not yet implemented for the anelastic module.')
!
endsubroutine dynamical_diffusion
!***********************************************************************
subroutine impose_density_floor(f)
!
real, dimension (mx,my,mz,mfarray), intent(in) :: f
!
call keep_compiler_quiet(f)
!
endsubroutine impose_density_floor
!***********************************************************************
subroutine boussinesq(f)
!
! 23-mar-2012/dintrans: coded
! dummy routine for the Boussinesq approximation
!
real, dimension (mx,my,mz,mfarray) :: f
!
call keep_compiler_quiet(f)
!
endsubroutine boussinesq
!***********************************************************************
function mean_density(f)
!
! Calculate mean density from f-array. With lreference_state=T it is the mean
! density deviation.
!
! 3-mar-14/MR: coded
!
use Mpicomm, only: mpiallreduce_sum
!
real :: mean_density
!
real, dimension (mx,my,mz,mfarray) :: f
intent(in) :: f
!
integer :: n,m,irhoxx
real :: tmp
!
if (lanelastic_lin) then
irhoxx=irho_b
else
irhoxx=irho
endif
!
mean_density=0.
!
do n=n1,n2
tmp=0.
do m=m1,m2
tmp=tmp+sum(f(l1:l2,m,n,irhoxx)*dVol_x(l1:l2))*dVol_y(m)
enddo
mean_density=mean_density+tmp*dVol_z(n)
enddo
!
mean_density = mean_density/box_volume
!
if (ncpus>1) then
call mpiallreduce_sum(mean_density,tmp)
mean_density=tmp
endif
!
endfunction mean_density
!***********************************************************************
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 write_z_stratification(f)
real, dimension (mx,my,mz,mfarray) :: f
call keep_compiler_quiet(f)
endsubroutine write_z_stratification
!***********************************************************************
subroutine update_char_vel_density(f)
real, dimension (mx,my,mz,mfarray) :: f
call keep_compiler_quiet(f)
endsubroutine update_char_vel_density
!***********************************************************************
subroutine impose_density_ceiling(f)
!
real, dimension (mx,my,mz,mfarray) :: f
call keep_compiler_quiet(f)
!
endsubroutine impose_density_ceiling
!***********************************************************************
subroutine pushpars2c(p_par)
use Syscalls, only: copy_addr
integer, parameter :: n_pars=0
integer(KIND=ikind8), dimension(n_pars) :: p_par
endsubroutine pushpars2c
!***********************************************************************s
endmodule Density