! $Id$
!
! Module for calculating Cosmic Ray Flux.
!
!** 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 :: lcosmicrayflux = .true.
!
! MVAR CONTRIBUTION 3
! MAUX CONTRIBUTION 0
!
!***************************************************************
module Cosmicrayflux
use Cparam
use Cdata
use General, only: keep_compiler_quiet
use Messages
implicit none
include 'cosmicrayflux.h'
character(len=labellen) :: initfcr='zero'
real, pointer :: K_para, K_perp
real :: kpara_t=0., kperp_t=0.
real :: tau=0., tau1=0.
real :: subgrid_c1=0., subgrid_c2=0.
real :: subgrid_brms=1.
real :: subgrid_bmin=0.
real :: subgrid_s=5./3.
real :: subgrid_k=0.
real :: ratio_kpara_kperp=0.
real, dimension(nx) :: vKperp, vKpara
real, dimension(nx) :: b_exp
logical :: lsubgrid_cr = .false.
logical :: lcosmicrayflux_diffus_dt = .false.
logical :: ladvect_fcr=.false., lupw_fcr=.false.
namelist /cosmicrayflux_init_pars/ &
tau, lsubgrid_cr, &
lcosmicrayflux_diffus_dt, &
subgrid_brms, subgrid_bmin, &
subgrid_c1,subgrid_c2, &
subgrid_s, subgrid_k, &
ratio_kpara_kperp
namelist /cosmicrayflux_run_pars/ &
tau, lsubgrid_cr, &
lcosmicrayflux_diffus_dt, ladvect_fcr, lupw_fcr, &
subgrid_brms, subgrid_bmin, &
subgrid_c1,subgrid_c2, &
subgrid_s, subgrid_k, &
ratio_kpara_kperp
contains
!*******************************************************************************
subroutine register_cosmicrayflux()
! Initialise variables which should know that we solve for the vector
! potential: ifcr, etc; increase nvar accordingly
!
! 1-may-02/wolf: coded
use FArrayManager
call farray_register_pde('fcr',ifcr,vector=3)
ifcrx = ifcr
ifcry = ifcr+1
ifcrz = ifcr+2
! Identify version number.
if (lroot) call svn_id("$Id$")
! Writing files for use with IDL
if (lroot) then
if (maux == 0) then
if (nvar < mvar) write (4,*) ',fcr $'
if (nvar == mvar) write (4,*) ',fcr'
else
write (4,*) ',fcr $'
endif
write (15,*) 'fcr = fltarr(mx,my,mz,3)*one'
endif
!
endsubroutine register_cosmicrayflux
!*******************************************************************************
subroutine initialize_cosmicrayflux(f)
! Perform any post-parameter-read initialization
use Messages, only: warning
use SharedVariables, only: get_shared_variable
real, dimension(mx,my,mz,mfarray) :: f
if (tau /= 0.) then
tau1 = 1./tau
else
call warning('initialize_cosmicrayflux', &
'parameter tau was set to zero!')
endif
! Reads shared diffusivities from the cosmicray module
call get_shared_variable('K_para',K_para)
call get_shared_variable('K_perp',K_perp)
kpara_t = K_para * tau1
kperp_t = K_perp * tau1
call keep_compiler_quiet(f)
endsubroutine initialize_cosmicrayflux
!*******************************************************************************
subroutine init_fcr(f)
! initialise magnetic field; called from start.f90
! AB: maybe we should here call different routines (such as rings)
! AB: and others, instead of accummulating all this in a huge routine.
! We have an init parameter (initaa) to stear magnetic i.c. independently.
!
! 7-nov-2001/wolf: coded
use Mpicomm
use Sub
use Initcond
use InitialCondition, only: initial_condition_fcr
real, dimension(mx,my,mz,mfarray) :: f
select case (initfcr)
case ('zero', '0')
f(:,:,:,ifcrx:ifcrz) = 0.
! probably no more cases needed for fcr
case default
! Catch unknown values
if (lroot) print*, 'init_fcr: No such such value for initfcr: ', trim(initfcr)
call stop_it(" ")
endselect
! Interface for user's own initial condition
if (linitial_condition) call initial_condition_fcr(f)
endsubroutine init_fcr
!*******************************************************************************
subroutine pencil_criteria_cosmicrayflux()
! All pencils that the Magnetic module depends on are specified here.
!
! 19-11-04/anders: coded
lpenc_requested(i_gecr) = .true.
lpenc_requested(i_bb) = .true.
endsubroutine pencil_criteria_cosmicrayflux
!*******************************************************************************
subroutine pencil_interdep_cosmicrayflux(lpencil_in)
logical, dimension(npencils) :: lpencil_in
call keep_compiler_quiet(lpencil_in)
endsubroutine pencil_interdep_cosmicrayflux
!*******************************************************************************
subroutine calc_pencils_cosmicrayflux(f,p)
! Calculate Cosmicray Flux pencils - to be done
real, dimension(mx,my,mz,mfarray) :: f
type (pencil_case) :: p
intent(in) :: f
intent(inout) :: p
call keep_compiler_quiet(f)
call keep_compiler_quiet(p)
! fcr
!if (lpencil(i_fcr)) p%fcr=f(l1:l2,m,n,ifcrx:ifcrz)
endsubroutine calc_pencils_cosmicrayflux
!*******************************************************************************
subroutine dfcr_dt(f,df,p)
! Cosmicray Flux evolution
! 08-mar-05/snod: adapted from daa_dt
! 12-jan-15/luiz: diffusion/advection contribution to the timestep
! 22-jun-17/luiz: subgrid prescription for the cosmic ray diffusion
use Sub
use Debug_IO, only: output_pencil
use Mpicomm, only: stop_it
use General, only: notanumber
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(mx,my,mz,mvar) :: df
real, dimension(nx,3) :: BuiBujgecr, bunit
real, dimension(nx) :: b2, b21, b_brms
real, dimension(nx) :: tmp, diffus_cr,advec_kfcr
real, dimension (nx,3,3) :: gfcr
real, dimension (nx,3) :: ugfcr
integer :: i, j
type (pencil_case) :: p
intent(in) :: f
intent(inout) :: df
! Identifies module and boundary conditions.
if (headtt .or. ldebug) print*,'dfcr_dt: SOLVE'
if (headtt) then
call identify_bcs('Fecx',ifcrx)
call identify_bcs('Fecy',ifcry)
call identify_bcs('Fecz',ifcrz)
endif
call dot2_mn(p%bb,b2)
! with frequency omega_Bz_ext
b21 = 1./max(tini,b2)
call multsv_mn(sqrt(b21),p%bb,bunit)
do i=1,3
tmp = 0.
do j=1, 3
tmp = tmp + bunit(:,i)*bunit(:,j)*p%gecr(:,j)
enddo
BuiBujgecr(:,i) = tmp
enddo
! Cosmic Ray Flux equation.
if (.not.lsubgrid_cr) then
! If not using the subgrid prescription, the diffusivities are
! constants
df(l1:l2,m,n,ifcrx:ifcrz) = df(l1:l2,m,n,ifcrx:ifcrz) &
- tau1*f(l1:l2,m,n,ifcrx:ifcrz) &
- kperp_t*p%gecr &
- (kpara_t - kperp_t)*BuiBujgecr
else
! Subgrid prescription for parallel diffusion:
! kpara = kpara_0 * [1 + c1*(B/Brms)^s + c2*(Brms/B)]
! Brms (subgrid_brms) is the rms value of the assumed background field
! and s is its spectral index.
! (subgrid_)c1 and c2 depend on properties of background field
vKpara(:) = kpara_t
b_brms = sqrt(b2)/subgrid_brms
if (subgrid_c1 /= 0.0) then
vKpara(:) = vKpara(:) + kpara_t*subgrid_c1*b_brms**subgrid_s
endif
if (subgrid_c2 /= 0.0) then
vKpara(:) = vKpara(:) + kpara_t*subgrid_c2*subgrid_brms*sqrt(b21)
endif
! Subgrid prescription for perpendicular diffusion:
if (ratio_kpara_kperp /= 0.0) then
! If a ration k_para/k_perp was specified, ignore the value of
! k_perp and computes the perpendicular diffusivity accordingly
! kperp = kpara/r * (B/Brms)^(-k)
! where r = ratio_kpara_kperp) and k = subgrid_k
vKperp(:) = vKpara(:)/ratio_kpara_kperp
if (subgrid_k /= 0.0) then
vKperp(:) = vKperp(:) * (b_brms)**(-subgrid_k)
endif
else
! Otherwise, use constant perpendicular diffusivity previously
! specified
vKperp(:) = kperp_t
endif
do i=1,3
df(l1:l2,m,n,ifcrx+i-1) = df(l1:l2,m,n,ifcrx+i-1) &
- tau1*f(l1:l2,m,n,ifcrx+i-1) &
- vKperp*p%gecr(:,i) &
- (vKpara - vKperp)*BuiBujgecr(:,i)
enddo
endif
! Allows optional use of advection term for fcr.
if (ladvect_fcr) then
call gij(f,ifcr,gfcr,1)
call u_dot_grad(f,ifcr,gfcr,p%uu,ugfcr,UPWIND=lupw_fcr)
df(l1:l2,m,n,ifcrx:ifcrz) = df(l1:l2,m,n,ifcrx:ifcrz) &
- ugfcr(1:nx,1:3)
endif
! For the timestep calculation, needs maximum diffusion or advection.
! Unless the switch lcosmicrayflux_diffus_dt is used, kpara_t and kperp_t
! are treated as an advection contribution. Otherwise,
! tau*kperp_t/tau*kpara_t are used as cosmic ray diffusivities.
if (lfirst .and. ldt) then
if (lcosmicrayflux_diffus_dt) then
if (lsubgrid_cr) then
diffus_cr = max(maxval(vKperp)*tau*dxyz_2, &
maxval(vKpara)*tau*dxyz_2)
else
diffus_cr = max(kperp_t*tau*dxyz_2,kpara_t*tau*dxyz_2)
endif
maxdiffus=max(maxdiffus,diffus_cr)
else
if (lsubgrid_cr) then
advec_kfcr = max(maxval(vKperp)*dxyz_2, &
maxval(vKpara)*dxyz_2)
else
advec_kfcr = max(kperp_t*dxyz_2, kpara_t*dxyz_2)
endif
if (notanumber(advec_kfcr)) print*, 'advec_kfcr =',advec_kfcr
advec2=advec2+advec_kfcr
endif
endif
endsubroutine dfcr_dt
!*******************************************************************************
subroutine read_cosmicrayflux_init_pars(iostat)
use File_io, only: parallel_unit
integer, intent(out) :: iostat
read(parallel_unit, NML=cosmicrayflux_init_pars, IOSTAT=iostat)
endsubroutine read_cosmicrayflux_init_pars
!*******************************************************************************
subroutine write_cosmicrayflux_init_pars(unit)
integer, intent(in) :: unit
write(unit, NML=cosmicrayflux_init_pars)
endsubroutine write_cosmicrayflux_init_pars
!*******************************************************************************
subroutine read_cosmicrayflux_run_pars(iostat)
use File_io, only: parallel_unit
integer, intent(out) :: iostat
read(parallel_unit, NML=cosmicrayflux_run_pars, IOSTAT=iostat)
endsubroutine read_cosmicrayflux_run_pars
!*******************************************************************************
subroutine write_cosmicrayflux_run_pars(unit)
integer, intent(in) :: unit
write(unit, NML=cosmicrayflux_run_pars)
endsubroutine write_cosmicrayflux_run_pars
!*******************************************************************************
subroutine rprint_cosmicrayflux(lreset,lwrite)
! Reads and registers print parameters relevant for cosmicrayflux.
! 3-may-02/axel: coded
! 27-may-02/axel: added possibility to reset list
use Sub
integer :: iname, inamez, ixy, irz
logical :: lreset
logical, optional :: lwrite
! Reset everything in case of RELOAD.
! (this needs to be consistent with what is defined above!)
if (lreset) then
!idiag_b2m=0; idiag_bm2=0; idiag_j2m=0; idiag_jm2=0; idiag_abm=0
endif
! Check for those quantities that we want to evaluate online.
do iname=1,nname
! call parse_name(iname,cname(iname),cform(iname),'dteta',idiag_dteta)
enddo
! Check for those quantities for which we want xy-averages.
do inamez=1,nnamez
! call parse_name(inamez,cnamez(inamez),cformz(inamez),'bxmz',idiag_bxmz)
enddo
! Check for those quantities for which we want z-averages.
do ixy=1,nnamexy
! call parse_name(ixy,cnamexy(ixy),cformxy(ixy),'bxmxy',idiag_bxmxy)
enddo
! Check for those quantities for which we want phi-averages.
do irz=1,nnamerz
! call parse_name(irz,cnamerz(irz),cformrz(irz),'brmphi',idiag_brmphi)
enddo
endsubroutine rprint_cosmicrayflux
!*******************************************************************************
endmodule Cosmicrayflux